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
37int interp3Dcomp_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 /* 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
203void interp3Dcomp_init_spline(interp3D_data* str, real* c,
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
259int interp3Dcomp_setup(interp3D_data* str, real* f,
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
288a5err interp3Dcomp_eval_f(real* f, interp3D_data* str, real x, real y, real z) {
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
448a5err interp3Dcomp_eval_df(real* f_df, interp3D_data* str,
449 real x, real y, real z) {
450
451 /* Make sure periodic coordinates are within [min, max] region. */
452 if(str->bc_x == PERIODICBC) {
453 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
454 x = x + (x < str->x_min) * (str->x_max - str->x_min);
455 }
456 if(str->bc_y == PERIODICBC) {
457 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
458 y = y + (y < str->y_min) * (str->y_max - str->y_min);
459 }
460 if(str->bc_z == PERIODICBC) {
461 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
462 z = z + (z < str->z_min) * (str->z_max - str->z_min);
463 }
464
465 /* Index for x variable. The -1 needed at exactly grid end. */
466 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
467 /* Normalized x coordinate in current cell */
468 real dx = ( x - (str->x_min + i_x*str->x_grid) ) / str->x_grid;
469 /* Helper variables */
470 real dx3 = dx*dx*dx - dx;
471 real dx3dx = 3*dx*dx - 1.0;
472 real dxi = 1.0 - dx;
473 real dxi3 = dxi*dxi*dxi - dxi;
474 real dxi3dx = -3*dxi*dxi + 1.0;
475 real xg = str->x_grid;
476 real xg2 = xg*xg;
477 real xgi = 1.0 / xg;
478
479 /* Index for y variable. The -1 needed at exactly grid end. */
480 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
481 /* Normalized y coordinate in current cell */
482 real dy = ( y - (str->y_min + i_y*str->y_grid) ) / str->y_grid;
483 /* Helper variables */
484 real dy3 = dy*dy*dy-dy;
485 real dy3dy = 3*dy*dy - 1.0;
486 real dyi = 1.0 - dy;
487 real dyi3 = dyi*dyi*dyi - dyi;
488 real dyi3dy = -3*dyi*dyi + 1.0;
489 real yg = str->y_grid;
490 real yg2 = yg*yg;
491 real ygi = 1.0 / yg;
492
493 /* Index for z variable. The -1 needed at exactly grid end. */
494 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
495 /* Normalized z coordinate in current cell */
496 real dz = ( z - (str->z_min + i_z*str->z_grid) ) / str->z_grid;
497 /* Helper variables */
498 real dz3 = dz*dz*dz - dz;
499 real dz3dz = 3*dz*dz - 1.0;
500 real dzi = 1.0 - dz;
501 real dzi3 = dzi*dzi*dzi - dzi;
502 real dzi3dz = -3*dzi*dzi + 1.0;
503 real zg = str->z_grid;
504 real zg2 = zg*zg;
505 real zgi = 1.0 / zg;
506
507 /* Index jump to cell */
508 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
509 int x1 = 8; /* Index jump one x forward */
510 int y1 = str->n_x*8; /* Index jump one y forward */
511 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
512
513 int err = 0;
514
515 /* Enforce periodic BC or check that the coordinate is within the domain. */
516 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
517 x1 = -(str->n_x-1)*x1;
518 }
519 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
520 err = 1;
521 }
522 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
523 y1 = -(str->n_y-1)*y1;
524 }
525 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
526 err = 1;
527 }
528 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
529 z1 = -(str->n_z-1)*z1;
530 }
531 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
532 err = 1;
533 }
534
535 if(!err) {
536
537 /* Fetch coefficients explicitly to fetch those that are adjacent
538 subsequently and to store in temporary variables coefficients that
539 will be used multiple times. This is to decrease computational time,
540 by exploiting simultaneous extraction of adjacent memory and by
541 avoiding going through the long str->c array repeatedly. */
542 real c0000 = str->c[n+0];
543 real c0001 = str->c[n+1];
544 real c0002 = str->c[n+2];
545 real c0003 = str->c[n+3];
546 real c0004 = str->c[n+4];
547 real c0005 = str->c[n+5];
548 real c0006 = str->c[n+6];
549 real c0007 = str->c[n+7];
550
551 real c0010 = str->c[n+x1+0];
552 real c0011 = str->c[n+x1+1];
553 real c0012 = str->c[n+x1+2];
554 real c0013 = str->c[n+x1+3];
555 real c0014 = str->c[n+x1+4];
556 real c0015 = str->c[n+x1+5];
557 real c0016 = str->c[n+x1+6];
558 real c0017 = str->c[n+x1+7];
559
560 real c0100 = str->c[n+y1+0];
561 real c0101 = str->c[n+y1+1];
562 real c0102 = str->c[n+y1+2];
563 real c0103 = str->c[n+y1+3];
564 real c0104 = str->c[n+y1+4];
565 real c0105 = str->c[n+y1+5];
566 real c0106 = str->c[n+y1+6];
567 real c0107 = str->c[n+y1+7];
568
569 real c1000 = str->c[n+z1+0];
570 real c1001 = str->c[n+z1+1];
571 real c1002 = str->c[n+z1+2];
572 real c1003 = str->c[n+z1+3];
573 real c1004 = str->c[n+z1+4];
574 real c1005 = str->c[n+z1+5];
575 real c1006 = str->c[n+z1+6];
576 real c1007 = str->c[n+z1+7];
577
578 real c0110 = str->c[n+y1+x1+0];
579 real c0111 = str->c[n+y1+x1+1];
580 real c0112 = str->c[n+y1+x1+2];
581 real c0113 = str->c[n+y1+x1+3];
582 real c0114 = str->c[n+y1+x1+4];
583 real c0115 = str->c[n+y1+x1+5];
584 real c0116 = str->c[n+y1+x1+6];
585 real c0117 = str->c[n+y1+x1+7];
586
587 real c1010 = str->c[n+z1+x1+0];
588 real c1011 = str->c[n+z1+x1+1];
589 real c1012 = str->c[n+z1+x1+2];
590 real c1013 = str->c[n+z1+x1+3];
591 real c1014 = str->c[n+z1+x1+4];
592 real c1015 = str->c[n+z1+x1+5];
593 real c1016 = str->c[n+z1+x1+6];
594 real c1017 = str->c[n+z1+x1+7];
595
596 real c1100 = str->c[n+z1+y1+0];
597 real c1101 = str->c[n+z1+y1+1];
598 real c1102 = str->c[n+z1+y1+2];
599 real c1103 = str->c[n+z1+y1+3];
600 real c1104 = str->c[n+z1+y1+4];
601 real c1105 = str->c[n+z1+y1+5];
602 real c1106 = str->c[n+z1+y1+6];
603 real c1107 = str->c[n+z1+y1+7];
604
605 real c1110 = str->c[n+z1+y1+x1+0];
606 real c1111 = str->c[n+z1+y1+x1+1];
607 real c1112 = str->c[n+z1+y1+x1+2];
608 real c1113 = str->c[n+z1+y1+x1+3];
609 real c1114 = str->c[n+z1+y1+x1+4];
610 real c1115 = str->c[n+z1+y1+x1+5];
611 real c1116 = str->c[n+z1+y1+x1+6];
612 real c1117 = str->c[n+z1+y1+x1+7];
613
614 /* Evaluate spline values */
615
616 /* f */
617 f_df[0] = (
618 dzi*(
619 dxi*(dyi*c0000+dy*c0100)
620 +dx*(dyi*c0010+dy*c0110))
621 +dz*(
622 dxi*(dyi*c1000+dy*c1100)
623 +dx*(dyi*c1010+dy*c1110)))
624 +xg2/6*(
625 dzi*(
626 dxi3*(dyi*c0001+dy*c0101)
627 +dx3*(dyi*c0011+dy*c0111))
628 +dz*(
629 dxi3*(dyi*c1001+dy*c1101)
630 +dx3*(dyi*c1011+dy*c1111)))
631 +yg2/6*(
632 dzi*(
633 dxi*(dyi3*c0002+dy3*c0102)
634 +dx*(dyi3*c0012+dy3*c0112))
635 +dz*(
636 dxi*(dyi3*c1002+dy3*c1102)
637 +dx*(dyi3*c1012+dy3*c1112)))
638 +zg2/6*(
639 dzi3*(
640 dxi*(dyi*c0003+dy*c0103)
641 +dx*(dyi*c0013+dy*c0113))
642 +dz3*(
643 dxi*(dyi*c1003+dy*c1103)
644 +dx*(dyi*c1013+dy*c1113)))
645 +xg2*yg2/36*(
646 dzi*(
647 dxi3*(dyi3*c0004+dy3*c0104)
648 +dx3*(dyi3*c0014+dy3*c0114))
649 +dz*(
650 dxi3*(dyi3*c1004+dy3*c1104)
651 +dx3*(dyi3*c1014+dy3*c1114)))
652 +xg2*zg2/36*(
653 dzi3*(
654 dxi3*(dyi*c0005+dy*c0105)
655 +dx3*(dyi*c0015+dy*c0115))
656 +dz3*(
657 dxi3*(dyi*c1005+dy*c1105)
658 +dx3*(dyi*c1015+dy*c1115)))
659 +yg2*zg2/36*(
660 dzi3*(
661 dxi*(dyi3*c0006+dy3*c0106)
662 +dx*(dyi3*c0016+dy3*c0116))
663 +dz3*(
664 dxi*(dyi3*c1006+dy3*c1106)
665 +dx*(dyi3*c1016+dy3*c1116)))
666 +xg2*yg2*zg2/216*(
667 dzi3*(
668 dxi3*(dyi3*c0007+dy3*c0107)
669 +dx3*(dyi3*c0017+dy3*c0117))
670 +dz3*(
671 dxi3*(dyi3*c1007+dy3*c1107)
672 +dx3*(dyi3*c1017+dy3*c1117)));
673
674 /* df/dx */
675 f_df[1] = xgi*(
676 dzi*(
677 -(dyi*c0000+dy*c0100)
678 +(dyi*c0010+dy*c0110))
679 +dz*(
680 -(dyi*c1000+dy*c1100)
681 +(dyi*c1010+dy*c1110)))
682 +xg/6*(
683 dzi*(
684 dxi3dx*(dyi*c0001+dy*c0101)
685 +dx3dx*(dyi*c0011+dy*c0111))
686 +dz*(
687 dxi3dx*(dyi*c1001 +dy*c1101)
688 +dx3dx*(dyi*c1011+dy*c1111)))
689 +xgi*yg2/6*(
690 dzi*(
691 -(dyi3*c0002+dy3*c0102)
692 +(dyi3*c0012+dy3*c0112))
693 +dz*(
694 -(dyi3*c1002+dy3*c1102)
695 +(dyi3*c1012+dy3*c1112)))
696 +xgi*zg2/6*(
697 dzi3*(
698 -(dyi*c0003+dy*c0103)
699 +(dyi*c0013+dy*c0113))
700 +dz3*(
701 -(dyi*c1003+dy*c1103)
702 +(dyi*c1013+dy*c1113)))
703 +xg*yg2/36*(
704 dzi*(
705 dxi3dx*(dyi3*c0004+dy3*c0104)
706 +dx3dx*(dyi3*c0014+dy3*c0114))
707 +dz*(
708 dxi3dx*(dyi3*c1004+dy3*c1104)
709 +dx3dx*(dyi3*c1014+dy3*c1114)))
710 +xg*zg2/36*(
711 dzi3*(
712 dxi3dx*(dyi*c0005+dy*c0105)
713 +dx3dx*(dyi*c0015+dy*c0115))
714 +dz3*(
715 dxi3dx*(dyi*c1005+dy*c1105)
716 +dx3dx*(dyi*c1015+dy*c1115)))
717 +xgi*yg2*zg2/36*(
718 dzi3*(
719 -(dyi3*c0006+dy3*c0106)
720 +(dyi3*c0016+dy3*c0116))
721 +dz3*(
722 -(dyi3*c1006+dy3*c1106)
723 +(dyi3*c1016+dy3*c1116)))
724 +xg*yg2*zg2/216*(
725 dzi3*(
726 dxi3dx*(dyi3*c0007+dy3*c0107)
727 +dx3dx*(dyi3*c0017+dy3*c0117))
728 +dz3*(
729 dxi3dx*(dyi3*c1007+dy3*c1107)
730 +dx3dx*(dyi3*c1017+dy3*c1117)));
731
732 /* df/dy */
733 f_df[2] = ygi*(
734 dzi*(
735 dxi*(-c0000+c0100)
736 +dx*(-c0010+c0110))
737 +dz*(
738 dxi*(-c1000+c1100)
739 +dx*(-c1010+c1110)))
740 +ygi*xg2/6*(
741 dzi*(
742 dxi3*(-c0001+c0101)
743 +dx3*(-c0011+c0111))
744 +dz*(
745 dxi3*(-c1001+c1101)
746 +dx3*(-c1011+c1111)))
747 +yg/6*(
748 dzi*(
749 dxi*(dyi3dy*c0002+dy3dy*c0102)
750 +dx*(dyi3dy*c0012+dy3dy*c0112))
751 +dz*(
752 dxi*(dyi3dy*c1002+dy3dy*c1102)
753 +dx*(dyi3dy*c1012+dy3dy*c1112)))
754 +ygi*zg2/6*(
755 dzi3*(
756 dxi*(-c0003+c0103)
757 +dx*(-c0013+c0113))
758 +dz3*(
759 dxi*(-c1003+c1103)
760 +dx*(-c1013+c1113)))
761 +xg2*yg/36*(
762 dzi*(
763 dxi3*(dyi3dy*c0004+dy3dy*c0104)
764 +dx3*(dyi3dy*c0014+dy3dy*c0114))
765 +dz*(
766 dxi3*(dyi3dy*c1004+dy3dy*c1104)
767 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
768 +ygi*xg2*zg2/36*(
769 dzi3*(
770 dxi3*(-c0005+c0105)
771 +dx3*(-c0015+c0115))
772 +dz3*(
773 dxi3*(-c1005+c1105)
774 +dx3*(-c1015+c1115)))
775 +yg*zg2/36*(
776 dzi3*(
777 dxi*(dyi3dy*c0006+dy3dy*c0106)
778 +dx*(dyi3dy*c0016+dy3dy*c0116))
779 +dz3*(
780 dxi*(dyi3dy*c1006+dy3dy*c1106)
781 +dx*(dyi3dy*c1016+dy3dy*c1116)))
782 +xg2*yg*zg2/216*(
783 dzi3*(
784 dxi3*(dyi3dy*c0007+dy3dy*c0107)
785 +dx3*(dyi3dy*c0017+dy3dy*c0117))
786 +dz3*(
787 dxi3*(dyi3dy*c1007+dy3dy*c1107)
788 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
789
790 /* df/dz */
791 f_df[3] = zgi*(
792 -(
793 dxi*(dyi*c0000+dy*c0100)
794 +dx*(dyi*c0010+dy*c0110))
795 +(
796 dxi*(dyi*c1000+dy*c1100)
797 +dx*(dyi*c1010+dy*c1110)))
798 +xg2*zgi/6*(
799 -(
800 dxi3*(dyi*c0001+dy*c0101)
801 +dx3*(dyi*c0011+dy*c0111))
802 +(
803 dxi3*(dyi*c1001+dy*c1101)
804 +dx3*(dyi*c1011+dy*c1111)))
805 +yg2*zgi/6*(
806 -(
807 dxi*(dyi3*c0002+dy3*c0102)
808 +dx*(dyi3*c0012+dy3*c0112))
809 +(
810 dxi*(dyi3*c1002+dy3*c1102)
811 +dx*(dyi3*c1012+dy3*c1112)))
812 +zg/6*(
813 dzi3dz*(
814 dxi*(dyi*c0003+dy*c0103)
815 +dx*(dyi*c0013+dy*c0113))
816 +dz3dz*(
817 dxi*(dyi*c1003+dy*c1103)
818 +dx*(dyi*c1013+dy*c1113)))
819 +xg2*yg2*zgi/36*(
820 -(
821 dxi3*(dyi3*c0004+dy3*c0104)
822 +dx3*(dyi3*c0014+dy3*c0114))
823 +(
824 dxi3*(dyi3*c1004+dy3*c1104)
825 +dx3*(dyi3*c1014+dy3*c1114)))
826 +xg2*zg/36*(
827 dzi3dz*(
828 dxi3*(dyi*c0005+dy*c0105)
829 +dx3*(dyi*c0015+dy*c0115))
830 +dz3dz*(
831 dxi3*(dyi*c1005+dy*c1105)
832 +dx3*(dyi*c1015+dy*c1115)))
833 +yg2*zg/36*(
834 dzi3dz*(
835 dxi*(dyi3*c0006+dy3*c0106)
836 +dx*(dyi3*c0016+dy3*c0116))
837 +dz3dz*(
838 dxi*(dyi3*c1006+dy3*c1106)
839 +dx*(dyi3*c1016+dy3*c1116)))
840 +xg2*yg2*zg/216*(
841 dzi3dz*(
842 dxi3*(dyi3*c0007+dy3*c0107)
843 +dx3*(dyi3*c0017+dy3*c0117))
844 +dz3dz*(
845 dxi3*(dyi3*c1007+dy3*c1107)
846 +dx3*(dyi3*c1017+dy3*c1117)));
847 }
848
849 return err;
850}
851
879a5err interp3Dcomp_eval_ddf(real* f_df, interp3D_data* str,
880 real x, real y, real z) {
881
882 /* Make sure periodic coordinates are within [min, max] region. */
883 if(str->bc_x == PERIODICBC) {
884 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
885 x = x + (x < str->x_min) * (str->x_max - str->x_min);
886 }
887 if(str->bc_y == PERIODICBC) {
888 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
889 y = y + (y < str->y_min) * (str->y_max - str->y_min);
890 }
891 if(str->bc_z == PERIODICBC) {
892 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
893 z = z + (z < str->z_min) * (str->z_max - str->z_min);
894 }
895
896 /* Index for x variable. The -1 needed at exactly grid end. */
897 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
898 /* Normalized x coordinate in current cell */
899 real dx = ( x - (str->x_min + i_x*str->x_grid) ) / str->x_grid;
900 /* Helper variables */
901 real dx3 = dx*dx*dx - dx;
902 real dx3dx = 3*dx*dx - 1.0;
903 real dxi = 1.0 - dx;
904 real dxi3 = dxi*dxi*dxi - dxi;
905 real dxi3dx = -3*dxi*dxi + 1.0;
906 real xg = str->x_grid;
907 real xg2 = xg*xg;
908 real xgi = 1.0 / xg;
909
910 /* Index for y variable. The -1 needed at exactly grid end. */
911 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
912 /* Normalized y coordinate in current cell */
913 real dy = ( y - (str->y_min + i_y*str->y_grid) ) / str->y_grid;
914 /* Helper variables */
915 real dy3 = dy*dy*dy-dy;
916 real dy3dy = 3*dy*dy - 1.0;
917 real dyi = 1.0 - dy;
918 real dyi3 = dyi*dyi*dyi - dyi;
919 real dyi3dy = -3*dyi*dyi + 1.0;
920 real yg = str->y_grid;
921 real yg2 = yg*yg;
922 real ygi = 1.0 / yg;
923
924 /* Index for z variable. The -1 needed at exactly grid end. */
925 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
926 /* Normalized z coordinate in current cell */
927 real dz = ( z - (str->z_min + i_z*str->z_grid) ) / str->z_grid;
928 /* Helper variables */
929 real dz3 = dz*dz*dz - dz;
930 real dz3dz = 3*dz*dz - 1.0;
931 real dzi = 1.0 - dz;
932 real dzi3 = dzi*dzi*dzi - dzi;
933 real dzi3dz = -3*dzi*dzi + 1.0;
934 real zg = str->z_grid;
935 real zg2 = zg*zg;
936 real zgi = 1.0 / zg;
937
938 /* Index jump to cell */
939 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
940 int x1 = 8; /* Index jump one x forward */
941 int y1 = str->n_x*8; /* Index jump one y forward */
942 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
943
944 int err = 0;
945
946 /* Enforce periodic BC or check that the coordinate is within the domain. */
947 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
948 x1 = -(str->n_x-1)*x1;
949 }
950 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
951 err = 1;
952 }
953 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
954 y1 = -(str->n_y-1)*y1;
955 }
956 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
957 err = 1;
958 }
959 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
960 z1 = -(str->n_z-1)*z1;
961 }
962 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
963 err = 1;
964 }
965
966 if(!err) {
967
968 /* Fetch coefficients explicitly to fetch those that are adjacent
969 subsequently and to store in temporary variables coefficients that
970 will be used multiple times. This is to decrease computational time,
971 by exploiting simultaneous extraction of adjacent memory and by
972 avoiding going through the long str->c array repeatedly. */
973 real c0000 = str->c[n+0];
974 real c0001 = str->c[n+1];
975 real c0002 = str->c[n+2];
976 real c0003 = str->c[n+3];
977 real c0004 = str->c[n+4];
978 real c0005 = str->c[n+5];
979 real c0006 = str->c[n+6];
980 real c0007 = str->c[n+7];
981
982 real c0010 = str->c[n+x1+0];
983 real c0011 = str->c[n+x1+1];
984 real c0012 = str->c[n+x1+2];
985 real c0013 = str->c[n+x1+3];
986 real c0014 = str->c[n+x1+4];
987 real c0015 = str->c[n+x1+5];
988 real c0016 = str->c[n+x1+6];
989 real c0017 = str->c[n+x1+7];
990
991 real c0100 = str->c[n+y1+0];
992 real c0101 = str->c[n+y1+1];
993 real c0102 = str->c[n+y1+2];
994 real c0103 = str->c[n+y1+3];
995 real c0104 = str->c[n+y1+4];
996 real c0105 = str->c[n+y1+5];
997 real c0106 = str->c[n+y1+6];
998 real c0107 = str->c[n+y1+7];
999
1000 real c1000 = str->c[n+z1+0];
1001 real c1001 = str->c[n+z1+1];
1002 real c1002 = str->c[n+z1+2];
1003 real c1003 = str->c[n+z1+3];
1004 real c1004 = str->c[n+z1+4];
1005 real c1005 = str->c[n+z1+5];
1006 real c1006 = str->c[n+z1+6];
1007 real c1007 = str->c[n+z1+7];
1008
1009 real c0110 = str->c[n+y1+x1+0];
1010 real c0111 = str->c[n+y1+x1+1];
1011 real c0112 = str->c[n+y1+x1+2];
1012 real c0113 = str->c[n+y1+x1+3];
1013 real c0114 = str->c[n+y1+x1+4];
1014 real c0115 = str->c[n+y1+x1+5];
1015 real c0116 = str->c[n+y1+x1+6];
1016 real c0117 = str->c[n+y1+x1+7];
1017
1018 real c1010 = str->c[n+z1+x1+0];
1019 real c1011 = str->c[n+z1+x1+1];
1020 real c1012 = str->c[n+z1+x1+2];
1021 real c1013 = str->c[n+z1+x1+3];
1022 real c1014 = str->c[n+z1+x1+4];
1023 real c1015 = str->c[n+z1+x1+5];
1024 real c1016 = str->c[n+z1+x1+6];
1025 real c1017 = str->c[n+z1+x1+7];
1026
1027 real c1100 = str->c[n+z1+y1+0];
1028 real c1101 = str->c[n+z1+y1+1];
1029 real c1102 = str->c[n+z1+y1+2];
1030 real c1103 = str->c[n+z1+y1+3];
1031 real c1104 = str->c[n+z1+y1+4];
1032 real c1105 = str->c[n+z1+y1+5];
1033 real c1106 = str->c[n+z1+y1+6];
1034 real c1107 = str->c[n+z1+y1+7];
1035
1036 real c1110 = str->c[n+z1+y1+x1+0];
1037 real c1111 = str->c[n+z1+y1+x1+1];
1038 real c1112 = str->c[n+z1+y1+x1+2];
1039 real c1113 = str->c[n+z1+y1+x1+3];
1040 real c1114 = str->c[n+z1+y1+x1+4];
1041 real c1115 = str->c[n+z1+y1+x1+5];
1042 real c1116 = str->c[n+z1+y1+x1+6];
1043 real c1117 = str->c[n+z1+y1+x1+7];
1044
1045 /* Evaluate spline values */
1046
1047 /* f */
1048 f_df[0] = (
1049 dzi*(
1050 dxi*(dyi*c0000+dy*c0100)
1051 +dx*(dyi*c0010+dy*c0110))
1052 +dz*(
1053 dxi*(dyi*c1000+dy*c1100)
1054 +dx*(dyi*c1010+dy*c1110)))
1055 +xg2/6*(
1056 dzi*(
1057 dxi3*(dyi*c0001+dy*c0101)
1058 +dx3*(dyi*c0011+dy*c0111))
1059 +dz*(
1060 dxi3*(dyi*c1001+dy*c1101)
1061 +dx3*(dyi*c1011+dy*c1111)))
1062 +yg2/6*(
1063 dzi*(
1064 dxi*(dyi3*c0002+dy3*c0102)
1065 +dx*(dyi3*c0012+dy3*c0112))
1066 +dz*(
1067 dxi*(dyi3*c1002+dy3*c1102)
1068 +dx*(dyi3*c1012+dy3*c1112)))
1069 +zg2/6*(
1070 dzi3*(
1071 dxi*(dyi*c0003+dy*c0103)
1072 +dx*(dyi*c0013+dy*c0113))
1073 +dz3*(
1074 dxi*(dyi*c1003+dy*c1103)
1075 +dx*(dyi*c1013+dy*c1113)))
1076 +xg2*yg2/36*(
1077 dzi*(
1078 dxi3*(dyi3*c0004+dy3*c0104)
1079 +dx3*(dyi3*c0014+dy3*c0114))
1080 +dz*(
1081 dxi3*(dyi3*c1004+dy3*c1104)
1082 +dx3*(dyi3*c1014+dy3*c1114)))
1083 +xg2*zg2/36*(
1084 dzi3*(
1085 dxi3*(dyi*c0005+dy*c0105)
1086 +dx3*(dyi*c0015+dy*c0115))
1087 +dz3*(
1088 dxi3*(dyi*c1005+dy*c1105)
1089 +dx3*(dyi*c1015+dy*c1115)))
1090 +yg2*zg2/36*(
1091 dzi3*(
1092 dxi*(dyi3*c0006+dy3*c0106)
1093 +dx*(dyi3*c0016+dy3*c0116))
1094 +dz3*(
1095 dxi*(dyi3*c1006+dy3*c1106)
1096 +dx*(dyi3*c1016+dy3*c1116)))
1097 +xg2*yg2*zg2/216*(
1098 dzi3*(
1099 dxi3*(dyi3*c0007+dy3*c0107)
1100 +dx3*(dyi3*c0017+dy3*c0117))
1101 +dz3*(
1102 dxi3*(dyi3*c1007+dy3*c1107)
1103 +dx3*(dyi3*c1017+dy3*c1117)));
1104
1105 /* df/dx */
1106 f_df[1] = xgi*(
1107 dzi*(
1108 -(dyi*c0000+dy*c0100)
1109 +(dyi*c0010+dy*c0110))
1110 +dz*(
1111 -(dyi*c1000+dy*c1100)
1112 +(dyi*c1010+dy*c1110)))
1113 +xg/6*(
1114 dzi*(
1115 dxi3dx*(dyi*c0001+dy*c0101)
1116 +dx3dx*(dyi*c0011+dy*c0111))
1117 +dz*(
1118 dxi3dx*(dyi*c1001 +dy*c1101)
1119 +dx3dx*(dyi*c1011+dy*c1111)))
1120 +xgi*yg2/6*(
1121 dzi*(
1122 -(dyi3*c0002+dy3*c0102)
1123 +(dyi3*c0012+dy3*c0112))
1124 +dz*(
1125 -(dyi3*c1002+dy3*c1102)
1126 +(dyi3*c1012+dy3*c1112)))
1127 +xgi*zg2/6*(
1128 dzi3*(
1129 -(dyi*c0003+dy*c0103)
1130 +(dyi*c0013+dy*c0113))
1131 +dz3*(
1132 -(dyi*c1003+dy*c1103)
1133 +(dyi*c1013+dy*c1113)))
1134 +xg*yg2/36*(
1135 dzi*(
1136 dxi3dx*(dyi3*c0004+dy3*c0104)
1137 +dx3dx*(dyi3*c0014+dy3*c0114))
1138 +dz*(
1139 dxi3dx*(dyi3*c1004+dy3*c1104)
1140 +dx3dx*(dyi3*c1014+dy3*c1114)))
1141 +xg*zg2/36*(
1142 dzi3*(
1143 dxi3dx*(dyi*c0005+dy*c0105)
1144 +dx3dx*(dyi*c0015+dy*c0115))
1145 +dz3*(
1146 dxi3dx*(dyi*c1005+dy*c1105)
1147 +dx3dx*(dyi*c1015+dy*c1115)))
1148 +xgi*yg2*zg2/36*(
1149 dzi3*(
1150 -(dyi3*c0006+dy3*c0106)
1151 +(dyi3*c0016+dy3*c0116))
1152 +dz3*(
1153 -(dyi3*c1006+dy3*c1106)
1154 +(dyi3*c1016+dy3*c1116)))
1155 +xg*yg2*zg2/216*(
1156 dzi3*(
1157 dxi3dx*(dyi3*c0007+dy3*c0107)
1158 +dx3dx*(dyi3*c0017+dy3*c0117))
1159 +dz3*(
1160 dxi3dx*(dyi3*c1007+dy3*c1107)
1161 +dx3dx*(dyi3*c1017+dy3*c1117)));
1162
1163 /* df/dy */
1164 f_df[2] = ygi*(
1165 dzi*(
1166 dxi*(-c0000+c0100)
1167 +dx*(-c0010+c0110))
1168 +dz*(
1169 dxi*(-c1000+c1100)
1170 +dx*(-c1010+c1110)))
1171 +ygi*xg2/6*(
1172 dzi*(
1173 dxi3*(-c0001+c0101)
1174 +dx3*(-c0011+c0111))
1175 +dz*(
1176 dxi3*(-c1001+c1101)
1177 +dx3*(-c1011+c1111)))
1178 +yg/6*(
1179 dzi*(
1180 dxi*(dyi3dy*c0002+dy3dy*c0102)
1181 +dx*(dyi3dy*c0012+dy3dy*c0112))
1182 +dz*(
1183 dxi*(dyi3dy*c1002+dy3dy*c1102)
1184 +dx*(dyi3dy*c1012+dy3dy*c1112)))
1185 +ygi*zg2/6*(
1186 dzi3*(
1187 dxi*(-c0003+c0103)
1188 +dx*(-c0013+c0113))
1189 +dz3*(
1190 dxi*(-c1003+c1103)
1191 +dx*(-c1013+c1113)))
1192 +xg2*yg/36*(
1193 dzi*(
1194 dxi3*(dyi3dy*c0004+dy3dy*c0104)
1195 +dx3*(dyi3dy*c0014+dy3dy*c0114))
1196 +dz*(
1197 dxi3*(dyi3dy*c1004+dy3dy*c1104)
1198 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
1199 +ygi*xg2*zg2/36*(
1200 dzi3*(
1201 dxi3*(-c0005+c0105)
1202 +dx3*(-c0015+c0115))
1203 +dz3*(
1204 dxi3*(-c1005+c1105)
1205 +dx3*(-c1015+c1115)))
1206 +yg*zg2/36*(
1207 dzi3*(
1208 dxi*(dyi3dy*c0006+dy3dy*c0106)
1209 +dx*(dyi3dy*c0016+dy3dy*c0116))
1210 +dz3*(
1211 dxi*(dyi3dy*c1006+dy3dy*c1106)
1212 +dx*(dyi3dy*c1016+dy3dy*c1116)))
1213 +xg2*yg*zg2/216*(
1214 dzi3*(
1215 dxi3*(dyi3dy*c0007+dy3dy*c0107)
1216 +dx3*(dyi3dy*c0017+dy3dy*c0117))
1217 +dz3*(
1218 dxi3*(dyi3dy*c1007+dy3dy*c1107)
1219 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
1220
1221 /* df/dz */
1222 f_df[3] = zgi*(
1223 -(
1224 dxi*(dyi*c0000+dy*c0100)
1225 +dx*(dyi*c0010+dy*c0110))
1226 +(
1227 dxi*(dyi*c1000+dy*c1100)
1228 +dx*(dyi*c1010+dy*c1110)))
1229 +xg2*zgi/6*(
1230 -(
1231 dxi3*(dyi*c0001+dy*c0101)
1232 +dx3*(dyi*c0011+dy*c0111))
1233 +(
1234 dxi3*(dyi*c1001+dy*c1101)
1235 +dx3*(dyi*c1011+dy*c1111)))
1236 +yg2*zgi/6*(
1237 -(
1238 dxi*(dyi3*c0002+dy3*c0102)
1239 +dx*(dyi3*c0012+dy3*c0112))
1240 +(
1241 dxi*(dyi3*c1002+dy3*c1102)
1242 +dx*(dyi3*c1012+dy3*c1112)))
1243 +zg/6*(
1244 dzi3dz*(
1245 dxi*(dyi*c0003+dy*c0103)
1246 +dx*(dyi*c0013+dy*c0113))
1247 +dz3dz*(
1248 dxi*(dyi*c1003+dy*c1103)
1249 +dx*(dyi*c1013+dy*c1113)))
1250 +xg2*yg2*zgi/36*(
1251 -(
1252 dxi3*(dyi3*c0004+dy3*c0104)
1253 +dx3*(dyi3*c0014+dy3*c0114))
1254 +(
1255 dxi3*(dyi3*c1004+dy3*c1104)
1256 +dx3*(dyi3*c1014+dy3*c1114)))
1257 +xg2*zg/36*(
1258 dzi3dz*(
1259 dxi3*(dyi*c0005+dy*c0105)
1260 +dx3*(dyi*c0015+dy*c0115))
1261 +dz3dz*(
1262 dxi3*(dyi*c1005+dy*c1105)
1263 +dx3*(dyi*c1015+dy*c1115)))
1264 +yg2*zg/36*(
1265 dzi3dz*(
1266 dxi*(dyi3*c0006+dy3*c0106)
1267 +dx*(dyi3*c0016+dy3*c0116))
1268 +dz3dz*(
1269 dxi*(dyi3*c1006+dy3*c1106)
1270 +dx*(dyi3*c1016+dy3*c1116)))
1271 +xg2*yg2*zg/216*(
1272 dzi3dz*(
1273 dxi3*(dyi3*c0007+dy3*c0107)
1274 +dx3*(dyi3*c0017+dy3*c0117))
1275 +dz3dz*(
1276 dxi3*(dyi3*c1007+dy3*c1107)
1277 +dx3*(dyi3*c1017+dy3*c1117)));
1278
1279 /* d2f/dx2 */
1280 f_df[4] = (
1281 dzi*(
1282 dxi*(dyi*c0001+dy*c0101)
1283 +dx*(dyi*c0011+dy*c0111))
1284 +dz*(
1285 dxi*(dyi*c1001+dy*c1101)
1286 +dx*(dyi*c1011+dy*c1111)))
1287 +yg2/6*(
1288 dzi*(
1289 dxi*(dyi3*c0004+dy3*c0104)
1290 +dx*(dyi3*c0014+dy3*c0114))
1291 +dz*(
1292 dxi*(dyi3*c1004+dy3*c1104)
1293 +dx*(dyi3*c1014+dy3*c1114)))
1294 +zg2/6*(
1295 dzi3*(
1296 dxi*(dyi*c0005+dy*c0105)
1297 +dx*(dyi*c0015+dy*c0115))
1298 +dz3*(
1299 dxi*(dyi*c1005+dy*c1105)
1300 +dx*(dyi*c1015+dy*c1115)))
1301 +yg2*zg2/36*(
1302 dzi3*(
1303 dxi*(dyi3*c0007+dy3*c0107)
1304 +dx*(dyi3*c0017+dy3*c0117))
1305 +dz3*(
1306 dxi*(dyi3*c1007+dy3*c1107)
1307 +dx*(dyi3*c1017+dy3*c1117)));
1308
1309 /* d2f/dy2 */
1310 f_df[5] = (
1311 dzi*(
1312 dxi*(dyi*c0002+dy*c0102)
1313 +dx*(dyi*c0012+dy*c0112))
1314 +dz*(
1315 dxi*(dyi*c1002+dy*c1102)
1316 +dx*(dyi*c1012+dy*c1112)))
1317 +xg2/6*(
1318 dzi*(
1319 dxi3*(dyi*c0004+dy*c0104)
1320 +dx3*(dyi*c0014+dy*c0114))
1321 +dz*(
1322 dxi3*(dyi*c1004+dy*c1104)
1323 +dx3*(dyi*c1014+dy*c1114)))
1324 +zg2/6*(
1325 dzi3*(
1326 dxi*(dyi*c0006+dy*c0106)
1327 +dx*(dyi*c0016+dy*c0116))
1328 +dz3*(
1329 dxi*(dyi*c1006+dy*c1106)
1330 +dx*(dyi*c1016+dy*c1116)))
1331 +xg2*zg2/36*(
1332 dzi3*(
1333 dxi3*(dyi*c0007+dy*c0107)
1334 +dx3*(dyi*c0017+dy*c0117))
1335 +dz3*(
1336 dxi3*(dyi*c1007+dy*c1107)
1337 +dx3*(dyi*c1017+dy*c1117)));
1338
1339 /* d2f/dz2 */
1340 f_df[6] = (
1341 dzi*(
1342 dxi*(dyi*c0003+dy*c0103)
1343 +dx*(dyi*c0013+dy*c0113))
1344 +dz*(
1345 dxi*(dyi*c1003+dy*c1103)
1346 +dx*(dyi*c1013+dy*c1113)))
1347 +xg2/6*(
1348 dzi*(
1349 dxi3*(dyi*c0005+dy*c0105)
1350 +dx3*(dyi*c0015+dy*c0115))
1351 +dz*(
1352 dxi3*(dyi*c1005+dy*c1105)
1353 +dx3*(dyi*c1015+dy*c1115)))
1354 +yg2/6*(
1355 dzi*(
1356 dxi*(dyi3*c0006+dy3*c0106)
1357 +dx*(dyi3*c0016+dy3*c0116))
1358 +dz*(
1359 dxi*(dyi3*c1006+dy3*c1106)
1360 +dx*(dyi3*c1016+dy3*c1116)))
1361 +xg2*yg2/36*(
1362 dzi*(
1363 dxi3*(dyi3*c0007+dy3*c0107)
1364 +dx3*(dyi3*c0017+dy3*c0117))
1365 +dz*(
1366 dxi3*(dyi3*c1007+dy3*c1107)
1367 +dx3*(dyi3*c1017+dy3*c1117)));
1368
1369 /* d2f/dxdy */
1370 f_df[7] = xgi*ygi*(
1371 dzi*(
1372 (c0000 -c0100)
1373 -(c0010-c0110))
1374 +dz*(
1375 (c1000 -c1100)
1376 -(c1010-c1110)))
1377 +ygi*xg/6*(
1378 dzi*(
1379 dxi3dx*(-c0001+c0101)
1380 +dx3dx*(-c0011+c0111))
1381 +dz*(
1382 dxi3dx*(-c1001+c1101)
1383 +dx3dx*(-c1011+c1111)))
1384 +xgi*yg/6*(
1385 dzi*(
1386 -(dyi3dy*c0002+dy3dy*c0102)
1387 +(dyi3dy*c0012+dy3dy*c0112))
1388 +dz*(
1389 -(dyi3dy*c1002+dy3dy*c1102)
1390 +(dyi3dy*c1012+dy3dy*c1112)))
1391 +xgi*ygi*zg2/6*(
1392 dzi3*(
1393 (c0003 -c0103)
1394 -(c0013-c0113))
1395 +dz3*(
1396 (c1003 -c1103)
1397 -(c1013-c1113)))
1398 +xg*yg/36*(
1399 dzi*(
1400 dxi3dx*(dyi3dy*c0004+dy3dy*c0104)
1401 +dx3dx*(dyi3dy*c0014+dy3dy*c0114))
1402 +dz*(
1403 dxi3dx*(dyi3dy*c1004+dy3dy*c1104)
1404 +dx3dx*(dyi3dy*c1014+dy3dy*c1114)))
1405 +ygi*xg*zg2/36*(
1406 dzi3*(
1407 dxi3dx*(-c0005+c0105)
1408 +dx3dx*(-c0015+c0115))
1409 +dz3*(
1410 dxi3dx*(-c1005+c1105)
1411 +dx3dx*(-c1015+c1115)))
1412 +xgi*yg*zg2/36*(
1413 dzi3*(
1414 -(dyi3dy*c0006+dy3dy*c0106)
1415 +(dyi3dy*c0016+dy3dy*c0116))
1416 +dz3*(
1417 -(dyi3dy*c1006+dy3dy*c1106)
1418 +(dyi3dy*c1016+dy3dy*c1116)))
1419 +xg*yg*zg2/216*(
1420 dzi3*(
1421 dxi3dx*(dyi3dy*c0007+dy3dy*c0107)
1422 +dx3dx*(dyi3dy*c0017+dy3dy*c0117))
1423 +dz3*(
1424 dxi3dx*(dyi3dy*c1007+dy3dy*c1107)
1425 +dx3dx*(dyi3dy*c1017+dy3dy*c1117)));
1426
1427 /* d2f/dxdz */
1428 f_df[8] = xgi*zgi*(
1429 (
1430 (dyi*c0000+dy*c0100)
1431 -(dyi*c0010+dy*c0110))
1432 -(
1433 (dyi*c1000+dy*c1100)
1434 -(dyi*c1010+dy*c1110)))
1435 +xg*zgi/6*(
1436 -(
1437 dxi3dx*(dyi*c0001+dy*c0101)
1438 +dx3dx*(dyi*c0011+dy*c0111))
1439 +(
1440 dxi3dx*(dyi*c1001+dy*c1101)
1441 +dx3dx*(dyi*c1011+dy*c1111)))
1442 +xgi*yg2*zgi/6*(
1443 (
1444 (dyi3*c0002+dy3*c0102)
1445 -(dyi3*c0012+dy3*c0112))
1446 -(
1447 (dyi3*c1002+dy3*c1102)
1448 -(dyi3*c1012+dy3*c1112)))
1449 +xgi*zg/6*(
1450 dzi3dz*(
1451 -(dyi*c0003+dy*c0103)
1452 +(dyi*c0013+dy*c0113))
1453 +dz3dz*(
1454 -(dyi*c1003+dy*c1103)
1455 +(dyi*c1013+dy*c1113)))
1456 +xg*yg2*zgi/36*(
1457 -(
1458 dxi3dx*(dyi3*c0004+dy3*c0104)
1459 +dx3dx*(dyi3*c0014+dy3*c0114))
1460 +(
1461 dxi3dx*(dyi3*c1004+dy3*c1104)
1462 +dx3dx*(dyi3*c1014+dy3*c1114)))
1463 +xg*zg/36*(
1464 dzi3dz*(
1465 dxi3dx*(dyi*c0005+dy*c0105)
1466 +dx3dx*(dyi*c0015+dy*c0115))
1467 +dz3dz*(
1468 dxi3dx*(dyi*c1005+dy*c1105)
1469 +dx3dx*(dyi*c1015+dy*c1115)))
1470 +xgi*yg2*zg/36*(
1471 dzi3dz*(
1472 -(dyi3*c0006+dy3*c0106)
1473 +(dyi3*c0016+dy3*c0116))
1474 +dz3dz*(
1475 -(dyi3*c1006+dy3*c1106)
1476 +(dyi3*c1016+dy3*c1116)))
1477 +xg*yg2*zg/216*(
1478 dzi3dz*(
1479 dxi3dx*(dyi3*c0007+dy3*c0107)
1480 +dx3dx*(dyi3*c0017+dy3*c0117))
1481 +dz3dz*(
1482 dxi3dx*(dyi3*c1007+dy3*c1107)
1483 +dx3dx*(dyi3*c1017+dy3*c1117)));
1484
1485 /* d2f/dydz */
1486 f_df[9] = ygi*zgi*(
1487 (
1488 dxi*(c0000 -c0100)
1489 +dx*(c0010-c0110))
1490 -(
1491 dxi*(c1000 -c1100)
1492 +dx*(c1010-c1110)))
1493 +ygi*xg2*zgi/6*(
1494 (
1495 dxi3*(c0001 -c0101)
1496 +dx3*(c0011-c0111))
1497 -(
1498 dxi3*(c1001 -c1101)
1499 +dx3*(c1011-c1111)))
1500 +yg*zgi/6*(
1501 -(
1502 dxi*(dyi3dy*c0002+dy3dy*c0102)
1503 +dx*(dyi3dy*c0012+dy3dy*c0112))
1504 +(
1505 dxi*(dyi3dy*c1002+dy3dy*c1102)
1506 +dx*(dyi3dy*c1012+dy3dy*c1112)))
1507 +ygi*zg/6*(
1508 dzi3dz*(
1509 dxi*(-c0003+c0103)
1510 +dx*(-c0013+c0113))
1511 +dz3dz*(
1512 dxi*(-c1003+c1103)
1513 +dx*(-c1013+c1113)))
1514 +xg2*yg*zgi/36*(
1515 -(
1516 dxi3*(dyi3dy*c0004+dy3dy*c0104)
1517 +dx3*(dyi3dy*c0014+dy3dy*c0114))
1518 +(
1519 dxi3*(dyi3dy*c1004+dy3dy*c1104)
1520 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
1521 +ygi*xg2*zg/36*(
1522 dzi3dz*(
1523 dxi3*(-c0005+c0105)
1524 +dx3*(-c0015+c0115))
1525 +dz3dz*(
1526 dxi3*(-c1005+c1105)
1527 +dx3*(-c1015+c1115)))
1528 +yg*zg/36*(
1529 dzi3dz*(
1530 dxi*(dyi3dy*c0006+dy3dy*c0106)
1531 +dx*(dyi3dy*c0016+dy3dy*c0116))
1532 +dz3dz*(
1533 dxi*(dyi3dy*c1006+dy3dy*c1106)
1534 +dx*(dyi3dy*c1016+dy3dy*c1116)))
1535 +xg2*yg*zg/216*(
1536 dzi3dz*(
1537 dxi3*(dyi3dy*c0007+dy3dy*c0107)
1538 +dx3*(dyi3dy*c0017+dy3dy*c0117))
1539 +dz3dz*(
1540 dxi3*(dyi3dy*c1007+dy3dy*c1107)
1541 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
1542 }
1543
1544 return err;
1545}
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
interp.h
Spline interpolation library.
interp3Dcomp_eval_f
DECLARE_TARGET_END a5err interp3Dcomp_eval_f(real *f, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D scalar field.
Definition interp3Dcomp.c:288
interp3Dcomp_setup
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.
Definition interp3Dcomp.c:259
interp3Dcomp_init_coeff
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.
Definition interp3Dcomp.c:37
NATURALBC
@ NATURALBC
Definition interp.h:37
PERIODICBC
@ PERIODICBC
Definition interp.h:38
interp3Dcomp_init_spline
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.
Definition interp3Dcomp.c:203
interp3Dcomp_eval_df
DECLARE_TARGET_END 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 1st derivatives.
Definition interp3Dcomp.c:448
interp3Dcomp_eval_ddf
DECLARE_TARGET_END a5err interp3Dcomp_eval_ddf(real *f_df, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D field and 1st and 2nd derivatives.
Definition interp3Dcomp.c:879
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.
splinecomp
void splinecomp(real *f, int n, int bc, real *c)
Calculate compact cubic spline interpolation coefficients in 1D.
Definition splinecomp.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 for ASCOT5 by doxygen 1.18.0