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
250a5err interp3Dcomp_eval_f(real* f, interp3D_data* str, real x, real y, real z) {
251
252 /* Make sure periodic coordinates are within [min, max] region. */
253 if(str->bc_x == PERIODICBC) {
254 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
255 x = x + (x < str->x_min) * (str->x_max - str->x_min);
256 }
257 if(str->bc_y == PERIODICBC) {
258 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
259 y = y + (y < str->y_min) * (str->y_max - str->y_min);
260 }
261 if(str->bc_z == PERIODICBC) {
262 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
263 z = z + (z < str->z_min) * (str->z_max - str->z_min);
264 }
265
266 /* Index for x variable. The -1 needed at exactly grid end. */
267 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
268 /* Normalized x coordinate in current cell */
269 real dx = (x - (str->x_min + i_x*str->x_grid)) / str->x_grid;
270 /* Helper varibles */
271 real dxi = 1.0 - dx;
272 real dx3 = dx*dx*dx - dx;
273 real dxi3 = (1.0 - dx) * (1.0 - dx) * (1.0 - dx) - (1.0 - dx);
274 real xg2 = str->x_grid*str->x_grid;
275
276 /* Index for y variable. The -1 needed at exactly grid end. */
277 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
278 /* Normalized y coordinate in current cell */
279 real dy = (y - (str->y_min + i_y*str->y_grid)) / str->y_grid;
280 /* Helper varibles */
281 real dyi = 1.0 - dy;
282 real dy3 = dy*dy*dy - dy;
283 real dyi3 = (1.0 - dy) * (1.0 - dy) * (1.0 - dy) - (1.0 - dy);
284 real yg2 = str->y_grid*str->y_grid;
285
286 /* Index for z variable. The -1 needed at exactly grid end. */
287 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
288 /* Normalized z coordinate in current cell */
289 real dz = (z - (str->z_min + i_z*str->z_grid)) / str->z_grid;
290 /* Helper varibles */
291 real dzi = 1.0 - dz;
292 real dz3 = dz*dz*dz - dz;
293 real dzi3 = (1.0 - dz) * (1.0 - dz) * (1.0 - dz) - (1.0-dz);
294 real zg2 = str->z_grid*str->z_grid;
295
297 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
298 int x1 = 8; /* Index jump one x forward */
299 int y1 = str->n_x*8; /* Index jump one y forward */
300 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
301
302 int err = 0;
303
304 /* Enforce periodic BC or check that the coordinate is within the domain. */
305 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
306 x1 = -(str->n_x-1)*x1;
307 }
308 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
309 err = 1;
310 }
311 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
312 y1 = -(str->n_y-1)*y1;
313 }
314 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
315 err = 1;
316 }
317 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
318 z1 = -(str->n_z-1)*z1;
319 }
320 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
321 err = 1;
322 }
323
324 if(!err) {
325
326 /* Evaluate spline value */
327 *f = (
328 dzi*(
329 dxi*(dyi*str->c[n+0]+dy*str->c[n+y1+0])
330 +dx*(dyi*str->c[n+x1+0]+dy*str->c[n+y1+x1+0]))
331 +dz*(
332 dxi*(dyi*str->c[n+z1+0]+dy*str->c[n+y1+z1+0])
333 +dx*(dyi*str->c[n+x1+z1+0]+dy*str->c[n+y1+z1+x1+0])))
334 +xg2/6*(
335 dzi*(
336 dxi3*(dyi*str->c[n+1]+dy*str->c[n+y1+1])
337 +dx3*(dyi*str->c[n+x1+1]+dy*str->c[n+y1+x1+1]))
338 +dz*(
339 dxi3*(dyi*str->c[n+z1+1]+dy*str->c[n+y1+z1+1])
340 +dx3*(dyi*str->c[n+x1+z1+1]+dy*str->c[n+y1+z1+x1+1])))
341 +yg2/6*(
342 dzi*(
343 dxi*(dyi3*str->c[n+2]+dy3*str->c[n+y1+2])
344 +dx*(dyi3*str->c[n+x1+2]+dy3*str->c[n+y1+x1+2]))
345 +dz*(
346 dxi*(dyi3*str->c[n+z1+2]+dy3*str->c[n+y1+z1+2])
347 +dx*(dyi3*str->c[n+x1+z1+2]+dy3*str->c[n+y1+z1+x1+2])))
348 +zg2/6*(
349 dzi3*(
350 dxi*(dyi*str->c[n+3]+dy*str->c[n+y1+3])
351 +dx*(dyi*str->c[n+x1+3]+dy*str->c[n+y1+x1+3]))
352 +dz3*(
353 dxi*(dyi*str->c[n+z1+3]+dy*str->c[n+y1+z1+3])
354 +dx*(dyi*str->c[n+x1+z1+3]+dy*str->c[n+y1+z1+x1+3])))
355 +xg2*yg2/36*(
356 dzi*(
357 dxi3*(dyi3*str->c[n+4]+dy3*str->c[n+y1+4])
358 +dx3*(dyi3*str->c[n+x1+4]+dy3*str->c[n+y1+x1+4]))
359 +dz*(
360 dxi3*(dyi3*str->c[n+z1+4]+dy3*str->c[n+y1+z1+4])
361 +dx3*(dyi3*str->c[n+x1+z1+4]+dy3*str->c[n+y1+z1+x1+4])))
362 +xg2*zg2/36*(
363 dzi3*(
364 dxi3*(dyi*str->c[n+5]+dy*str->c[n+y1+5])
365 +dx3*(dyi*str->c[n+x1+5]+dy*str->c[n+y1+x1+5]))
366 +dz3*(
367 dxi3*(dyi*str->c[n+z1+5]+dy*str->c[n+y1+z1+5])
368 +dx3*(dyi*str->c[n+x1+z1+5]+dy*str->c[n+y1+z1+x1+5])))
369 +yg2*zg2/36*(
370 dzi3*(
371 dxi*(dyi3*str->c[n+6]+dy3*str->c[n+y1+6])
372 +dx*(dyi3*str->c[n+x1+6]+dy3*str->c[n+y1+x1+6]))
373 +dz3*(
374 dxi*(dyi3*str->c[n+z1+6]+dy3*str->c[n+y1+z1+6])
375 +dx*(dyi3*str->c[n+x1+z1+6]+dy3*str->c[n+y1+z1+x1+6])))
376 +xg2*yg2*zg2/216*(
377 dzi3*(
378 dxi3*(dyi3*str->c[n+7]+dy3*str->c[n+y1+7])
379 +dx3*(dyi3*str->c[n+x1+7]+dy3*str->c[n+y1+x1+7]))
380 +dz3*(
381 dxi3*(dyi3*str->c[n+z1+7]+dy3*str->c[n+y1+z1+7])
382 +dx3*(dyi3*str->c[n+x1+z1+7]+dy3*str->c[n+y1+z1+x1+7])));
383
384 }
385
386 return err;
387}
388
416a5err interp3Dcomp_eval_df(real* f_df, interp3D_data* str,
417 real x, real y, real z) {
418
419 /* Make sure periodic coordinates are within [min, max] region. */
420 if(str->bc_x == PERIODICBC) {
421 x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
422 x = x + (x < str->x_min) * (str->x_max - str->x_min);
423 }
424 if(str->bc_y == PERIODICBC) {
425 y = fmod(y - str->y_min, str->y_max - str->y_min) + str->y_min;
426 y = y + (y < str->y_min) * (str->y_max - str->y_min);
427 }
428 if(str->bc_z == PERIODICBC) {
429 z = fmod(z - str->z_min, str->z_max - str->z_min) + str->z_min;
430 z = z + (z < str->z_min) * (str->z_max - str->z_min);
431 }
432
433 /* Index for x variable. The -1 needed at exactly grid end. */
434 int i_x = (x - str->x_min) / str->x_grid - 1*(x==str->x_max);
435 /* Normalized x coordinate in current cell */
436 real dx = ( x - (str->x_min + i_x*str->x_grid) ) / str->x_grid;
437 /* Helper variables */
438 real dx3 = dx*dx*dx - dx;
439 real dx3dx = 3*dx*dx - 1.0;
440 real dxi = 1.0 - dx;
441 real dxi3 = dxi*dxi*dxi - dxi;
442 real dxi3dx = -3*dxi*dxi + 1.0;
443 real xg = str->x_grid;
444 real xg2 = xg*xg;
445 real xgi = 1.0 / xg;
446
447 /* Index for y variable. The -1 needed at exactly grid end. */
448 int i_y = (y - str->y_min) / str->y_grid - 1*(y==str->y_max);
449 /* Normalized y coordinate in current cell */
450 real dy = ( y - (str->y_min + i_y*str->y_grid) ) / str->y_grid;
451 /* Helper variables */
452 real dy3 = dy*dy*dy-dy;
453 real dy3dy = 3*dy*dy - 1.0;
454 real dyi = 1.0 - dy;
455 real dyi3 = dyi*dyi*dyi - dyi;
456 real dyi3dy = -3*dyi*dyi + 1.0;
457 real yg = str->y_grid;
458 real yg2 = yg*yg;
459 real ygi = 1.0 / yg;
460
461 /* Index for z variable. The -1 needed at exactly grid end. */
462 int i_z = (z - str->z_min) / str->z_grid - 1*(z==str->z_max);
463 /* Normalized z coordinate in current cell */
464 real dz = ( z - (str->z_min + i_z*str->z_grid) ) / str->z_grid;
465 /* Helper variables */
466 real dz3 = dz*dz*dz - dz;
467 real dz3dz = 3*dz*dz - 1.0;
468 real dzi = 1.0 - dz;
469 real dzi3 = dzi*dzi*dzi - dzi;
470 real dzi3dz = -3*dzi*dzi + 1.0;
471 real zg = str->z_grid;
472 real zg2 = zg*zg;
473 real zgi = 1.0 / zg;
474
475 /* Index jump to cell */
476 int n = i_z*str->n_y*str->n_x*8 + i_y*str->n_x*8 + i_x*8;
477 int x1 = 8; /* Index jump one x forward */
478 int y1 = str->n_x*8; /* Index jump one y forward */
479 int z1 = str->n_y*str->n_x*8; /* Index jump one z forward */
480
481 int err = 0;
482
483 /* Enforce periodic BC or check that the coordinate is within the domain. */
484 if( str->bc_x == PERIODICBC && i_x == str->n_x-1 ) {
485 x1 = -(str->n_x-1)*x1;
486 }
487 else if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
488 err = 1;
489 }
490 if( str->bc_y == PERIODICBC && i_y == str->n_y-1 ) {
491 y1 = -(str->n_y-1)*y1;
492 }
493 else if( str->bc_y == NATURALBC && !(y >= str->y_min && y <= str->y_max) ) {
494 err = 1;
495 }
496 if( str->bc_z == PERIODICBC && i_z == str->n_z-1 ) {
497 z1 = -(str->n_z-1)*z1;
498 }
499 else if( str->bc_z == NATURALBC && !(z >= str->z_min && z <= str->z_max) ) {
500 err = 1;
501 }
502
503 if(!err) {
504
505 /* Fetch coefficients explicitly to fetch those that are adjacent
506 subsequently and to store in temporary variables coefficients that
507 will be used multiple times. This is to decrease computational time,
508 by exploiting simultaneous extraction of adjacent memory and by
509 avoiding going through the long str->c array repeatedly. */
510 real c0000 = str->c[n+0];
511 real c0001 = str->c[n+1];
512 real c0002 = str->c[n+2];
513 real c0003 = str->c[n+3];
514 real c0004 = str->c[n+4];
515 real c0005 = str->c[n+5];
516 real c0006 = str->c[n+6];
517 real c0007 = str->c[n+7];
518
519 real c0010 = str->c[n+x1+0];
520 real c0011 = str->c[n+x1+1];
521 real c0012 = str->c[n+x1+2];
522 real c0013 = str->c[n+x1+3];
523 real c0014 = str->c[n+x1+4];
524 real c0015 = str->c[n+x1+5];
525 real c0016 = str->c[n+x1+6];
526 real c0017 = str->c[n+x1+7];
527
528 real c0100 = str->c[n+y1+0];
529 real c0101 = str->c[n+y1+1];
530 real c0102 = str->c[n+y1+2];
531 real c0103 = str->c[n+y1+3];
532 real c0104 = str->c[n+y1+4];
533 real c0105 = str->c[n+y1+5];
534 real c0106 = str->c[n+y1+6];
535 real c0107 = str->c[n+y1+7];
536
537 real c1000 = str->c[n+z1+0];
538 real c1001 = str->c[n+z1+1];
539 real c1002 = str->c[n+z1+2];
540 real c1003 = str->c[n+z1+3];
541 real c1004 = str->c[n+z1+4];
542 real c1005 = str->c[n+z1+5];
543 real c1006 = str->c[n+z1+6];
544 real c1007 = str->c[n+z1+7];
545
546 real c0110 = str->c[n+y1+x1+0];
547 real c0111 = str->c[n+y1+x1+1];
548 real c0112 = str->c[n+y1+x1+2];
549 real c0113 = str->c[n+y1+x1+3];
550 real c0114 = str->c[n+y1+x1+4];
551 real c0115 = str->c[n+y1+x1+5];
552 real c0116 = str->c[n+y1+x1+6];
553 real c0117 = str->c[n+y1+x1+7];
554
555 real c1010 = str->c[n+z1+x1+0];
556 real c1011 = str->c[n+z1+x1+1];
557 real c1012 = str->c[n+z1+x1+2];
558 real c1013 = str->c[n+z1+x1+3];
559 real c1014 = str->c[n+z1+x1+4];
560 real c1015 = str->c[n+z1+x1+5];
561 real c1016 = str->c[n+z1+x1+6];
562 real c1017 = str->c[n+z1+x1+7];
563
564 real c1100 = str->c[n+z1+y1+0];
565 real c1101 = str->c[n+z1+y1+1];
566 real c1102 = str->c[n+z1+y1+2];
567 real c1103 = str->c[n+z1+y1+3];
568 real c1104 = str->c[n+z1+y1+4];
569 real c1105 = str->c[n+z1+y1+5];
570 real c1106 = str->c[n+z1+y1+6];
571 real c1107 = str->c[n+z1+y1+7];
572
573 real c1110 = str->c[n+z1+y1+x1+0];
574 real c1111 = str->c[n+z1+y1+x1+1];
575 real c1112 = str->c[n+z1+y1+x1+2];
576 real c1113 = str->c[n+z1+y1+x1+3];
577 real c1114 = str->c[n+z1+y1+x1+4];
578 real c1115 = str->c[n+z1+y1+x1+5];
579 real c1116 = str->c[n+z1+y1+x1+6];
580 real c1117 = str->c[n+z1+y1+x1+7];
581
582 /* Evaluate spline values */
583
584 /* f */
585 f_df[0] = (
586 dzi*(
587 dxi*(dyi*c0000+dy*c0100)
588 +dx*(dyi*c0010+dy*c0110))
589 +dz*(
590 dxi*(dyi*c1000+dy*c1100)
591 +dx*(dyi*c1010+dy*c1110)))
592 +xg2/6*(
593 dzi*(
594 dxi3*(dyi*c0001+dy*c0101)
595 +dx3*(dyi*c0011+dy*c0111))
596 +dz*(
597 dxi3*(dyi*c1001+dy*c1101)
598 +dx3*(dyi*c1011+dy*c1111)))
599 +yg2/6*(
600 dzi*(
601 dxi*(dyi3*c0002+dy3*c0102)
602 +dx*(dyi3*c0012+dy3*c0112))
603 +dz*(
604 dxi*(dyi3*c1002+dy3*c1102)
605 +dx*(dyi3*c1012+dy3*c1112)))
606 +zg2/6*(
607 dzi3*(
608 dxi*(dyi*c0003+dy*c0103)
609 +dx*(dyi*c0013+dy*c0113))
610 +dz3*(
611 dxi*(dyi*c1003+dy*c1103)
612 +dx*(dyi*c1013+dy*c1113)))
613 +xg2*yg2/36*(
614 dzi*(
615 dxi3*(dyi3*c0004+dy3*c0104)
616 +dx3*(dyi3*c0014+dy3*c0114))
617 +dz*(
618 dxi3*(dyi3*c1004+dy3*c1104)
619 +dx3*(dyi3*c1014+dy3*c1114)))
620 +xg2*zg2/36*(
621 dzi3*(
622 dxi3*(dyi*c0005+dy*c0105)
623 +dx3*(dyi*c0015+dy*c0115))
624 +dz3*(
625 dxi3*(dyi*c1005+dy*c1105)
626 +dx3*(dyi*c1015+dy*c1115)))
627 +yg2*zg2/36*(
628 dzi3*(
629 dxi*(dyi3*c0006+dy3*c0106)
630 +dx*(dyi3*c0016+dy3*c0116))
631 +dz3*(
632 dxi*(dyi3*c1006+dy3*c1106)
633 +dx*(dyi3*c1016+dy3*c1116)))
634 +xg2*yg2*zg2/216*(
635 dzi3*(
636 dxi3*(dyi3*c0007+dy3*c0107)
637 +dx3*(dyi3*c0017+dy3*c0117))
638 +dz3*(
639 dxi3*(dyi3*c1007+dy3*c1107)
640 +dx3*(dyi3*c1017+dy3*c1117)));
641
642 /* df/dx */
643 f_df[1] = xgi*(
644 dzi*(
645 -(dyi*c0000+dy*c0100)
646 +(dyi*c0010+dy*c0110))
647 +dz*(
648 -(dyi*c1000+dy*c1100)
649 +(dyi*c1010+dy*c1110)))
650 +xg/6*(
651 dzi*(
652 dxi3dx*(dyi*c0001+dy*c0101)
653 +dx3dx*(dyi*c0011+dy*c0111))
654 +dz*(
655 dxi3dx*(dyi*c1001 +dy*c1101)
656 +dx3dx*(dyi*c1011+dy*c1111)))
657 +xgi*yg2/6*(
658 dzi*(
659 -(dyi3*c0002+dy3*c0102)
660 +(dyi3*c0012+dy3*c0112))
661 +dz*(
662 -(dyi3*c1002+dy3*c1102)
663 +(dyi3*c1012+dy3*c1112)))
664 +xgi*zg2/6*(
665 dzi3*(
666 -(dyi*c0003+dy*c0103)
667 +(dyi*c0013+dy*c0113))
668 +dz3*(
669 -(dyi*c1003+dy*c1103)
670 +(dyi*c1013+dy*c1113)))
671 +xg*yg2/36*(
672 dzi*(
673 dxi3dx*(dyi3*c0004+dy3*c0104)
674 +dx3dx*(dyi3*c0014+dy3*c0114))
675 +dz*(
676 dxi3dx*(dyi3*c1004+dy3*c1104)
677 +dx3dx*(dyi3*c1014+dy3*c1114)))
678 +xg*zg2/36*(
679 dzi3*(
680 dxi3dx*(dyi*c0005+dy*c0105)
681 +dx3dx*(dyi*c0015+dy*c0115))
682 +dz3*(
683 dxi3dx*(dyi*c1005+dy*c1105)
684 +dx3dx*(dyi*c1015+dy*c1115)))
685 +xgi*yg2*zg2/36*(
686 dzi3*(
687 -(dyi3*c0006+dy3*c0106)
688 +(dyi3*c0016+dy3*c0116))
689 +dz3*(
690 -(dyi3*c1006+dy3*c1106)
691 +(dyi3*c1016+dy3*c1116)))
692 +xg*yg2*zg2/216*(
693 dzi3*(
694 dxi3dx*(dyi3*c0007+dy3*c0107)
695 +dx3dx*(dyi3*c0017+dy3*c0117))
696 +dz3*(
697 dxi3dx*(dyi3*c1007+dy3*c1107)
698 +dx3dx*(dyi3*c1017+dy3*c1117)));
699
700 /* df/dy */
701 f_df[2] = ygi*(
702 dzi*(
703 dxi*(-c0000+c0100)
704 +dx*(-c0010+c0110))
705 +dz*(
706 dxi*(-c1000+c1100)
707 +dx*(-c1010+c1110)))
708 +ygi*xg2/6*(
709 dzi*(
710 dxi3*(-c0001+c0101)
711 +dx3*(-c0011+c0111))
712 +dz*(
713 dxi3*(-c1001+c1101)
714 +dx3*(-c1011+c1111)))
715 +yg/6*(
716 dzi*(
717 dxi*(dyi3dy*c0002+dy3dy*c0102)
718 +dx*(dyi3dy*c0012+dy3dy*c0112))
719 +dz*(
720 dxi*(dyi3dy*c1002+dy3dy*c1102)
721 +dx*(dyi3dy*c1012+dy3dy*c1112)))
722 +ygi*zg2/6*(
723 dzi3*(
724 dxi*(-c0003+c0103)
725 +dx*(-c0013+c0113))
726 +dz3*(
727 dxi*(-c1003+c1103)
728 +dx*(-c1013+c1113)))
729 +xg2*yg/36*(
730 dzi*(
731 dxi3*(dyi3dy*c0004+dy3dy*c0104)
732 +dx3*(dyi3dy*c0014+dy3dy*c0114))
733 +dz*(
734 dxi3*(dyi3dy*c1004+dy3dy*c1104)
735 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
736 +ygi*xg2*zg2/36*(
737 dzi3*(
738 dxi3*(-c0005+c0105)
739 +dx3*(-c0015+c0115))
740 +dz3*(
741 dxi3*(-c1005+c1105)
742 +dx3*(-c1015+c1115)))
743 +yg*zg2/36*(
744 dzi3*(
745 dxi*(dyi3dy*c0006+dy3dy*c0106)
746 +dx*(dyi3dy*c0016+dy3dy*c0116))
747 +dz3*(
748 dxi*(dyi3dy*c1006+dy3dy*c1106)
749 +dx*(dyi3dy*c1016+dy3dy*c1116)))
750 +xg2*yg*zg2/216*(
751 dzi3*(
752 dxi3*(dyi3dy*c0007+dy3dy*c0107)
753 +dx3*(dyi3dy*c0017+dy3dy*c0117))
754 +dz3*(
755 dxi3*(dyi3dy*c1007+dy3dy*c1107)
756 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
757
758 /* df/dz */
759 f_df[3] = zgi*(
760 -(
761 dxi*(dyi*c0000+dy*c0100)
762 +dx*(dyi*c0010+dy*c0110))
763 +(
764 dxi*(dyi*c1000+dy*c1100)
765 +dx*(dyi*c1010+dy*c1110)))
766 +xg2*zgi/6*(
767 -(
768 dxi3*(dyi*c0001+dy*c0101)
769 +dx3*(dyi*c0011+dy*c0111))
770 +(
771 dxi3*(dyi*c1001+dy*c1101)
772 +dx3*(dyi*c1011+dy*c1111)))
773 +yg2*zgi/6*(
774 -(
775 dxi*(dyi3*c0002+dy3*c0102)
776 +dx*(dyi3*c0012+dy3*c0112))
777 +(
778 dxi*(dyi3*c1002+dy3*c1102)
779 +dx*(dyi3*c1012+dy3*c1112)))
780 +zg/6*(
781 dzi3dz*(
782 dxi*(dyi*c0003+dy*c0103)
783 +dx*(dyi*c0013+dy*c0113))
784 +dz3dz*(
785 dxi*(dyi*c1003+dy*c1103)
786 +dx*(dyi*c1013+dy*c1113)))
787 +xg2*yg2*zgi/36*(
788 -(
789 dxi3*(dyi3*c0004+dy3*c0104)
790 +dx3*(dyi3*c0014+dy3*c0114))
791 +(
792 dxi3*(dyi3*c1004+dy3*c1104)
793 +dx3*(dyi3*c1014+dy3*c1114)))
794 +xg2*zg/36*(
795 dzi3dz*(
796 dxi3*(dyi*c0005+dy*c0105)
797 +dx3*(dyi*c0015+dy*c0115))
798 +dz3dz*(
799 dxi3*(dyi*c1005+dy*c1105)
800 +dx3*(dyi*c1015+dy*c1115)))
801 +yg2*zg/36*(
802 dzi3dz*(
803 dxi*(dyi3*c0006+dy3*c0106)
804 +dx*(dyi3*c0016+dy3*c0116))
805 +dz3dz*(
806 dxi*(dyi3*c1006+dy3*c1106)
807 +dx*(dyi3*c1016+dy3*c1116)))
808 +xg2*yg2*zg/216*(
809 dzi3dz*(
810 dxi3*(dyi3*c0007+dy3*c0107)
811 +dx3*(dyi3*c0017+dy3*c0117))
812 +dz3dz*(
813 dxi3*(dyi3*c1007+dy3*c1107)
814 +dx3*(dyi3*c1017+dy3*c1117)));
815
816 /* d2f/dx2 */
817 f_df[4] = (
818 dzi*(
819 dxi*(dyi*c0001+dy*c0101)
820 +dx*(dyi*c0011+dy*c0111))
821 +dz*(
822 dxi*(dyi*c1001+dy*c1101)
823 +dx*(dyi*c1011+dy*c1111)))
824 +yg2/6*(
825 dzi*(
826 dxi*(dyi3*c0004+dy3*c0104)
827 +dx*(dyi3*c0014+dy3*c0114))
828 +dz*(
829 dxi*(dyi3*c1004+dy3*c1104)
830 +dx*(dyi3*c1014+dy3*c1114)))
831 +zg2/6*(
832 dzi3*(
833 dxi*(dyi*c0005+dy*c0105)
834 +dx*(dyi*c0015+dy*c0115))
835 +dz3*(
836 dxi*(dyi*c1005+dy*c1105)
837 +dx*(dyi*c1015+dy*c1115)))
838 +yg2*zg2/36*(
839 dzi3*(
840 dxi*(dyi3*c0007+dy3*c0107)
841 +dx*(dyi3*c0017+dy3*c0117))
842 +dz3*(
843 dxi*(dyi3*c1007+dy3*c1107)
844 +dx*(dyi3*c1017+dy3*c1117)));
845
846 /* d2f/dy2 */
847 f_df[5] = (
848 dzi*(
849 dxi*(dyi*c0002+dy*c0102)
850 +dx*(dyi*c0012+dy*c0112))
851 +dz*(
852 dxi*(dyi*c1002+dy*c1102)
853 +dx*(dyi*c1012+dy*c1112)))
854 +xg2/6*(
855 dzi*(
856 dxi3*(dyi*c0004+dy*c0104)
857 +dx3*(dyi*c0014+dy*c0114))
858 +dz*(
859 dxi3*(dyi*c1004+dy*c1104)
860 +dx3*(dyi*c1014+dy*c1114)))
861 +zg2/6*(
862 dzi3*(
863 dxi*(dyi*c0006+dy*c0106)
864 +dx*(dyi*c0016+dy*c0116))
865 +dz3*(
866 dxi*(dyi*c1006+dy*c1106)
867 +dx*(dyi*c1016+dy*c1116)))
868 +xg2*zg2/36*(
869 dzi3*(
870 dxi3*(dyi*c0007+dy*c0107)
871 +dx3*(dyi*c0017+dy*c0117))
872 +dz3*(
873 dxi3*(dyi*c1007+dy*c1107)
874 +dx3*(dyi*c1017+dy*c1117)));
875
876 /* d2f/dz2 */
877 f_df[6] = (
878 dzi*(
879 dxi*(dyi*c0003+dy*c0103)
880 +dx*(dyi*c0013+dy*c0113))
881 +dz*(
882 dxi*(dyi*c1003+dy*c1103)
883 +dx*(dyi*c1013+dy*c1113)))
884 +xg2/6*(
885 dzi*(
886 dxi3*(dyi*c0005+dy*c0105)
887 +dx3*(dyi*c0015+dy*c0115))
888 +dz*(
889 dxi3*(dyi*c1005+dy*c1105)
890 +dx3*(dyi*c1015+dy*c1115)))
891 +yg2/6*(
892 dzi*(
893 dxi*(dyi3*c0006+dy3*c0106)
894 +dx*(dyi3*c0016+dy3*c0116))
895 +dz*(
896 dxi*(dyi3*c1006+dy3*c1106)
897 +dx*(dyi3*c1016+dy3*c1116)))
898 +xg2*yg2/36*(
899 dzi*(
900 dxi3*(dyi3*c0007+dy3*c0107)
901 +dx3*(dyi3*c0017+dy3*c0117))
902 +dz*(
903 dxi3*(dyi3*c1007+dy3*c1107)
904 +dx3*(dyi3*c1017+dy3*c1117)));
905
906 /* d2f/dxdy */
907 f_df[7] = xgi*ygi*(
908 dzi*(
909 (c0000 -c0100)
910 -(c0010-c0110))
911 +dz*(
912 (c1000 -c1100)
913 -(c1010-c1110)))
914 +ygi*xg/6*(
915 dzi*(
916 dxi3dx*(-c0001+c0101)
917 +dx3dx*(-c0011+c0111))
918 +dz*(
919 dxi3dx*(-c1001+c1101)
920 +dx3dx*(-c1011+c1111)))
921 +xgi*yg/6*(
922 dzi*(
923 -(dyi3dy*c0002+dy3dy*c0102)
924 +(dyi3dy*c0012+dy3dy*c0112))
925 +dz*(
926 -(dyi3dy*c1002+dy3dy*c1102)
927 +(dyi3dy*c1012+dy3dy*c1112)))
928 +xgi*ygi*zg2/6*(
929 dzi3*(
930 (c0003 -c0103)
931 -(c0013-c0113))
932 +dz3*(
933 (c1003 -c1103)
934 -(c1013-c1113)))
935 +xg*yg/36*(
936 dzi*(
937 dxi3dx*(dyi3dy*c0004+dy3dy*c0104)
938 +dx3dx*(dyi3dy*c0014+dy3dy*c0114))
939 +dz*(
940 dxi3dx*(dyi3dy*c1004+dy3dy*c1104)
941 +dx3dx*(dyi3dy*c1014+dy3dy*c1114)))
942 +ygi*xg*zg2/36*(
943 dzi3*(
944 dxi3dx*(-c0005+c0105)
945 +dx3dx*(-c0015+c0115))
946 +dz3*(
947 dxi3dx*(-c1005+c1105)
948 +dx3dx*(-c1015+c1115)))
949 +xgi*yg*zg2/36*(
950 dzi3*(
951 -(dyi3dy*c0006+dy3dy*c0106)
952 +(dyi3dy*c0016+dy3dy*c0116))
953 +dz3*(
954 -(dyi3dy*c1006+dy3dy*c1106)
955 +(dyi3dy*c1016+dy3dy*c1116)))
956 +xg*yg*zg2/216*(
957 dzi3*(
958 dxi3dx*(dyi3dy*c0007+dy3dy*c0107)
959 +dx3dx*(dyi3dy*c0017+dy3dy*c0117))
960 +dz3*(
961 dxi3dx*(dyi3dy*c1007+dy3dy*c1107)
962 +dx3dx*(dyi3dy*c1017+dy3dy*c1117)));
963
964 /* d2f/dxdz */
965 f_df[8] = xgi*zgi*(
966 (
967 (dyi*c0000+dy*c0100)
968 -(dyi*c0010+dy*c0110))
969 -(
970 (dyi*c1000+dy*c1100)
971 -(dyi*c1010+dy*c1110)))
972 +xg*zgi/6*(
973 -(
974 dxi3dx*(dyi*c0001+dy*c0101)
975 +dx3dx*(dyi*c0011+dy*c0111))
976 +(
977 dxi3dx*(dyi*c1001+dy*c1101)
978 +dx3dx*(dyi*c1011+dy*c1111)))
979 +xgi*yg2*zgi/6*(
980 (
981 (dyi3*c0002+dy3*c0102)
982 -(dyi3*c0012+dy3*c0112))
983 -(
984 (dyi3*c1002+dy3*c1102)
985 -(dyi3*c1012+dy3*c1112)))
986 +xgi*zg/6*(
987 dzi3dz*(
988 -(dyi*c0003+dy*c0103)
989 +(dyi*c0013+dy*c0113))
990 +dz3dz*(
991 -(dyi*c1003+dy*c1103)
992 +(dyi*c1013+dy*c1113)))
993 +xg*yg2*zgi/36*(
994 -(
995 dxi3dx*(dyi3*c0004+dy3*c0104)
996 +dx3dx*(dyi3*c0014+dy3*c0114))
997 +(
998 dxi3dx*(dyi3*c1004+dy3*c1104)
999 +dx3dx*(dyi3*c1014+dy3*c1114)))
1000 +xg*zg/36*(
1001 dzi3dz*(
1002 dxi3dx*(dyi*c0005+dy*c0105)
1003 +dx3dx*(dyi*c0015+dy*c0115))
1004 +dz3dz*(
1005 dxi3dx*(dyi*c1005+dy*c1105)
1006 +dx3dx*(dyi*c1015+dy*c1115)))
1007 +xgi*yg2*zg/36*(
1008 dzi3dz*(
1009 -(dyi3*c0006+dy3*c0106)
1010 +(dyi3*c0016+dy3*c0116))
1011 +dz3dz*(
1012 -(dyi3*c1006+dy3*c1106)
1013 +(dyi3*c1016+dy3*c1116)))
1014 +xg*yg2*zg/216*(
1015 dzi3dz*(
1016 dxi3dx*(dyi3*c0007+dy3*c0107)
1017 +dx3dx*(dyi3*c0017+dy3*c0117))
1018 +dz3dz*(
1019 dxi3dx*(dyi3*c1007+dy3*c1107)
1020 +dx3dx*(dyi3*c1017+dy3*c1117)));
1021
1022 /* d2f/dydz */
1023 f_df[9] = ygi*zgi*(
1024 (
1025 dxi*(c0000 -c0100)
1026 +dx*(c0010-c0110))
1027 -(
1028 dxi*(c1000 -c1100)
1029 +dx*(c1010-c1110)))
1030 +ygi*xg2*zgi/6*(
1031 (
1032 dxi3*(c0001 -c0101)
1033 +dx3*(c0011-c0111))
1034 -(
1035 dxi3*(c1001 -c1101)
1036 +dx3*(c1011-c1111)))
1037 +yg*zgi/6*(
1038 -(
1039 dxi*(dyi3dy*c0002+dy3dy*c0102)
1040 +dx*(dyi3dy*c0012+dy3dy*c0112))
1041 +(
1042 dxi*(dyi3dy*c1002+dy3dy*c1102)
1043 +dx*(dyi3dy*c1012+dy3dy*c1112)))
1044 +ygi*zg/6*(
1045 dzi3dz*(
1046 dxi*(-c0003+c0103)
1047 +dx*(-c0013+c0113))
1048 +dz3dz*(
1049 dxi*(-c1003+c1103)
1050 +dx*(-c1013+c1113)))
1051 +xg2*yg*zgi/36*(
1052 -(
1053 dxi3*(dyi3dy*c0004+dy3dy*c0104)
1054 +dx3*(dyi3dy*c0014+dy3dy*c0114))
1055 +(
1056 dxi3*(dyi3dy*c1004+dy3dy*c1104)
1057 +dx3*(dyi3dy*c1014+dy3dy*c1114)))
1058 +ygi*xg2*zg/36*(
1059 dzi3dz*(
1060 dxi3*(-c0005+c0105)
1061 +dx3*(-c0015+c0115))
1062 +dz3dz*(
1063 dxi3*(-c1005+c1105)
1064 +dx3*(-c1015+c1115)))
1065 +yg*zg/36*(
1066 dzi3dz*(
1067 dxi*(dyi3dy*c0006+dy3dy*c0106)
1068 +dx*(dyi3dy*c0016+dy3dy*c0116))
1069 +dz3dz*(
1070 dxi*(dyi3dy*c1006+dy3dy*c1106)
1071 +dx*(dyi3dy*c1016+dy3dy*c1116)))
1072 +xg2*yg*zg/216*(
1073 dzi3dz*(
1074 dxi3*(dyi3dy*c0007+dy3dy*c0107)
1075 +dx3*(dyi3dy*c0017+dy3dy*c0117))
1076 +dz3dz*(
1077 dxi3*(dyi3dy*c1007+dy3dy*c1107)
1078 +dx3*(dyi3dy*c1017+dy3dy*c1117)));
1079 }
1080
1081 return err;
1082}
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
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
interp3Dcomp_eval_df
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.
Definition interp3Dcomp.c:416
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_f
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:250
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.
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 Wed Sep 25 2024 07:10:18 for ASCOT5 by doxygen 1.10.0