ASCOT5
Loading...
Searching...
No Matches
interp3Dexpl.c
Go to the documentation of this file.
1
5#include <stdlib.h>
6#include <math.h>
7#include "../ascot5.h"
8#include "../math.h"
9#include "../consts.h"
10#include "interp.h"
11#include "spline.h"
12
37int interp3Dexpl_init_coeff(real* c, real* f,
38 int n_x, int n_y, int n_z,
39 int bc_x, int bc_y, int bc_z,
40 real x_min, real x_max,
41 real y_min, real y_max,
42 real z_min, real z_max) {
43
44 /* Allocate helper quantities */
45 real* f_x = malloc(n_x*sizeof(real));
46 real* f_y = malloc(n_y*sizeof(real));
47 real* f_z = malloc(n_z*sizeof(real));
48 real* c_x = malloc((n_x-1*(bc_x==NATURALBC))*NSIZE_EXPL1D*sizeof(real));
49 real* c_y = malloc((n_y-1*(bc_y==NATURALBC))*NSIZE_EXPL1D*sizeof(real));
50 real* c_z = malloc((n_z-1*(bc_z==NATURALBC))*NSIZE_EXPL1D*sizeof(real));
51 int i_ct;
52
53 if(f_x == NULL || f_y == NULL || f_z == NULL ||
54 c_x == NULL || c_y == NULL || c_z == NULL) {
55 return 1;
56 }
57
58 /* Calculate tricubic spline volume coefficients. For each grid cell
59 (i_x, i_y, i_z), there are 64 coefficients, one for each variable
60 product dx^p_x*dy^p_y*dz^p_z in the evaluation formula, where p_x, p_y,
61 p_z = 0, 1, 2, 3. Note how we account for normalized grid. */
62
63 /* Bicubic spline surfaces over xy-grid for each z */
64 for(int i_z=0; i_z<n_z; i_z++) {
65
66 /* Cubic spline along x for each y, using f values to get a total of
67 four coefficients */
68 for(int i_y=0; i_y<n_y; i_y++) {
69 for(int i_x=0; i_x<n_x; i_x++) {
70 f_x[i_x] = f[i_z*n_y*n_x+i_y*n_x+i_x];
71 }
72 splineexpl(f_x, n_x, bc_x, c_x);
73 for(int i_x=0; i_x<n_x-1*(bc_x==NATURALBC); i_x++) {
74 for(int i_c=0; i_c<4; i_c++) {
75 c[i_z*n_y*n_x*64+i_y*n_x*64+i_x*64+i_c] = c_x[i_x*4+i_c];
76 }
77 }
78 }
79
80 /* Four cubic splines along y for each x, using the above calulated four
81 coefficient values to get a total of 16 coefficients */
82 for(int i_x=0; i_x<n_x-1*(bc_x==NATURALBC); i_x++) {
83 for(int i_s=0; i_s<4; i_s++) {
84 for(int i_y=0; i_y<n_y; i_y++) {
85 f_y[i_y] = c[i_z*n_y*n_x*64+i_y*n_x*64+i_x*64+i_s];
86 }
87 splineexpl(f_y,n_y,bc_y,c_y);
88 for(int i_y=0; i_y<n_y-1*(bc_y==NATURALBC); i_y++) {
89 i_ct = 0;
90 for(int i_c=i_s; i_c<16; i_c=i_c+4) {
91 c[i_z*n_y*n_x*64+i_y*n_x*64+i_x*64+i_c]
92 = c_y[i_y*4+i_ct];
93 i_ct++;
94 }
95 }
96 }
97 }
98 }
99
100 /* Cubic splines along z for each xy-pair, using the above calculated 16
101 coefficient values to get a total of 64 coefficients */
102 for(int i_y=0; i_y<n_y-1*(bc_y==NATURALBC); i_y++) {
103 for(int i_x=0; i_x<n_x-1*(bc_x==NATURALBC); i_x++) {
104 for(int i_ss=0; i_ss<4; i_ss++) {
105 for(int i_s=0; i_s<4; i_s++) {
106 for(int i_z=0; i_z<n_z; i_z++) {
107 f_z[i_z] = c[i_z*n_y*n_x*64+i_y*n_x*64
108 +i_x*64+(i_ss*4+i_s)];
109 }
110 splineexpl(f_z,n_z,bc_z,c_z);
111 for(int i_z=0; i_z<n_z-1*(bc_z==NATURALBC); i_z++) {
112 i_ct = 0;
113 for(int i_c=4*i_ss+i_s; i_c<64; i_c=i_c+16) {
114 c[i_z*n_y*n_x*64+i_y*n_x*64+i_x*64+i_c]
115 = c_z[i_z*4+i_ct];
116 i_ct++;
117 }
118 }
119 }
120 }
121 }
122 }
123
124 /* Free allocated memory */
125 free(f_x);
126 free(f_y);
127 free(f_z);
128 free(c_x);
129 free(c_y);
130 free(c_z);
131
132 return 0;
133}
134
153void interp3Dexpl_init_spline(interp3D_data* str, real* c,
154 int n_x, int n_y, int n_z,
155 int bc_x, int bc_y, int bc_z,
156 real x_min, real x_max,
157 real y_min, real y_max,
158 real z_min, real z_max) {
159
160 /* Calculate grid intervals. For periodic boundary condition, grid maximum
161 value and the last data point are not the same. Take this into account
162 in grid intervals. */
163 real x_grid = (x_max - x_min) / ( n_x - 1 * (bc_x == NATURALBC) );
164 real y_grid = (y_max - y_min) / ( n_y - 1 * (bc_y == NATURALBC) );
165 real z_grid = (z_max - z_min) / ( n_z - 1 * (bc_z == NATURALBC) );
166
167 /* Initialize the interp3D_data struct */
168 str->n_x = n_x;
169 str->n_y = n_y;
170 str->n_z = n_z;
171 str->bc_x = bc_x;
172 str->bc_y = bc_y;
173 str->bc_z = bc_z;
174 str->x_min = x_min;
175 str->x_max = x_max;
176 str->x_grid = x_grid;
177 str->y_min = y_min;
178 str->y_max = y_max;
179 str->y_grid = y_grid;
180 str->z_min = z_min;
181 str->z_max = z_max;
182 str->z_grid = z_grid;
183 str->c = c;
184}
185
200a5err interp3Dexpl_eval_f(real* f, interp3D_data* str, real x, real y, real z) {
201
202 /* Make sure periodic coordinates are within [min, max] region. */
203 if(str->bc_x == PERIODICBC) {
204 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
205 x = x + (x < str->x_min) * (str->x_max - str->x_min);
206 }
207 if(str->bc_y == PERIODICBC) {
208 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
209 y = y + (y < str->y_min) * (str->y_max - str->y_min);
210 }
211 if(str->bc_z == PERIODICBC) {
212 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
213 z = z + (z < str->z_min) * (str->z_max - str->z_min);
214 }
215
216 /* Index for x variable. The -1 needed at exactly grid end. */
217 int i_x = (x-str->x_min)/str->x_grid - 1*(x==str->x_max);
218 /* Normalized x coordinate in current cell */
219 real dx = (x-(str->x_min+i_x*str->x_grid))/str->x_grid;
220 /* Helper varibles */
221 real dx2 = dx*dx;
222 real dx3 = dx2*dx;
223
224 /* Index for y variable. The -1 needed at exactly grid end. */
225 int i_y = (y-str->y_min)/str->y_grid - 1*(y==str->y_max);
226 /* Normalized y coordinate in current cell */
227 real dy = (y-(str->y_min+i_y*str->y_grid))/str->y_grid;
228 /* Helper varibles */
229 real dy2 = dy*dy;
230 real dy3 = dy2*dy;
231
232 /* Index for z variable. The -1 needed at exactly grid end. */
233 int i_z = (z-str->z_min)/str->z_grid - 1*(z==str->z_max);
234 /* Normalized z coordinate in current cell */
235 real dz = (z-(str->z_min+i_z*str->z_grid))/str->z_grid;
236 /* Helper varibles */
237 real dz2 = dz*dz;
238 real dz3 = dz2*dz;
239
241 int n = i_z*str->n_y*str->n_x*64+i_y*str->n_x*64+i_x*64;
242
243 int err = 0;
244
245 /* Check that the coordinate is within the domain. */
246 if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
247 err = 1;
248 }
249 if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
250 err = 1;
251 }
252 if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
253 err = 1;
254 }
255
256 if(!err) {
257
258 /* Evaluate spline value */
259 *f = (
260 str->c[n+ 0]+dx*str->c[n+ 1]
261 +dx2*str->c[n+ 2]+dx3*str->c[n+ 3]
262 +dy*(
263 str->c[n+ 4]+dx*str->c[n+ 5]
264 +dx2*str->c[n+ 6]+dx3*str->c[n+ 7])
265 +dy2*(
266 str->c[n+ 8]+dx*str->c[n+ 9]
267 +dx2*str->c[n+10]+dx3*str->c[n+11])
268 +dy3*(
269 str->c[n+12]+dx*str->c[n+13]
270 +dx2*str->c[n+14]+dx3*str->c[n+15]))
271 +dz*(
272 str->c[n+16]+dx*str->c[n+17]
273 +dx2*str->c[n+18]+dx3*str->c[n+19]
274 +dy*(
275 str->c[n+20]+dx*str->c[n+21]
276 +dx2*str->c[n+22]+dx3*str->c[n+23])
277 +dy2*(
278 str->c[n+24]+dx*str->c[n+25]
279 +dx2*str->c[n+26]+dx3*str->c[n+27])
280 +dy3*(
281 str->c[n+28]+dx*str->c[n+29]
282 +dx2*str->c[n+30]+dx3*str->c[n+31]))
283 +dz2*(
284 str->c[n+32]+dx*str->c[n+33]
285 +dx2*str->c[n+34]+dx3*str->c[n+35]
286 +dy*(
287 str->c[n+36]+dx*str->c[n+37]
288 +dx2*str->c[n+38]+dx3*str->c[n+39])
289 +dy2*(
290 str->c[n+40]+dx*str->c[n+41]
291 +dx2*str->c[n+42]+dx3*str->c[n+43])
292 +dy3*(
293 str->c[n+44]+dx*str->c[n+45]
294 +dx2*str->c[n+46]+dx3*str->c[n+47]))
295 +dz3*(
296 str->c[n+48]+dx*str->c[n+49]
297 +dx2*str->c[n+50]+dx3*str->c[n+51]
298 +dy*(
299 str->c[n+52]+dx*str->c[n+53]
300 +dx2*str->c[n+54]+dx3*str->c[n+55])
301 +dy2*(
302 str->c[n+56]+dx*str->c[n+57]
303 +dx2*str->c[n+58]+dx3*str->c[n+59])
304 +dy3*(
305 str->c[n+60]+dx*str->c[n+61]
306 +dx2*str->c[n+62]+dx3*str->c[n+63]));
307
308 }
309
310 return err;
311}
312
340a5err interp3Dexpl_eval_df(real* f_df, interp3D_data* str, real x, real y, real z) {
341
342 /* Make sure periodic coordinates are within [min, max] region. */
343 if(str->bc_x == PERIODICBC) {
344 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
345 x = x + (x < str->x_min) * (str->x_max - str->x_min);
346 }
347 if(str->bc_y == PERIODICBC) {
348 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
349 y = y + (y < str->y_min) * (str->y_max - str->y_min);
350 }
351 if(str->bc_z == PERIODICBC) {
352 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
353 z = z + (z < str->z_min) * (str->z_max - str->z_min);
354 }
355
356 /* Index for x variable. The -1 needed at exactly grid end. */
357 int i_x = (x-str->x_min)/str->x_grid - 1*(x==str->x_max);
358 /* Normalized x coordinate in current cell */
359 real dx = (x-(str->x_min+i_x*str->x_grid))/str->x_grid;
360 /* Helper variables */
361 real dx2 = dx*dx;
362 real dx3 = dx2*dx;
363 real xgi = 1.0/str->x_grid;
364
365 /* Index for y variable. The -1 needed at exactly grid end. */
366 int i_y = (y-str->y_min)/str->y_grid - 1*(y==str->y_max);
367 /* Normalized y coordinate in current cell */
368 real dy = (y-(str->y_min+i_y*str->y_grid))/str->y_grid;
369 /* Helper variables */
370 real dy2 = dy*dy;
371 real dy3 = dy2*dy;
372 real ygi = 1.0/str->y_grid;
373
374 /* Index for z variable. The -1 needed at exactly grid end. */
375 int i_z = (z-str->z_min)/str->z_grid - 1*(z==str->z_max);
376 /* Normalized z coordinate in current cell */
377 real dz = (z-(str->z_min+i_z*str->z_grid))/str->z_grid;
378 /* Helper variables */
379 real dz2 = dz*dz;
380 real dz3 = dz2*dz;
381 real zgi = 1.0/str->z_grid;
382
383 /* Index jump to cell */
384 int n = i_z*str->n_y*str->n_x*64+i_y*str->n_x*64+i_x*64;
385
386 int err = 0;
387
388 /* Check that the coordinate is within the domain. */
389 if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
390 err = 1;
391 }
392 if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
393 err = 1;
394 }
395 if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
396 err = 1;
397 }
398
399 if(!err) {
400
401 /* Fetch coefficients explicitly to fetch those that are adjacent
402 subsequently and to store in temporary variables coefficients that
403 will be used multiple times. This is to decrease computational time,
404 by exploiting simultaneous extraction of adjacent memory and by
405 avoiding going through the long str->c array repeatedly. */
406 real c000 = str->c[n+0];
407 real c001 = str->c[n+1];
408 real c002 = str->c[n+2];
409 real c003 = str->c[n+3];
410 real c010 = str->c[n+4];
411 real c011 = str->c[n+5];
412 real c012 = str->c[n+6];
413 real c013 = str->c[n+7];
414
415 real c020 = str->c[n+8];
416 real c021 = str->c[n+9];
417 real c022 = str->c[n+10];
418 real c023 = str->c[n+11];
419 real c030 = str->c[n+12];
420 real c031 = str->c[n+13];
421 real c032 = str->c[n+14];
422 real c033 = str->c[n+15];
423
424 real c100 = str->c[n+16];
425 real c101 = str->c[n+17];
426 real c102 = str->c[n+18];
427 real c103 = str->c[n+19];
428 real c110 = str->c[n+20];
429 real c111 = str->c[n+21];
430 real c112 = str->c[n+22];
431 real c113 = str->c[n+23];
432
433 real c120 = str->c[n+24];
434 real c121 = str->c[n+25];
435 real c122 = str->c[n+26];
436 real c123 = str->c[n+27];
437 real c130 = str->c[n+28];
438 real c131 = str->c[n+29];
439 real c132 = str->c[n+30];
440 real c133 = str->c[n+31];
441
442 real c200 = str->c[n+32];
443 real c201 = str->c[n+33];
444 real c202 = str->c[n+34];
445 real c203 = str->c[n+35];
446 real c210 = str->c[n+36];
447 real c211 = str->c[n+37];
448 real c212 = str->c[n+38];
449 real c213 = str->c[n+39];
450
451 real c220 = str->c[n+40];
452 real c221 = str->c[n+41];
453 real c222 = str->c[n+42];
454 real c223 = str->c[n+43];
455 real c230 = str->c[n+44];
456 real c231 = str->c[n+45];
457 real c232 = str->c[n+46];
458 real c233 = str->c[n+47];
459
460 real c300 = str->c[n+48];
461 real c301 = str->c[n+49];
462 real c302 = str->c[n+50];
463 real c303 = str->c[n+51];
464 real c310 = str->c[n+52];
465 real c311 = str->c[n+53];
466 real c312 = str->c[n+54];
467 real c313 = str->c[n+55];
468
469 real c320 = str->c[n+56];
470 real c321 = str->c[n+57];
471 real c322 = str->c[n+58];
472 real c323 = str->c[n+59];
473 real c330 = str->c[n+60];
474 real c331 = str->c[n+61];
475 real c332 = str->c[n+62];
476 real c333 = str->c[n+63];
477
478 /* Evaluate spline value */
479
480 /* f */
481 f_df[0] = (
482 c000+dx*c001+dx2*c002+dx3*c003
483 +dy*(c010+dx*c011+dx2*c012+dx3*c013)
484 +dy2*(c020+dx*c021+dx2*c022+dx3*c023)
485 +dy3*(c030+dx*c031+dx2*c032+dx3*c033))
486 +dz*(
487 c100+dx*c101+dx2*c102+dx3*c103
488 +dy*(c110+dx*c111+dx2*c112+dx3*c113)
489 +dy2*(c120+dx*c121+dx2*c122+dx3*c123)
490 +dy3*(c130+dx*c131+dx2*c132+dx3*c133))
491 +dz2*(
492 c200+dx*c201+dx2*c202+dx3*c203
493 +dy*(c210+dx*c211+dx2*c212+dx3*c213)
494 +dy2*(c220+dx*c221+dx2*c222+dx3*c223)
495 +dy3*(c230+dx*c231+dx2*c232+dx3*c233))
496 +dz3*(
497 c300+dx*c301+dx2*c302+dx3*c303
498 +dy*(c310+dx*c311+dx2*c312+dx3*c313)
499 +dy2*(c320+dx*c321+dx2*c322+dx3*c323)
500 +dy3*(c330+dx*c331+dx2*c332+dx3*c333));
501
502 /* df/dx */
503 f_df[1] = xgi*(
504 (
505 c001+2*dx*c002+3*dx2*c003
506 +dy*(c011+2*dx*c012+3*dx2*c013)
507 +dy2*(c021+2*dx*c022+3*dx2*c023)
508 +dy3*(c031+2*dx*c032+3*dx2*c033))
509 +dz*(
510 c101+2*dx*c102+3*dx2*c103
511 +dy*(c111+2*dx*c112+3*dx2*c113)
512 +dy2*(c121+2*dx*c122+3*dx2*c123)
513 +dy3*(c131+2*dx*c132+3*dx2*c133))
514 +dz2*(
515 c201+2*dx*c202+3*dx2*c203
516 +dy*(c211+2*dx*c212+3*dx2*c213)
517 +dy2*(c221+2*dx*c222+3*dx2*c223)
518 +dy3*(c231+2*dx*c232+3*dx2*c233))
519 +dz3*(
520 c301+2*dx*c302+3*dx2*c303
521 +dy*(c311+2*dx*c312+3*dx2*c313)
522 +dy2*(c321+2*dx*c322+3*dx2*c323)
523 +dy3*(c331+2*dx*c332+3*dx2*c333)));
524
525 /* df/dy */
526 f_df[2] = ygi*(
527 (
528 c010+dx*c011+dx2*c012+dx3*c013
529 +2*dy*(c020+dx*c021+dx2*c022+dx3*c023)
530 +3*dy2*(c030+dx*c031+dx2*c032+dx3*c033))
531 +dz*(
532 c110+dx*c111+dx2*c112+dx3*c113
533 +2*dy*(c120+dx*c121+dx2*c122+dx3*c123)
534 +3*dy2*(c130+dx*c131+dx2*c132+dx3*c133))
535 +dz2*(
536 c210+dx*c211+dx2*c212+dx3*c213
537 +2*dy*(c220+dx*c221+dx2*c222+dx3*c223)
538 +3*dy2*(c230+dx*c231+dx2*c232+dx3*c233))
539 +dz3*(
540 c310+dx*c311+dx2*c312+dx3*c313
541 +2*dy*(c320+dx*c321+dx2*c322+dx3*c323)
542 +3*dy2*(c330+dx*c331+dx2*c332+dx3*c333)));
543
544 /* df/dz */
545 f_df[3] = zgi*(
546 (
547 c100+dx*c101+dx2*c102+dx3*c103
548 +dy*(c110+dx*c111+dx2*c112+dx3*c113)
549 +dy2*(c120+dx*c121+dx2*c122+dx3*c123)
550 +dy3*(c130+dx*c131+dx2*c132+dx3*c133))
551 +2*dz*(
552 c200+dx*c201+dx2*c202+dx3*c203
553 +dy*(c210+dx*c211+dx2*c212+dx3*c213)
554 +dy2*(c220+dx*c221+dx2*c222+dx3*c223)
555 +dy3*(c230+dx*c231+dx2*c232+dx3*c233))
556 +3*dz2*(
557 c300+dx*c301+dx2*c302+dx3*c303
558 +dy*(c310+dx*c311+dx2*c312+dx3*c313)
559 +dy2*(c320+dx*c321+dx2*c322+dx3*c323)
560 +dy3*(c330+dx*c331+dx2*c332+dx3*c333)));
561
562 /* d2f/dx2 */
563 f_df[4] = xgi*xgi*(
564 (
565 2*c002+6*dx*c003
566 +dy*(2*c012+6*dx*c013)
567 +dy2*(2*c022+6*dx*c023)
568 +dy3*(2*c032+6*dx*c033))
569 +dz*(
570 2*c102+6*dx*c103
571 +dy*(2*c112+6*dx*c113)
572 +dy2*(2*c122+6*dx*c123)
573 +dy3*(2*c132+6*dx*c133))
574 +dz2*(
575 2*c202+6*dx*c203
576 +dy*(2*c212+6*dx*c213)
577 +dy2*(2*c222+6*dx*c223)
578 +dy3*(2*c232+6*dx*c233))
579 +dz3*(
580 2*c302+6*dx*c303
581 +dy*(2*c312+6*dx*c313)
582 +dy2*(2*c322+6*dx*c323)
583 +dy3*(2*c332+6*dx*c333)));
584
585 /* d2f/dy2 */
586 f_df[5] = ygi*ygi*(
587 (
588 2*(c020+dx*c021+dx2*c022+dx3*c023)
589 +6*dy*(c030+dx*c031+dx2*c032+dx3*c033))
590 +dz*(
591 2*(c120+dx*c121+dx2*c122+dx3*c123)
592 +6*dy*(c130+dx*c131+dx2*c132+dx3*c133))
593 +dz2*(
594 2*(c220+dx*c221+dx2*c222+dx3*c223)
595 +6*dy*(c230+dx*c231+dx2*c232+dx3*c233))
596 +dz3*(
597 2*(c320+dx*c321+dx2*c322+dx3*c323)
598 +6*dy*(c330+dx*c331+dx2*c332+dx3*c333)));
599
600 /* d2f/dz2 */
601 f_df[6] = zgi*zgi*(
602 2*(
603 c200+dx*c201+dx2*c202+dx3*c203
604 +dy*(c210+dx*c211+dx2*c212+dx3*c213)
605 +dy2*(c220+dx*c221+dx2*c222+dx3*c223)
606 +dy3*(c230+dx*c231+dx2*c232+dx3*c233))
607 +6*dz*(
608 c300+dx*c301+dx2*c302+dx3*c303
609 +dy*(c310+dx*c311+dx2*c312+dx3*c313)
610 +dy2*(c320+dx*c321+dx2*c322+dx3*c323)
611 +dy3*(c330+dx*c331+dx2*c332+dx3*c333)));
612
613 /* d2f/dydx */
614 f_df[7] = xgi*ygi*(
615 (
616 c011+2*dx*c012+3*dx2*c013
617 +2*dy*(c021+2*dx*c022+3*dx2*c023)
618 +3*dy2*(c031+2*dx*c032+3*dx2*c033))
619 +dz*(
620 c111+2*dx*c112+3*dx2*c113
621 +2*dy*(c121+2*dx*c122+3*dx2*c123)
622 +3*dy2*(c131+2*dx*c132+3*dx2*c133))
623 +dz2*(
624 c211+2*dx*c212+3*dx2*c213
625 +2*dy*(c221+2*dx*c222+3*dx2*c223)
626 +3*dy2*(c231+2*dx*c232+3*dx2*c233))
627 +dz3*(
628 c311+2*dx*c312+3*dx2*c313
629 +2*dy*(c321+2*dx*c322+3*dx2*c323)
630 +3*dy2*(c331+2*dx*c332+3*dx2*c333)));
631
632 /* d2f/dzdx */
633 f_df[8] = xgi*zgi*(
634 (
635 c101+2*dx*c102+3*dx2*c103
636 +dy*(c111+2*dx*c112+3*dx2*c113)
637 +dy2*(c121+2*dx*c122+3*dx2*c123)
638 +dy3*(c131+2*dx*c132+3*dx2*c133))
639 +2*dz*(
640 c201+2*dx*c202+3*dx2*c203
641 +dy*(c211+2*dx*c212+3*dx2*c213)
642 +dy2*(c221+2*dx*c222+3*dx2*c223)
643 +dy3*(c231+2*dx*c232+3*dx2*c233))
644 +3*dz2*(
645 c301+2*dx*c302+3*dx2*c303
646 +dy*(c311+2*dx*c312+3*dx2*c313)
647 +dy2*(c321+2*dx*c322+3*dx2*c323)
648 +dy3*(c331+2*dx*c332+3*dx2*c333)));
649
650 /* d2f/dzdy */
651 f_df[9] = ygi*zgi*(
652 (
653 c110+dx*c111+dx2*c112+dx3*c113
654 +2*dy*(c120+dx*c121+dx2*c122+dx3*c123)
655 +3*dy2*(c130+dx*c131+dx2*c132+dx3*c133))
656 +2*dz*(
657 c210+dx*c211+dx2*c212+dx3*c213
658 +2*dy*(c220+dx*c221+dx2*c222+dx3*c223)
659 +3*dy2*(c230+dx*c231+dx2*c232+dx3*c233))
660 +3*dz2*(
661 c310+dx*c311+dx2*c312+dx3*c313
662 +2*dy*(c320+dx*c321+dx2*c322+dx3*c323)
663 +3*dy2*(c330+dx*c331+dx2*c332+dx3*c333)));
664 }
665
666 return err;
667}
ascot5.h
Main header file for ASCOT5.
real
double real
Definition ascot5.h:85
consts.h
Header file containing physical and mathematical constants.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
interp3Dexpl_init_coeff
int interp3Dexpl_init_coeff(real *c, real *f, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Calculate tricubic spline interpolation coefficients for 3D data.
Definition interp3Dexpl.c:37
interp3Dexpl_eval_df
a5err interp3Dexpl_eval_df(real *f_df, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D field and 1st and 2nd derivatives.
Definition interp3Dexpl.c:340
interp3Dexpl_eval_f
a5err interp3Dexpl_eval_f(real *f, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D scalar field.
Definition interp3Dexpl.c:200
interp3Dexpl_init_spline
void interp3Dexpl_init_spline(interp3D_data *str, real *c, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Initialize a tricubic spline.
Definition interp3Dexpl.c:153
interp.h
Spline interpolation library.
NATURALBC
@ NATURALBC
Definition interp.h:37
PERIODICBC
@ PERIODICBC
Definition interp.h:38
fmod
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
math.h
Header file for math.c.
spline.h
Header file for splineexpl.c and splinecomp.c.
splineexpl
void splineexpl(real *f, int n, int bc, real *c)
Calculate explicit cubic spline interpolation coefficients in 1D.
Definition splineexpl.c:22
interp3D_data
Tricubic interpolation struct.
Definition interp.h:85
interp3D_data::bc_x
int bc_x
Definition interp.h:89
interp3D_data::bc_z
int bc_z
Definition interp.h:91
interp3D_data::z_max
real z_max
Definition interp.h:99
interp3D_data::z_min
real z_min
Definition interp.h:98
interp3D_data::z_grid
real z_grid
Definition interp.h:100
interp3D_data::x_min
real x_min
Definition interp.h:92
interp3D_data::y_max
real y_max
Definition interp.h:96
interp3D_data::y_grid
real y_grid
Definition interp.h:97
interp3D_data::bc_y
int bc_y
Definition interp.h:90
interp3D_data::c
real * c
Definition interp.h:101
interp3D_data::n_y
int n_y
Definition interp.h:87
interp3D_data::n_z
int n_z
Definition interp.h:88
interp3D_data::x_grid
real x_grid
Definition interp.h:94
interp3D_data::x_max
real x_max
Definition interp.h:93
interp3D_data::y_min
real y_min
Definition interp.h:95
interp3D_data::n_x
int n_x
Definition interp.h:86
Generated on Wed Oct 9 2024 17:13:57 for ASCOT5 by doxygen 1.10.0