ASCOT5
Loading...
Searching...
No Matches
interp3Dcomp.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
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 /* Check boundary conditions and calculate grid intervals. Grid intervals
45 needed because we use normalized grid intervals. For periodic boundary
46 condition, grid maximum value and the last data point are not the same.
47 Take this into account in grid intervals. */
48 real x_grid, y_grid, z_grid;
49 if(bc_x == NATURALBC || bc_x == PERIODICBC) {
50 x_grid = (x_max - x_min) / ( n_x - 1 * (bc_x == NATURALBC) );
51 }
52 else {
53 return 1;
54 }
55
56 if(bc_y == NATURALBC || bc_y == PERIODICBC) {
57 y_grid = (y_max - y_min) / ( n_y - 1 * (bc_y == NATURALBC) );
58 }
59 else {
60 return 1;
61 }
62
63 if(bc_z == NATURALBC || bc_z == PERIODICBC) {
64 z_grid = (z_max - z_min) / ( n_z - 1 * (bc_z == NATURALBC) );
65 }
66 else {
67 return 1;
68 }
69
70 /* Allocate helper quantities */
71 real* f_x = malloc(n_x*sizeof(real));
72 real* f_y = malloc(n_y*sizeof(real));
73 real* f_z = malloc(n_z*sizeof(real));
74 real* c_x = malloc(n_x*NSIZE_COMP1D*sizeof(real));
75 real* c_y = malloc(n_y*NSIZE_COMP1D*sizeof(real));
76 real* c_z = malloc(n_z*NSIZE_COMP1D*sizeof(real));
77
78 if(f_x == NULL || f_y == NULL || f_z == NULL ||
79 c_x == NULL || c_y == NULL || c_z == NULL) {
80 return 1;
81 }
82
83 /* Calculate tricubic spline volume coefficients, i.e. second derivatives.
84 For each grid cell (i_x, i_y, i_z), there are eight coefficients:
85 [f, fxx, fyy, fzz, fxxyy, fxxzz, fyyzz, fxxyyzz]. Note how we account
86 for normalized grid intervals. */
87
88 /* Bicubic spline surfaces over xy-grid for each z */
89 for(int i_z=0; i_z<n_z; i_z++) {
90
91 /* Cubic spline along x for each y, using f values to get fxx */
92 for(int i_y=0; i_y<n_y; i_y++) {
93 /* fxx */
94 for(int i_x=0; i_x<n_x; i_x++) {
95 f_x[i_x] = f[i_z*n_y*n_x + i_y*n_x + i_x];
96 }
97 splinecomp(f_x, n_x, bc_x, c_x);
98 for(int i_x=0; i_x<n_x; i_x++) {
99 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 ] = c_x[i_x*2];
100 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 1] = c_x[i_x*2 + 1]
101 / (x_grid*x_grid);
102 }
103 }
104
105 /* Two cubic splines along y for each x, one using f values to
106 get fyy, and the other using fxx values to get fxxyy */
107 for(int i_x=0; i_x<n_x; i_x++) {
108 /* fyy */
109 for(int i_y=0; i_y<n_y; i_y++) {
110 f_y[i_y] = f[i_z*n_y*n_x + i_y*n_x + i_x];
111 }
112 splinecomp(f_y, n_y, bc_y, c_y);
113 for(int i_y=0; i_y<n_y; i_y++) {
114 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 2] = c_y[i_y*2 + 1]
115 / (y_grid*y_grid);
116 }
117 /* fxxyy */
118 for(int i_y=0; i_y<n_y; i_y++) {
119 f_y[i_y] = c[i_z*n_y*n_x*8 + i_y*n_x*8+i_x*8 + 1];
120 }
121 splinecomp(f_y, n_y, bc_y, c_y);
122 for(int i_y=0; i_y<n_y; i_y++) {
123 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 4] = c_y[i_y*2 + 1]
124 / (y_grid*y_grid);
125 }
126 }
127
128 }
129
130 /* Four cubic splines along z for each xy-pair, one using f values to get
131 fzz, one using fxx to get fxxzz, one using fyy to get fyyzz, and one
132 using fxxyy to get fxxyyzz */
133 for(int i_y=0; i_y<n_y; i_y++) {
134 for(int i_x=0; i_x<n_x; i_x++) {
135 /* fzz */
136 for(int i_z=0; i_z<n_z; i_z++) {
137 f_z[i_z] = f[i_z*n_y*n_x + i_y*n_x + i_x];
138 }
139 splinecomp(f_z, n_z, bc_z, c_z);
140 for(int i_z=0; i_z<n_z; i_z++) {
141 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 3] = c_z[i_z*2 + 1]
142 / (z_grid*z_grid);
143 }
144 /* fxxzz */
145 for(int i_z=0; i_z<n_z; i_z++) {
146 f_z[i_z] = c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 1];
147 }
148 splinecomp(f_z, n_z, bc_z, c_z);
149 for(int i_z=0; i_z<n_z; i_z++) {
150 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 5] = c_z[i_z*2 + 1]
151 / (z_grid*z_grid);
152 }
153 /* fyyzz */
154 for(int i_z=0; i_z<n_z; i_z++) {
155 f_z[i_z] = c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 2];
156 }
157 splinecomp(f_z, n_z, bc_z, c_z);
158 for(int i_z=0; i_z<n_z; i_z++) {
159 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 6] = c_z[i_z*2+1]
160 / (z_grid*z_grid);
161 }
162 /* fxxyyzz */
163 for(int i_z=0; i_z<n_z; i_z++) {
164 f_z[i_z] = c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 4];
165 }
166 splinecomp(f_z, n_z, bc_z, c_z);
167 for(int i_z=0; i_z<n_z; i_z++) {
168 c[i_z*n_y*n_x*8 + i_y*n_x*8 + i_x*8 + 7] = c_z[i_z*2+1]
169 / (z_grid*z_grid);
170 }
171 }
172 }
173
174 /* Free allocated memory */
175 free(f_x);
176 free(f_y);
177 free(f_z);
178 free(c_x);
179 free(c_y);
180 free(c_z);
181
182 return 0;
183}
184
204 int n_x, int n_y, int n_z,
205 int bc_x, int bc_y, int bc_z,
206 real x_min, real x_max,
207 real y_min, real y_max,
208 real z_min, real z_max) {
209
210 /* Calculate grid intervals. For periodic boundary condition, grid maximum
211 value and the last data point are not the same. Take this into account
212 in grid intervals. */
213 real x_grid = (x_max - x_min) / ( n_x - 1 * (bc_x == NATURALBC) );
214 real y_grid = (y_max - y_min) / ( n_y - 1 * (bc_y == NATURALBC) );
215 real z_grid = (z_max - z_min) / ( n_z - 1 * (bc_z == NATURALBC) );
216
217 /* Initialize the interp3D_data struct */
218 str->n_x = n_x;
219 str->n_y = n_y;
220 str->n_z = n_z;
221 str->bc_x = bc_x;
222 str->bc_y = bc_y;
223 str->bc_z = bc_z;
224 str->x_min = x_min;
225 str->x_max = x_max;
226 str->x_grid = x_grid;
227 str->y_min = y_min;
228 str->y_max = y_max;
229 str->y_grid = y_grid;
230 str->z_min = z_min;
231 str->z_max = z_max;
232 str->z_grid = z_grid;
233 str->c = c;
234}
235
260 int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z,
261 real x_min, real x_max, real y_min, real y_max,
262 real z_min, real z_max) {
263 real* c = (real*) malloc(n_z*n_y*n_x*NSIZE_COMP3D*sizeof(real));
264 int err = interp3Dcomp_init_coeff(c, f, n_x, n_y, n_z, bc_x, bc_y, bc_z,
265 x_min, x_max, y_min, y_max, z_min, z_max);
266 if(err) {
267 return err;
268 }
269 interp3Dcomp_init_spline(str, c, n_x, n_y, n_z, bc_x, bc_y, bc_z,
270 x_min, x_max, y_min, y_max, z_min, z_max);
271 return 0;
272}
273
289
290 /* Make sure periodic coordinates are within [min, max] region. */
291 if(str->bc_x == PERIODICBC) {
292 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
293 x = x + (x < str->x_min) * (str->x_max - str->x_min);
294 }
295 if(str->bc_y == PERIODICBC) {
296 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
297 y = y + (y < str->y_min) * (str->y_max - str->y_min);
298 }
299 if(str->bc_z == PERIODICBC) {
300 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
301 z = z + (z < str->z_min) * (str->z_max - str->z_min);
302 }
303
304 /* Index for x variable. The -1 needed at exactly grid end. */
305 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
306 /* Normalized x coordinate in current cell */
307 real dx = (x - (str->x_min + i_x*str->x_grid)) / str->x_grid;
308 /* Helper varibles */
309 real dxi = 1.0 - dx;
310 real dx3 = dx*dx*dx - dx;
311 real dxi3 = (1.0 - dx) * (1.0 - dx) * (1.0 - dx) - (1.0 - dx);
312 real xg2 = str->x_grid*str->x_grid;
313
314 /* Index for y variable. The -1 needed at exactly grid end. */
315 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
316 /* Normalized y coordinate in current cell */
317 real dy = (y - (str->y_min + i_y*str->y_grid)) / str->y_grid;
318 /* Helper varibles */
319 real dyi = 1.0 - dy;
320 real dy3 = dy*dy*dy - dy;
321 real dyi3 = (1.0 - dy) * (1.0 - dy) * (1.0 - dy) - (1.0 - dy);
322 real yg2 = str->y_grid*str->y_grid;
323
324 /* Index for z variable. The -1 needed at exactly grid end. */
325 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
326 /* Normalized z coordinate in current cell */
327 real dz = (z - (str->z_min + i_z*str->z_grid)) / str->z_grid;
328 /* Helper varibles */
329 real dzi = 1.0 - dz;
330 real dz3 = dz*dz*dz - dz;
331 real dzi3 = (1.0 - dz) * (1.0 - dz) * (1.0 - dz) - (1.0-dz);
332 real zg2 = str->z_grid*str->z_grid;
333
335 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
336 int x1 = 8; /* Index jump one x forward */
337 int y1 = str->n_x*8; /* Index jump one y forward */
338 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
339
340 int err = 0;
341
342 /* Enforce periodic BC or check that the coordinate is within the domain. */
343 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
344 x1 = -(str->n_x-1)*x1;
345 }
346 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
347 err = 1;
348 }
349 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
350 y1 = -(str->n_y-1)*y1;
351 }
352 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
353 err = 1;
354 }
355 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
356 z1 = -(str->n_z-1)*z1;
357 }
358 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
359 err = 1;
360 }
361
362 if(!err) {
363
364 /* Evaluate spline value */
365 *f = (
366 dzi*(
367 dxi*(dyi*str->c[n+0]+dy*str->c[n+y1+0])
368 +dx*(dyi*str->c[n+x1+0]+dy*str->c[n+y1+x1+0]))
369 +dz*(
370 dxi*(dyi*str->c[n+z1+0]+dy*str->c[n+y1+z1+0])
371 +dx*(dyi*str->c[n+x1+z1+0]+dy*str->c[n+y1+z1+x1+0])))
372 +xg2/6*(
373 dzi*(
374 dxi3*(dyi*str->c[n+1]+dy*str->c[n+y1+1])
375 +dx3*(dyi*str->c[n+x1+1]+dy*str->c[n+y1+x1+1]))
376 +dz*(
377 dxi3*(dyi*str->c[n+z1+1]+dy*str->c[n+y1+z1+1])
378 +dx3*(dyi*str->c[n+x1+z1+1]+dy*str->c[n+y1+z1+x1+1])))
379 +yg2/6*(
380 dzi*(
381 dxi*(dyi3*str->c[n+2]+dy3*str->c[n+y1+2])
382 +dx*(dyi3*str->c[n+x1+2]+dy3*str->c[n+y1+x1+2]))
383 +dz*(
384 dxi*(dyi3*str->c[n+z1+2]+dy3*str->c[n+y1+z1+2])
385 +dx*(dyi3*str->c[n+x1+z1+2]+dy3*str->c[n+y1+z1+x1+2])))
386 +zg2/6*(
387 dzi3*(
388 dxi*(dyi*str->c[n+3]+dy*str->c[n+y1+3])
389 +dx*(dyi*str->c[n+x1+3]+dy*str->c[n+y1+x1+3]))
390 +dz3*(
391 dxi*(dyi*str->c[n+z1+3]+dy*str->c[n+y1+z1+3])
392 +dx*(dyi*str->c[n+x1+z1+3]+dy*str->c[n+y1+z1+x1+3])))
393 +xg2*yg2/36*(
394 dzi*(
395 dxi3*(dyi3*str->c[n+4]+dy3*str->c[n+y1+4])
396 +dx3*(dyi3*str->c[n+x1+4]+dy3*str->c[n+y1+x1+4]))
397 +dz*(
398 dxi3*(dyi3*str->c[n+z1+4]+dy3*str->c[n+y1+z1+4])
399 +dx3*(dyi3*str->c[n+x1+z1+4]+dy3*str->c[n+y1+z1+x1+4])))
400 +xg2*zg2/36*(
401 dzi3*(
402 dxi3*(dyi*str->c[n+5]+dy*str->c[n+y1+5])
403 +dx3*(dyi*str->c[n+x1+5]+dy*str->c[n+y1+x1+5]))
404 +dz3*(
405 dxi3*(dyi*str->c[n+z1+5]+dy*str->c[n+y1+z1+5])
406 +dx3*(dyi*str->c[n+x1+z1+5]+dy*str->c[n+y1+z1+x1+5])))
407 +yg2*zg2/36*(
408 dzi3*(
409 dxi*(dyi3*str->c[n+6]+dy3*str->c[n+y1+6])
410 +dx*(dyi3*str->c[n+x1+6]+dy3*str->c[n+y1+x1+6]))
411 +dz3*(
412 dxi*(dyi3*str->c[n+z1+6]+dy3*str->c[n+y1+z1+6])
413 +dx*(dyi3*str->c[n+x1+z1+6]+dy3*str->c[n+y1+z1+x1+6])))
414 +xg2*yg2*zg2/216*(
415 dzi3*(
416 dxi3*(dyi3*str->c[n+7]+dy3*str->c[n+y1+7])
417 +dx3*(dyi3*str->c[n+x1+7]+dy3*str->c[n+y1+x1+7]))
418 +dz3*(
419 dxi3*(dyi3*str->c[n+z1+7]+dy3*str->c[n+y1+z1+7])
420 +dx3*(dyi3*str->c[n+x1+z1+7]+dy3*str->c[n+y1+z1+x1+7])));
421
422 }
423
424 return err;
425}
426
455 real x, real y, real z) {
456
457 /* Make sure periodic coordinates are within [min, max] region. */
458 if(str->bc_x == PERIODICBC) {
459 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
460 x = x + (x < str->x_min) * (str->x_max - str->x_min);
461 }
462 if(str->bc_y == PERIODICBC) {
463 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
464 y = y + (y < str->y_min) * (str->y_max - str->y_min);
465 }
466 if(str->bc_z == PERIODICBC) {
467 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
468 z = z + (z < str->z_min) * (str->z_max - str->z_min);
469 }
470
471 /* Index for x variable. The -1 needed at exactly grid end. */
472 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
473 /* Normalized x coordinate in current cell */
474 real dx = ( x - (str->x_min + i_x*str->x_grid) ) / str->x_grid;
475 /* Helper variables */
476 real dx3 = dx*dx*dx - dx;
477 real dx3dx = 3*dx*dx - 1.0;
478 real dxi = 1.0 - dx;
479 real dxi3 = dxi*dxi*dxi - dxi;
480 real dxi3dx = -3*dxi*dxi + 1.0;
481 real xg = str->x_grid;
482 real xg2 = xg*xg;
483 real xgi = 1.0 / xg;
484
485 /* Index for y variable. The -1 needed at exactly grid end. */
486 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
487 /* Normalized y coordinate in current cell */
488 real dy = ( y - (str->y_min + i_y*str->y_grid) ) / str->y_grid;
489 /* Helper variables */
490 real dy3 = dy*dy*dy-dy;
491 real dy3dy = 3*dy*dy - 1.0;
492 real dyi = 1.0 - dy;
493 real dyi3 = dyi*dyi*dyi - dyi;
494 real dyi3dy = -3*dyi*dyi + 1.0;
495 real yg = str->y_grid;
496 real yg2 = yg*yg;
497 real ygi = 1.0 / yg;
498
499 /* Index for z variable. The -1 needed at exactly grid end. */
500 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
501 /* Normalized z coordinate in current cell */
502 real dz = ( z - (str->z_min + i_z*str->z_grid) ) / str->z_grid;
503 /* Helper variables */
504 real dz3 = dz*dz*dz - dz;
505 real dz3dz = 3*dz*dz - 1.0;
506 real dzi = 1.0 - dz;
507 real dzi3 = dzi*dzi*dzi - dzi;
508 real dzi3dz = -3*dzi*dzi + 1.0;
509 real zg = str->z_grid;
510 real zg2 = zg*zg;
511 real zgi = 1.0 / zg;
512
513 /* Index jump to cell */
514 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
515 int x1 = 8; /* Index jump one x forward */
516 int y1 = str->n_x*8; /* Index jump one y forward */
517 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
518
519 int err = 0;
520
521 /* Enforce periodic BC or check that the coordinate is within the domain. */
522 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
523 x1 = -(str->n_x-1)*x1;
524 }
525 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
526 err = 1;
527 }
528 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
529 y1 = -(str->n_y-1)*y1;
530 }
531 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
532 err = 1;
533 }
534 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
535 z1 = -(str->n_z-1)*z1;
536 }
537 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
538 err = 1;
539 }
540
541 if(!err) {
542
543 /* Fetch coefficients explicitly to fetch those that are adjacent
544 subsequently and to store in temporary variables coefficients that
545 will be used multiple times. This is to decrease computational time,
546 by exploiting simultaneous extraction of adjacent memory and by
547 avoiding going through the long str->c array repeatedly. */
548 real c0000 = str->c[n+0];
549 real c0001 = str->c[n+1];
550 real c0002 = str->c[n+2];
551 real c0003 = str->c[n+3];
552 real c0004 = str->c[n+4];
553 real c0005 = str->c[n+5];
554 real c0006 = str->c[n+6];
555 real c0007 = str->c[n+7];
556
557 real c0010 = str->c[n+x1+0];
558 real c0011 = str->c[n+x1+1];
559 real c0012 = str->c[n+x1+2];
560 real c0013 = str->c[n+x1+3];
561 real c0014 = str->c[n+x1+4];
562 real c0015 = str->c[n+x1+5];
563 real c0016 = str->c[n+x1+6];
564 real c0017 = str->c[n+x1+7];
565
566 real c0100 = str->c[n+y1+0];
567 real c0101 = str->c[n+y1+1];
568 real c0102 = str->c[n+y1+2];
569 real c0103 = str->c[n+y1+3];
570 real c0104 = str->c[n+y1+4];
571 real c0105 = str->c[n+y1+5];
572 real c0106 = str->c[n+y1+6];
573 real c0107 = str->c[n+y1+7];
574
575 real c1000 = str->c[n+z1+0];
576 real c1001 = str->c[n+z1+1];
577 real c1002 = str->c[n+z1+2];
578 real c1003 = str->c[n+z1+3];
579 real c1004 = str->c[n+z1+4];
580 real c1005 = str->c[n+z1+5];
581 real c1006 = str->c[n+z1+6];
582 real c1007 = str->c[n+z1+7];
583
584 real c0110 = str->c[n+y1+x1+0];
585 real c0111 = str->c[n+y1+x1+1];
586 real c0112 = str->c[n+y1+x1+2];
587 real c0113 = str->c[n+y1+x1+3];
588 real c0114 = str->c[n+y1+x1+4];
589 real c0115 = str->c[n+y1+x1+5];
590 real c0116 = str->c[n+y1+x1+6];
591 real c0117 = str->c[n+y1+x1+7];
592
593 real c1010 = str->c[n+z1+x1+0];
594 real c1011 = str->c[n+z1+x1+1];
595 real c1012 = str->c[n+z1+x1+2];
596 real c1013 = str->c[n+z1+x1+3];
597 real c1014 = str->c[n+z1+x1+4];
598 real c1015 = str->c[n+z1+x1+5];
599 real c1016 = str->c[n+z1+x1+6];
600 real c1017 = str->c[n+z1+x1+7];
601
602 real c1100 = str->c[n+z1+y1+0];
603 real c1101 = str->c[n+z1+y1+1];
604 real c1102 = str->c[n+z1+y1+2];
605 real c1103 = str->c[n+z1+y1+3];
606 real c1104 = str->c[n+z1+y1+4];
607 real c1105 = str->c[n+z1+y1+5];
608 real c1106 = str->c[n+z1+y1+6];
609 real c1107 = str->c[n+z1+y1+7];
610
611 real c1110 = str->c[n+z1+y1+x1+0];
612 real c1111 = str->c[n+z1+y1+x1+1];
613 real c1112 = str->c[n+z1+y1+x1+2];
614 real c1113 = str->c[n+z1+y1+x1+3];
615 real c1114 = str->c[n+z1+y1+x1+4];
616 real c1115 = str->c[n+z1+y1+x1+5];
617 real c1116 = str->c[n+z1+y1+x1+6];
618 real c1117 = str->c[n+z1+y1+x1+7];
619
620 /* Evaluate spline values */
621
622 /* f */
623 f_df[0] = (
624 dzi*(
625 dxi*(dyi*c0000+dy*c0100)
626 +dx*(dyi*c0010+dy*c0110))
627 +dz*(
628 dxi*(dyi*c1000+dy*c1100)
629 +dx*(dyi*c1010+dy*c1110)))
630 +xg2/6*(
631 dzi*(
632 dxi3*(dyi*c0001+dy*c0101)
633 +dx3*(dyi*c0011+dy*c0111))
634 +dz*(
635 dxi3*(dyi*c1001+dy*c1101)
636 +dx3*(dyi*c1011+dy*c1111)))
637 +yg2/6*(
638 dzi*(
639 dxi*(dyi3*c0002+dy3*c0102)
640 +dx*(dyi3*c0012+dy3*c0112))
641 +dz*(
642 dxi*(dyi3*c1002+dy3*c1102)
643 +dx*(dyi3*c1012+dy3*c1112)))
644 +zg2/6*(
645 dzi3*(
646 dxi*(dyi*c0003+dy*c0103)
647 +dx*(dyi*c0013+dy*c0113))
648 +dz3*(
649 dxi*(dyi*c1003+dy*c1103)
650 +dx*(dyi*c1013+dy*c1113)))
651 +xg2*yg2/36*(
652 dzi*(
653 dxi3*(dyi3*c0004+dy3*c0104)
654 +dx3*(dyi3*c0014+dy3*c0114))
655 +dz*(
656 dxi3*(dyi3*c1004+dy3*c1104)
657 +dx3*(dyi3*c1014+dy3*c1114)))
658 +xg2*zg2/36*(
659 dzi3*(
660 dxi3*(dyi*c0005+dy*c0105)
661 +dx3*(dyi*c0015+dy*c0115))
662 +dz3*(
663 dxi3*(dyi*c1005+dy*c1105)
664 +dx3*(dyi*c1015+dy*c1115)))
665 +yg2*zg2/36*(
666 dzi3*(
667 dxi*(dyi3*c0006+dy3*c0106)
668 +dx*(dyi3*c0016+dy3*c0116))
669 +dz3*(
670 dxi*(dyi3*c1006+dy3*c1106)
671 +dx*(dyi3*c1016+dy3*c1116)))
672 +xg2*yg2*zg2/216*(
673 dzi3*(
674 dxi3*(dyi3*c0007+dy3*c0107)
675 +dx3*(dyi3*c0017+dy3*c0117))
676 +dz3*(
677 dxi3*(dyi3*c1007+dy3*c1107)
678 +dx3*(dyi3*c1017+dy3*c1117)));
679
680 /* df/dx */
681 f_df[1] = xgi*(
682 dzi*(
683 -(dyi*c0000+dy*c0100)
684 +(dyi*c0010+dy*c0110))
685 +dz*(
686 -(dyi*c1000+dy*c1100)
687 +(dyi*c1010+dy*c1110)))
688 +xg/6*(
689 dzi*(
690 dxi3dx*(dyi*c0001+dy*c0101)
691 +dx3dx*(dyi*c0011+dy*c0111))
692 +dz*(
693 dxi3dx*(dyi*c1001 +dy*c1101)
694 +dx3dx*(dyi*c1011+dy*c1111)))
695 +xgi*yg2/6*(
696 dzi*(
697 -(dyi3*c0002+dy3*c0102)
698 +(dyi3*c0012+dy3*c0112))
699 +dz*(
700 -(dyi3*c1002+dy3*c1102)
701 +(dyi3*c1012+dy3*c1112)))
702 +xgi*zg2/6*(
703 dzi3*(
704 -(dyi*c0003+dy*c0103)
705 +(dyi*c0013+dy*c0113))
706 +dz3*(
707 -(dyi*c1003+dy*c1103)
708 +(dyi*c1013+dy*c1113)))
709 +xg*yg2/36*(
710 dzi*(
711 dxi3dx*(dyi3*c0004+dy3*c0104)
712 +dx3dx*(dyi3*c0014+dy3*c0114))
713 +dz*(
714 dxi3dx*(dyi3*c1004+dy3*c1104)
715 +dx3dx*(dyi3*c1014+dy3*c1114)))
716 +xg*zg2/36*(
717 dzi3*(
718 dxi3dx*(dyi*c0005+dy*c0105)
719 +dx3dx*(dyi*c0015+dy*c0115))
720 +dz3*(
721 dxi3dx*(dyi*c1005+dy*c1105)
722 +dx3dx*(dyi*c1015+dy*c1115)))
723 +xgi*yg2*zg2/36*(
724 dzi3*(
725 -(dyi3*c0006+dy3*c0106)
726 +(dyi3*c0016+dy3*c0116))
727 +dz3*(
728 -(dyi3*c1006+dy3*c1106)
729 +(dyi3*c1016+dy3*c1116)))
730 +xg*yg2*zg2/216*(
731 dzi3*(
732 dxi3dx*(dyi3*c0007+dy3*c0107)
733 +dx3dx*(dyi3*c0017+dy3*c0117))
734 +dz3*(
735 dxi3dx*(dyi3*c1007+dy3*c1107)
736 +dx3dx*(dyi3*c1017+dy3*c1117)));
737
738 /* df/dy */
739 f_df[2] = ygi*(
740 dzi*(
741 dxi*(-c0000+c0100)
742 +dx*(-c0010+c0110))
743 +dz*(
744 dxi*(-c1000+c1100)
745 +dx*(-c1010+c1110)))
746 +ygi*xg2/6*(
747 dzi*(
748 dxi3*(-c0001+c0101)
749 +dx3*(-c0011+c0111))
750 +dz*(
751 dxi3*(-c1001+c1101)
752 +dx3*(-c1011+c1111)))
753 +yg/6*(
754 dzi*(
755 dxi*(dyi3dy*c0002+dy3dy*c0102)
756 +dx*(dyi3dy*c0012+dy3dy*c0112))
757 +dz*(
758 dxi*(dyi3dy*c1002+dy3dy*c1102)
759 +dx*(dyi3dy*c1012+dy3dy*c1112)))
760 +ygi*zg2/6*(
761 dzi3*(
762 dxi*(-c0003+c0103)
763 +dx*(-c0013+c0113))
764 +dz3*(
765 dxi*(-c1003+c1103)
766 +dx*(-c1013+c1113)))
767 +xg2*yg/36*(
768 dzi*(
769 dxi3*(dyi3dy*c0004+dy3dy*c0104)
770 +dx3*(dyi3dy*c0014+dy3dy*c0114))
771 +dz*(
772 dxi3*(dyi3dy*c1004+dy3dy*c1104)
773 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
774 +ygi*xg2*zg2/36*(
775 dzi3*(
776 dxi3*(-c0005+c0105)
777 +dx3*(-c0015+c0115))
778 +dz3*(
779 dxi3*(-c1005+c1105)
780 +dx3*(-c1015+c1115)))
781 +yg*zg2/36*(
782 dzi3*(
783 dxi*(dyi3dy*c0006+dy3dy*c0106)
784 +dx*(dyi3dy*c0016+dy3dy*c0116))
785 +dz3*(
786 dxi*(dyi3dy*c1006+dy3dy*c1106)
787 +dx*(dyi3dy*c1016+dy3dy*c1116)))
788 +xg2*yg*zg2/216*(
789 dzi3*(
790 dxi3*(dyi3dy*c0007+dy3dy*c0107)
791 +dx3*(dyi3dy*c0017+dy3dy*c0117))
792 +dz3*(
793 dxi3*(dyi3dy*c1007+dy3dy*c1107)
794 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
795
796 /* df/dz */
797 f_df[3] = zgi*(
798 -(
799 dxi*(dyi*c0000+dy*c0100)
800 +dx*(dyi*c0010+dy*c0110))
801 +(
802 dxi*(dyi*c1000+dy*c1100)
803 +dx*(dyi*c1010+dy*c1110)))
804 +xg2*zgi/6*(
805 -(
806 dxi3*(dyi*c0001+dy*c0101)
807 +dx3*(dyi*c0011+dy*c0111))
808 +(
809 dxi3*(dyi*c1001+dy*c1101)
810 +dx3*(dyi*c1011+dy*c1111)))
811 +yg2*zgi/6*(
812 -(
813 dxi*(dyi3*c0002+dy3*c0102)
814 +dx*(dyi3*c0012+dy3*c0112))
815 +(
816 dxi*(dyi3*c1002+dy3*c1102)
817 +dx*(dyi3*c1012+dy3*c1112)))
818 +zg/6*(
819 dzi3dz*(
820 dxi*(dyi*c0003+dy*c0103)
821 +dx*(dyi*c0013+dy*c0113))
822 +dz3dz*(
823 dxi*(dyi*c1003+dy*c1103)
824 +dx*(dyi*c1013+dy*c1113)))
825 +xg2*yg2*zgi/36*(
826 -(
827 dxi3*(dyi3*c0004+dy3*c0104)
828 +dx3*(dyi3*c0014+dy3*c0114))
829 +(
830 dxi3*(dyi3*c1004+dy3*c1104)
831 +dx3*(dyi3*c1014+dy3*c1114)))
832 +xg2*zg/36*(
833 dzi3dz*(
834 dxi3*(dyi*c0005+dy*c0105)
835 +dx3*(dyi*c0015+dy*c0115))
836 +dz3dz*(
837 dxi3*(dyi*c1005+dy*c1105)
838 +dx3*(dyi*c1015+dy*c1115)))
839 +yg2*zg/36*(
840 dzi3dz*(
841 dxi*(dyi3*c0006+dy3*c0106)
842 +dx*(dyi3*c0016+dy3*c0116))
843 +dz3dz*(
844 dxi*(dyi3*c1006+dy3*c1106)
845 +dx*(dyi3*c1016+dy3*c1116)))
846 +xg2*yg2*zg/216*(
847 dzi3dz*(
848 dxi3*(dyi3*c0007+dy3*c0107)
849 +dx3*(dyi3*c0017+dy3*c0117))
850 +dz3dz*(
851 dxi3*(dyi3*c1007+dy3*c1107)
852 +dx3*(dyi3*c1017+dy3*c1117)));
853
854 /* d2f/dx2 */
855 f_df[4] = (
856 dzi*(
857 dxi*(dyi*c0001+dy*c0101)
858 +dx*(dyi*c0011+dy*c0111))
859 +dz*(
860 dxi*(dyi*c1001+dy*c1101)
861 +dx*(dyi*c1011+dy*c1111)))
862 +yg2/6*(
863 dzi*(
864 dxi*(dyi3*c0004+dy3*c0104)
865 +dx*(dyi3*c0014+dy3*c0114))
866 +dz*(
867 dxi*(dyi3*c1004+dy3*c1104)
868 +dx*(dyi3*c1014+dy3*c1114)))
869 +zg2/6*(
870 dzi3*(
871 dxi*(dyi*c0005+dy*c0105)
872 +dx*(dyi*c0015+dy*c0115))
873 +dz3*(
874 dxi*(dyi*c1005+dy*c1105)
875 +dx*(dyi*c1015+dy*c1115)))
876 +yg2*zg2/36*(
877 dzi3*(
878 dxi*(dyi3*c0007+dy3*c0107)
879 +dx*(dyi3*c0017+dy3*c0117))
880 +dz3*(
881 dxi*(dyi3*c1007+dy3*c1107)
882 +dx*(dyi3*c1017+dy3*c1117)));
883
884 /* d2f/dy2 */
885 f_df[5] = (
886 dzi*(
887 dxi*(dyi*c0002+dy*c0102)
888 +dx*(dyi*c0012+dy*c0112))
889 +dz*(
890 dxi*(dyi*c1002+dy*c1102)
891 +dx*(dyi*c1012+dy*c1112)))
892 +xg2/6*(
893 dzi*(
894 dxi3*(dyi*c0004+dy*c0104)
895 +dx3*(dyi*c0014+dy*c0114))
896 +dz*(
897 dxi3*(dyi*c1004+dy*c1104)
898 +dx3*(dyi*c1014+dy*c1114)))
899 +zg2/6*(
900 dzi3*(
901 dxi*(dyi*c0006+dy*c0106)
902 +dx*(dyi*c0016+dy*c0116))
903 +dz3*(
904 dxi*(dyi*c1006+dy*c1106)
905 +dx*(dyi*c1016+dy*c1116)))
906 +xg2*zg2/36*(
907 dzi3*(
908 dxi3*(dyi*c0007+dy*c0107)
909 +dx3*(dyi*c0017+dy*c0117))
910 +dz3*(
911 dxi3*(dyi*c1007+dy*c1107)
912 +dx3*(dyi*c1017+dy*c1117)));
913
914 /* d2f/dz2 */
915 f_df[6] = (
916 dzi*(
917 dxi*(dyi*c0003+dy*c0103)
918 +dx*(dyi*c0013+dy*c0113))
919 +dz*(
920 dxi*(dyi*c1003+dy*c1103)
921 +dx*(dyi*c1013+dy*c1113)))
922 +xg2/6*(
923 dzi*(
924 dxi3*(dyi*c0005+dy*c0105)
925 +dx3*(dyi*c0015+dy*c0115))
926 +dz*(
927 dxi3*(dyi*c1005+dy*c1105)
928 +dx3*(dyi*c1015+dy*c1115)))
929 +yg2/6*(
930 dzi*(
931 dxi*(dyi3*c0006+dy3*c0106)
932 +dx*(dyi3*c0016+dy3*c0116))
933 +dz*(
934 dxi*(dyi3*c1006+dy3*c1106)
935 +dx*(dyi3*c1016+dy3*c1116)))
936 +xg2*yg2/36*(
937 dzi*(
938 dxi3*(dyi3*c0007+dy3*c0107)
939 +dx3*(dyi3*c0017+dy3*c0117))
940 +dz*(
941 dxi3*(dyi3*c1007+dy3*c1107)
942 +dx3*(dyi3*c1017+dy3*c1117)));
943
944 /* d2f/dxdy */
945 f_df[7] = xgi*ygi*(
946 dzi*(
947 (c0000 -c0100)
948 -(c0010-c0110))
949 +dz*(
950 (c1000 -c1100)
951 -(c1010-c1110)))
952 +ygi*xg/6*(
953 dzi*(
954 dxi3dx*(-c0001+c0101)
955 +dx3dx*(-c0011+c0111))
956 +dz*(
957 dxi3dx*(-c1001+c1101)
958 +dx3dx*(-c1011+c1111)))
959 +xgi*yg/6*(
960 dzi*(
961 -(dyi3dy*c0002+dy3dy*c0102)
962 +(dyi3dy*c0012+dy3dy*c0112))
963 +dz*(
964 -(dyi3dy*c1002+dy3dy*c1102)
965 +(dyi3dy*c1012+dy3dy*c1112)))
966 +xgi*ygi*zg2/6*(
967 dzi3*(
968 (c0003 -c0103)
969 -(c0013-c0113))
970 +dz3*(
971 (c1003 -c1103)
972 -(c1013-c1113)))
973 +xg*yg/36*(
974 dzi*(
975 dxi3dx*(dyi3dy*c0004+dy3dy*c0104)
976 +dx3dx*(dyi3dy*c0014+dy3dy*c0114))
977 +dz*(
978 dxi3dx*(dyi3dy*c1004+dy3dy*c1104)
979 +dx3dx*(dyi3dy*c1014+dy3dy*c1114)))
980 +ygi*xg*zg2/36*(
981 dzi3*(
982 dxi3dx*(-c0005+c0105)
983 +dx3dx*(-c0015+c0115))
984 +dz3*(
985 dxi3dx*(-c1005+c1105)
986 +dx3dx*(-c1015+c1115)))
987 +xgi*yg*zg2/36*(
988 dzi3*(
989 -(dyi3dy*c0006+dy3dy*c0106)
990 +(dyi3dy*c0016+dy3dy*c0116))
991 +dz3*(
992 -(dyi3dy*c1006+dy3dy*c1106)
993 +(dyi3dy*c1016+dy3dy*c1116)))
994 +xg*yg*zg2/216*(
995 dzi3*(
996 dxi3dx*(dyi3dy*c0007+dy3dy*c0107)
997 +dx3dx*(dyi3dy*c0017+dy3dy*c0117))
998 +dz3*(
999 dxi3dx*(dyi3dy*c1007+dy3dy*c1107)
1000 +dx3dx*(dyi3dy*c1017+dy3dy*c1117)));
1001
1002 /* d2f/dxdz */
1003 f_df[8] = xgi*zgi*(
1004 (
1005 (dyi*c0000+dy*c0100)
1006 -(dyi*c0010+dy*c0110))
1007 -(
1008 (dyi*c1000+dy*c1100)
1009 -(dyi*c1010+dy*c1110)))
1010 +xg*zgi/6*(
1011 -(
1012 dxi3dx*(dyi*c0001+dy*c0101)
1013 +dx3dx*(dyi*c0011+dy*c0111))
1014 +(
1015 dxi3dx*(dyi*c1001+dy*c1101)
1016 +dx3dx*(dyi*c1011+dy*c1111)))
1017 +xgi*yg2*zgi/6*(
1018 (
1019 (dyi3*c0002+dy3*c0102)
1020 -(dyi3*c0012+dy3*c0112))
1021 -(
1022 (dyi3*c1002+dy3*c1102)
1023 -(dyi3*c1012+dy3*c1112)))
1024 +xgi*zg/6*(
1025 dzi3dz*(
1026 -(dyi*c0003+dy*c0103)
1027 +(dyi*c0013+dy*c0113))
1028 +dz3dz*(
1029 -(dyi*c1003+dy*c1103)
1030 +(dyi*c1013+dy*c1113)))
1031 +xg*yg2*zgi/36*(
1032 -(
1033 dxi3dx*(dyi3*c0004+dy3*c0104)
1034 +dx3dx*(dyi3*c0014+dy3*c0114))
1035 +(
1036 dxi3dx*(dyi3*c1004+dy3*c1104)
1037 +dx3dx*(dyi3*c1014+dy3*c1114)))
1038 +xg*zg/36*(
1039 dzi3dz*(
1040 dxi3dx*(dyi*c0005+dy*c0105)
1041 +dx3dx*(dyi*c0015+dy*c0115))
1042 +dz3dz*(
1043 dxi3dx*(dyi*c1005+dy*c1105)
1044 +dx3dx*(dyi*c1015+dy*c1115)))
1045 +xgi*yg2*zg/36*(
1046 dzi3dz*(
1047 -(dyi3*c0006+dy3*c0106)
1048 +(dyi3*c0016+dy3*c0116))
1049 +dz3dz*(
1050 -(dyi3*c1006+dy3*c1106)
1051 +(dyi3*c1016+dy3*c1116)))
1052 +xg*yg2*zg/216*(
1053 dzi3dz*(
1054 dxi3dx*(dyi3*c0007+dy3*c0107)
1055 +dx3dx*(dyi3*c0017+dy3*c0117))
1056 +dz3dz*(
1057 dxi3dx*(dyi3*c1007+dy3*c1107)
1058 +dx3dx*(dyi3*c1017+dy3*c1117)));
1059
1060 /* d2f/dydz */
1061 f_df[9] = ygi*zgi*(
1062 (
1063 dxi*(c0000 -c0100)
1064 +dx*(c0010-c0110))
1065 -(
1066 dxi*(c1000 -c1100)
1067 +dx*(c1010-c1110)))
1068 +ygi*xg2*zgi/6*(
1069 (
1070 dxi3*(c0001 -c0101)
1071 +dx3*(c0011-c0111))
1072 -(
1073 dxi3*(c1001 -c1101)
1074 +dx3*(c1011-c1111)))
1075 +yg*zgi/6*(
1076 -(
1077 dxi*(dyi3dy*c0002+dy3dy*c0102)
1078 +dx*(dyi3dy*c0012+dy3dy*c0112))
1079 +(
1080 dxi*(dyi3dy*c1002+dy3dy*c1102)
1081 +dx*(dyi3dy*c1012+dy3dy*c1112)))
1082 +ygi*zg/6*(
1083 dzi3dz*(
1084 dxi*(-c0003+c0103)
1085 +dx*(-c0013+c0113))
1086 +dz3dz*(
1087 dxi*(-c1003+c1103)
1088 +dx*(-c1013+c1113)))
1089 +xg2*yg*zgi/36*(
1090 -(
1091 dxi3*(dyi3dy*c0004+dy3dy*c0104)
1092 +dx3*(dyi3dy*c0014+dy3dy*c0114))
1093 +(
1094 dxi3*(dyi3dy*c1004+dy3dy*c1104)
1095 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
1096 +ygi*xg2*zg/36*(
1097 dzi3dz*(
1098 dxi3*(-c0005+c0105)
1099 +dx3*(-c0015+c0115))
1100 +dz3dz*(
1101 dxi3*(-c1005+c1105)
1102 +dx3*(-c1015+c1115)))
1103 +yg*zg/36*(
1104 dzi3dz*(
1105 dxi*(dyi3dy*c0006+dy3dy*c0106)
1106 +dx*(dyi3dy*c0016+dy3dy*c0116))
1107 +dz3dz*(
1108 dxi*(dyi3dy*c1006+dy3dy*c1106)
1109 +dx*(dyi3dy*c1016+dy3dy*c1116)))
1110 +xg2*yg*zg/216*(
1111 dzi3dz*(
1112 dxi3*(dyi3dy*c0007+dy3dy*c0107)
1113 +dx3*(dyi3dy*c0017+dy3dy*c0117))
1114 +dz3dz*(
1115 dxi3*(dyi3dy*c1007+dy3dy*c1107)
1116 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
1117 }
1118
1119 return err;
1120}
Main header file for ASCOT5.
double real
Definition ascot5.h:85
Header file containing physical and mathematical constants.
unsigned long int a5err
Simulation error flag.
Definition error.h:17
int interp3Dcomp_setup(interp3D_data *str, 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)
Set up splines to interpolate 3D scalar data.
int interp3Dcomp_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.
a5err interp3Dcomp_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.
void interp3Dcomp_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.
a5err interp3Dcomp_eval_f(real *f, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D scalar field.
Spline interpolation library.
@ NATURALBC
Definition interp.h:37
@ PERIODICBC
Definition interp.h:38
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
Header file for math.c.
Header file for splineexpl.c and splinecomp.c.
void splinecomp(real *f, int n, int bc, real *c)
Calculate compact cubic spline interpolation coefficients in 1D.
Definition splinecomp.c:22
Tricubic interpolation struct.
Definition interp.h:85
real z_max
Definition interp.h:99
real z_min
Definition interp.h:98
real z_grid
Definition interp.h:100
real x_min
Definition interp.h:92
real y_max
Definition interp.h:96
real y_grid
Definition interp.h:97
real * c
Definition interp.h:101
real x_grid
Definition interp.h:94
real x_max
Definition interp.h:93
real y_min
Definition interp.h:95