36 #define ID_MAX_NVAR 256
68 int i, ip, nv, success;
70 char *sub_str, sline[256];
71 const char delimiters[] =
" \t\r\f\n";
76 print1 (
"> Input data:\n\n");
84 fp = fopen(grid_fname,
"r");
86 print1 (
"! InputDataSet: grid file %s not found\n",grid_fname);
92 sub_str = strtok(sline, delimiters);
93 while (sub_str != NULL){
94 if (!strcmp(sub_str,
"GEOMETRY:")) {
95 sub_str = strtok(NULL, delimiters);
99 sub_str = strtok(NULL, delimiters);
108 print1 (
"! InputDataSet: unknown geometry\n");
112 print1 (
" Input grid file: %s\n", grid_fname);
113 print1 (
" Input grid geometry: %s\n", sub_str);
122 fgetpos(fp, &file_pos);
124 if (sline[0] !=
'#') success = 1;
127 fsetpos(fp, &file_pos);
135 fscanf (fp,
"%d \n",&
id_nx1);
137 for (i = 0; i <
id_nx1; i++){
138 fscanf(fp,
"%d %lf %lf\n", &ip, &xl, &xr);
142 fscanf (fp,
"%d \n",&
id_nx2);
144 for (i = 0; i <
id_nx2; i++){
145 fscanf(fp,
"%d %lf %lf\n", &ip, &xl, &xr);
149 fscanf (fp,
"%d \n",&
id_nx3);
151 for (i = 0; i <
id_nx3; i++){
152 fscanf(fp,
"%d %lf %lf\n", &ip, &xl, &xr);
159 if (id_nx1 == 1)
id_x1[0] = 0.0;
160 if (id_nx2 == 1)
id_x2[0] = 0.0;
161 if (id_nx3 == 1)
id_x3[0] = 0.0;
163 print1 (
" Input grid extension: x1 = [%12.3e, %12.3e] (%d points)\n",
165 print1 (
"\t\t\t x2 = [%12.3e, %12.3e] (%d points)\n",
167 print1 (
"\t\t\t x3 = [%12.3e, %12.3e] (%d points)\n",
180 if (get_var[nv] != -1) {
207 int i,
j,
k, nv, swap_endian=
NO;
208 size_t dsize, dcount;
223 print1 (
" Input data file: %s (endianity: %s) \n",
224 data_fname, endianity);
230 dcount = strlen(data_fname);
231 for (i = 0; i < 3; i++) ext[i] = data_fname[dcount-3+i];
233 if (!strcmp(ext,
"dbl")){
234 print1 (
" Precision: (double)\n");
235 dsize =
sizeof(double);
236 }
else if (!strcmp(ext,
"flt")) {
237 print1 (
" Precision:\t\t (single)\n");
238 dsize =
sizeof(float);
240 print1 (
"! InputDataRead: unsupported data type '%s'\n",ext);
248 fp = fopen(data_fname,
"rb");
250 print1 (
"! InputDataRead: file %s does not exist\n");
253 for (nv = 0; nv <
id_nvar; nv++){
258 if (dsize ==
sizeof(
double)){
259 for (k = 0; k <
id_nx3; k++){
260 for (j = 0; j <
id_nx2; j++){
261 for (i = 0; i <
id_nx1; i++){
262 if (fread (&udbl, dsize, dcount, fp) != dcount){
263 print1 (
"! InputDataRead: error reading data %d.\n",nv);
270 for (k = 0; k <
id_nx3; k++){
271 for (j = 0; j <
id_nx2; j++){
272 for (i = 0; i <
id_nx1; i++){
273 if (fread (&uflt, dsize, dcount, fp) != dcount){
274 print1 (
"! InputDataRead: error reading data %d.\n",nv);
301 int il = 0, jl = 0, kl = 0;
304 int i,
j,
k, nv, inv;
312 #if GEOMETRY == CARTESIAN
319 R = sqrt(x1*x1 + x2*x2);
324 x1 = R; x2 = z; x3 = phi;
327 R = sqrt(x1*x1 + x2*x2);
332 x1 = R; x2 = phi; x3 = z;
334 double r, theta, phi;
335 r =
D_EXPAND(x1*x1, + x2*x2, + x3*x3);
340 if (theta < 0.0) theta += 2.0*
CONST_PI;
342 x1 = r; x2 = theta; x3 = phi;
344 print1 (
"! InputDataInterpolate: invalid or unsupported coordinate transformation.\n");
347 #elif GEOMETRY == CYLINDRICAL
353 double r, theta, phi;
354 r =
D_EXPAND(x1*x1, + x2*x2, + 0.0);
358 if (theta < 0.0) theta += 2.0*
CONST_PI;
360 x1 = r; x2 = theta; x3 = phi;
362 print1 (
"! InputDataInterpolate: invalid or unsupported coordinate transformation.\n");
365 #elif GEOMETRY == POLAR
376 x1 = x; x2 = y; x3 = z;
378 print1 (
"! InputDataInterpolate: invalid or unsupported coordinate transformation.\n");
381 #elif GEOMETRY == SPHERICAL
389 x = x1*sin(x2)*cos(x3);
390 y = x1*sin(x2)*sin(x3);
393 x1 = x; x2 = y; x3 = z;
395 print1 (
"! InputDataInterpolate: invalid or unsupported coordinate transformation.\n");
422 while (il != (ih-1)){
424 if (x1 <=
id_x1[im]) ih = im;
431 while (jl != (jh-1)){
433 if (x2 <=
id_x2[jm]) jh = jm;
441 while (kl != (kh - 1)){
443 if (x3 <=
id_x3[km]) kh = km;
465 for (nv = 0; nv <
id_nvar; nv++) {
468 vs[inv] = V[kl][jl][il]*(1.0 - xx)*(1.0 - yy)*(1.0 - zz)
469 + V[kl][jl][il+1]*xx*(1.0 - yy)*(1.0 - zz);
471 vs[inv] += V[kl][jl+1][il]*(1.0 - xx)*yy*(1.0 - zz)
472 + V[kl][jl+1][il+1]*xx*yy*(1.0 - zz);
475 vs[inv] += V[kl+1][jl][il]*(1.0 - xx)*(1.0 - yy)*zz
476 + V[kl+1][jl][il+1]*xx*(1.0 - yy)*zz
477 + V[kl+1][jl+1][il]*(1.0 - xx)*yy*zz
478 + V[kl+1][jl+1][il+1]*xx*yy*zz;
491 for (nv = 0; nv <
id_nvar; nv++){
492 free ((
char *)
Vin[nv][0][0]);
493 free ((
char *)
Vin[nv][0]);
494 free ((
char *)
Vin[nv]);
void print1(const char *fmt,...)
#define ARRAY_3D(nx, ny, nz, type)
#define ARRAY_1D(nx, type)
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
#define QUIT_PLUTO(e_code)