HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
rescale_interp.c
Go to the documentation of this file.
1 
9 #include "decs.h"
10 
11 
19 int rescale(int which, int dir, FTYPE *pr, struct of_geom *ptrgeom,FTYPE *p2interp)
20 {
21  FTYPE scale[NPR2INTERP],r,th,X[NDIM],V[NDIM];
22  int ii,jj,kk;
23  int pl,pliter;
24  struct of_state q;
25  extern int quasivsq_compute(FTYPE *pr, struct of_geom *geom, FTYPE *quasivsq);
26  FTYPE quasivsq;
27  extern int limit_quasivsq(FTYPE quasivsqnew, struct of_geom *geom, FTYPE *pr);
28  FTYPE vcon[NDIM],ucon[NDIM];
29  int j;
30  FTYPE newpr[NPR];
31  FTYPE uconrel[NDIM],ucovrel[NDIM];
32  FTYPE vconrel[NDIM];
33  FTYPE normuconrel,normuconrel_fromui;
34  FTYPE gamma,qsq;
35 
36 #if(VARTOINTERPFIELD==PULSARFIELD)
37  extern void getconsts(FTYPE *uconmetp, FTYPE *V, struct of_geom *ptrgeom, FTYPE (*dxdxp)[NDIM],FTYPE *uconconst);
38  extern void undoconsts(FTYPE *uconconst, FTYPE *V, struct of_geom *ptrgeom, FTYPE (*dxdxp)[NDIM],FTYPE *uconmetp);
39 #endif
40  FTYPE Bconin[NDIM],Bconout[NDIM];
41  FTYPE Uconin[NDIM],Uconout[NDIM];
42  FTYPE dxdxp[NDIM][NDIM];
43 
44 
45  int dopl[NPR2INTERP];
46  for(pl=0;pl<NPR2INTERP;pl++) dopl[pl]=0; // default don't do
47  PINTERPLOOP(pliter,pl) dopl[pl]=1; // only do if in loop range (do this instead of looping over actual rescalings below)
48 
49  ii=ptrgeom->i;
50  jj=ptrgeom->j;
51  kk=ptrgeom->k;
52 
53  coord_ijk(ii,jj,kk,ptrgeom->p,X);
54  bl_coord_ijk(ii,jj,kk,ptrgeom->p,V);
55  r=V[1]; th=V[2];
56 
57 #if(VARTOINTERP==PRIMTOINTERP)
58 
59  // dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if RESCALEINTERP=%d\n",VARTOINTERP,RESCALEINTERP);
60  // myexit(1);
61 
62  // allow multiple types of rescales, and so just identity if not doing this (i.e. PRIMTOINTERP)
63  if(which==DORESCALE){ // rescale before interpolation
64  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
65  }
66  else if(which==UNRESCALE){ // unrescale after interpolation
67  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
68  }
69  else{
70  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
71  myexit(100);
72  }
73 
74 
75 
76 
77 #elif(VARTOINTERP==PRIMTOINTERP_JONRESCALED1)
78 
79  if(dir==1){
80  // optimized for pole
81  if(rescaletype==0){
82  scale[RHO]=pow(r,-1.5);
83  }
84  else if(rescaletype==1){
85  scale[RHO]=pow(r,-2.7);
86  }
87  scale[UU]=scale[RHO]/r;
88  scale[U1]=scale[RHO];
89  scale[U2]=1.0;
90  scale[U3]=1.0/(r*r);
91  if(rescaletype==0){
92  scale[B1]=scale[U3];
93  }
94  else if(rescaletype==1){
95  scale[B1]=pow(r,-2.4);
96  }
97  // if(statpos[2]+jj < 0 || startpos[2]+jj >= totalsize[2]) scale[B1] *= -1. ;
98  scale[B2]=scale[B1];
99  scale[B3]=scale[B1];
100 
101  if(DOENTROPY) scale[ENTROPY]=1.0;
102  }
103  else if(dir==2){
104  scale[RHO]=1.0;
105  scale[UU]=1.0;
106  scale[U1]=1.0;
107  scale[U2]=1.0;
108  scale[U3]=1.0;
109  scale[B1]=1.0;
110  scale[B2]=1.0;
111  scale[B3]=1.0;
112  if(DOENTROPY) scale[ENTROPY]=1.0;
113  }
114  else{
115  dualfprintf(fail_file,"rescale(): no such direction! dir=%d\n",dir);
116  myexit(100);
117  }
118 
119 
120  if(which==DORESCALE){ // rescale before interpolation
121  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl]/scale[pl];
122  }
123  else if(which==UNRESCALE){ // unrescale after interpolation
124  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl]*scale[pl];
125  }
126  else{
127  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
128  myexit(100);
129  }
130 
131 
132 
133 #elif(VARTOINTERP==CONSTOINTERP)
134  // this doesn't work at all, even if no bug.
135  // doesn't work even if setup better guess, as in interpU code.
136 
137 
138  if(which==DORESCALE){ // rescale before interpolation
139  MYFUN(get_state(pr, ptrgeom, &q),"interp.c:rescale()", "get_state()", 1);
140  MYFUN(primtoU(UDIAG,pr, &q, ptrgeom, p2interp, NULL),"interp.c:rescale()", "primtoU()", 1);
141  }
142  else if(which==UNRESCALE){ // unrescale after interpolation
143  struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
144  int showmessages=1;
145  int allowlocalfailurefixandnoreport=1;
146  int eomtype=EOMDEFAULT;
147  FTYPE dissmeasure=-1.0; // assume energy inversion ok
148  int whichcap=CAPTYPEBASIC;
149  int whichmethod=MODEDEFAULT;
150  int modprim=0;
151  int checkoninversiongas=CHECKONINVERSION;
152  int checkoninversionrad=CHECKONINVERSIONRAD;
153  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,OTHERUTOPRIM,UDIAG,p2interp, NULL, ptrgeom, dissmeasure, pr, pr,&newtonstats),"interp.c:rescale()", "Utoprimgen", 1);
154  }
155  else{
156  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
157  myexit(100);
158  }
159 
160 
161 #elif(VARTOINTERP==PRIMTOINTERPLGDEN)
162 
163 
164  if(which==DORESCALE){ // rescale before interpolation
165  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
166  if(dopl[RHO]) p2interp[RHO]=log(pr[RHO]);
167  if(dopl[UU]) p2interp[UU]=log(pr[UU]);
168  }
169  else if(which==UNRESCALE){ // unrescale after interpolation
170  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
171  if(dopl[RHO]) pr[RHO]=exp(p2interp[RHO]);
172  if(dopl[UU]) pr[UU]=exp(p2interp[UU]);
173  }
174  else{
175  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
176  myexit(100);
177  }
178 
179 #elif(VARTOINTERP==PRIMTOINTERP_LGDEN_RHOU)
180  // unstable!
181  // probably because using LGDEN with RHOU.
182  // Between 1E-6 and 1, density would be 1E-3 in log.
183  // But v=0 and 1, so rho*v~0 and 1, but then resulting v=10^3 !
184 
185  if(which==DORESCALE){ // rescale before interpolation
186  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
187  if(dopl[RHO]) p2interp[RHO]=log(pr[RHO]);
188  if(dopl[UU]) p2interp[UU]=log(pr[UU]);
189 
190  if(dopl[U1]) p2interp[U1]=pr[RHO]*pr[U1];
191  if(dopl[U2]) p2interp[U2]=pr[RHO]*pr[U2];
192  if(dopl[U3]) p2interp[U3]=pr[RHO]*pr[U3];
193  }
194  else if(which==UNRESCALE){ // unrescale after interpolation
195  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
196  if(dopl[RHO]) pr[RHO]=exp(p2interp[RHO]);
197  if(dopl[UU]) pr[UU]=exp(p2interp[UU]);
198 
199  if(dopl[U1]) pr[U1]=p2interp[U1]/pr[RHO];
200  if(dopl[U2]) pr[U2]=p2interp[U2]/pr[RHO];
201  if(dopl[U3]) pr[U3]=p2interp[U3]/pr[RHO];
202 
203  }
204  else{
205  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
206  myexit(100);
207  }
208 
209 
210 
211 #elif(VARTOINTERP==PRIMTOINTERP_RHOU)
212  // kinda works, except for test=7 (peak problem)
213 
214  if(which==DORESCALE){ // rescale before interpolation
215  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
216 
217  if(dopl[U1]) p2interp[U1]=pr[RHO]*pr[U1];
218  if(dopl[U2]) p2interp[U2]=pr[RHO]*pr[U2];
219  if(dopl[U3]) p2interp[U3]=pr[RHO]*pr[U3];
220  }
221  else if(which==UNRESCALE){ // unrescale after interpolation
222  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
223 
224  if(dopl[U1]) pr[U1]=p2interp[U1]/pr[RHO];
225  if(dopl[U2]) pr[U2]=p2interp[U2]/pr[RHO];
226  if(dopl[U3]) pr[U3]=p2interp[U3]/pr[RHO];
227 
228  }
229  else{
230  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
231  myexit(100);
232  }
233 
234 
235 
236 
237 #elif(VARTOINTERP==PRIMTOINTERP_VSQ)
238 
239 
240 
241 #if(DOEXTRAINTERP==0)
242  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
243  myexit(1);
244 #endif
245 
246 
247  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
248  if(dopl[U1]){
249  quasivsq_compute(pr, ptrgeom, &quasivsq);
250  }
251  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
252 
253  if(dopl[U1]) p2interp[U1]=pr[RHO]*pr[U1];
254  if(dopl[U2]) p2interp[U2]=pr[RHO]*pr[U2];
255  if(dopl[U3]) p2interp[U3]=pr[RHO]*pr[U3];
256 
257  if(dopl[VSQ]){
258  // max helps limit oscillatory behavior for non-limiter schemes
259  p2interp[VSQ]=max(quasivsq,0.0);
260  //p2interp[VSQ]=log(quasivsq); // assumes positive
261  }
262  }
263  else if(which==UNRESCALE){ // after interpolation
264  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
265 
266  if(dopl[U1]) pr[U1]=p2interp[U1]/pr[RHO];
267  if(dopl[U2]) pr[U2]=p2interp[U2]/pr[RHO];
268  if(dopl[U3]) pr[U3]=p2interp[U3]/pr[RHO];
269 
270  if(dopl[VSQ]){
271  // now rescale velocities to agree with quasivsq
272  // max helps limit oscillatory behavior for non-limiter schemes
273  limit_quasivsq(max(p2interp[VSQ],0.0),ptrgeom,pr);
274  // limit_quasivsq(exp(p2interp[VSQ]),ptrgeom,pr);
275  }
276  }
277  else{
278  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
279  myexit(100);
280  }
281 
282 
283 
284 #elif(VARTOINTERP==PRIMTOINTERP_3VEL_GAMMA)
285 
286 
287 #if(DOEXTRAINTERP==0)
288  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
289  myexit(1);
290 #endif
291 
292 
293  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
294 
295  if(dopl[U1] && VSQ>=0){
296  pr2ucon(WHICHVEL,pr, ptrgeom ,ucon);
297  // 3-velocity
298  SLOOPA(j) vcon[j]=ucon[j]/ucon[TT];
299  }
300 
301  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
302 
303  for(pl=U1;pl<=U3;pl++) if(dopl[pl]) p2interp[pl]= vcon[pl-U1+1];
304 
305  if(dopl[VSQ] && VSQ>=0) p2interp[VSQ]=ucon[TT];
306 
307  }
308  else if(which==UNRESCALE){ // after interpolation
309 
310  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
311 
312  for(pl=U1;pl<=U3;pl++) if(dopl[pl]) vcon[pl-U1+1] = p2interp[pl];
313 
314  if(dopl[U1] && VSQ>=0){
315  SLOOPA(j) ucon[j] = vcon[j]*p2interp[VSQ];
316  ucon2pr(WHICHVEL,ucon,ptrgeom,pr);
317  }
318 
319  }
320  else{
321  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
322  myexit(100);
323  }
324 
325 
326 
327 #elif(VARTOINTERP==PRIMTOINTERP_RHOV_GAMMA)
328 
329 #if(DOEXTRAINTERP==0)
330  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
331  myexit(1);
332 #endif
333 
334 
335  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
336  if(dopl[U1] && dopl[VSQ] && VSQ>=0){
337  pr2ucon(WHICHVEL,pr, ptrgeom ,ucon);
338  // 3-velocity
339  SLOOPA(j) vcon[j]=ucon[j]/ucon[TT];
340  }
341 
342  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
343 
344  //comment this out if you do not want to interpolate rho * v
345  for(pl=U1;pl<=U3;pl++) if(dopl[pl]) p2interp[pl]= pr[RHO] * vcon[pl-U1+1];
346 
347  if(dopl[VSQ]) p2interp[VSQ]=ucon[TT];
348 
349  }
350  else if(which==UNRESCALE){ // after interpolation
351 
352  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
353 
354  //comment this out if you do notwant to interp. rhov
355  for(pl=U1;pl<=U3;pl++) if(dopl[pl]) vcon[pl-U1+1] = p2interp[pl]/pr[RHO];
356 
357  if(dopl[U1] && dopl[VSQ] && VSQ>=0){
358  SLOOPA(j) ucon[j] = vcon[j]*p2interp[VSQ];
359  ucon2pr(WHICHVEL,ucon,ptrgeom,pr);
360  }
361 
362  }
363  else{
364  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
365  myexit(100);
366  }
367 
368 
369 #elif(VARTOINTERP==PRIMTOINTERP_VELREL4SQ)
370 
371 
372 #if(DOEXTRAINTERP==0)
373  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
374  myexit(1);
375 #endif
376 
377 #if(WHICHVEL==VEL3)
378  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if WHICHVEL=%d\n",VARTOINTERP,WHICHVEL);
379  myexit(1);
380 #endif
381 
382 
383 
384  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
385 
386  if(dopl[U1]){
387  // get relative 4-velocity
388  if(WHICHVEL!=VELREL4){
389  pr2ucon(WHICHVEL,pr, ptrgeom ,ucon);
390  ucon2pr(VELREL4,ucon,ptrgeom,newpr);
391  uconrel[TT]=0.0;
392  SLOOPA(j) uconrel[j]=newpr[UU+j];
393  }
394  else{
395  uconrel[TT]=0.0;
396  SLOOPA(j) uconrel[j]=pr[UU+j];
397  }
398 
399  // get |\tilde{u}|
400  lower_vec(uconrel,ptrgeom,ucovrel);
401  //normuconrel=sqrt(dot(uconrel,ucovrel));
402  normuconrel=dot(uconrel,ucovrel);
403  }
404 
405  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
406  if(dopl[VSQ]) p2interp[VSQ]=normuconrel;
407 
408  }
409  else if(which==UNRESCALE){ // after interpolation
410 
411  // assign over everything, adjust velocity below
412  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
413 
414  if(dopl[U1]){
415  // renormalize by relative 4-velocity
416  if(WHICHVEL!=VELREL4){
417  // get relative 4-velocity from interpolated velocity
418  pr2ucon(WHICHVEL,p2interp, ptrgeom ,ucon);
419  ucon2pr(VELREL4,ucon,ptrgeom,newpr);
420  uconrel[TT]=0.0;
421  SLOOPA(j) uconrel[j]=newpr[UU+j];
422  }
423  else{
424  uconrel[TT]=0.0;
425  SLOOPA(j) uconrel[j]=p2interp[UU+j];
426  }
427 
428  // get |\tilde{u}| from interpolated \tilde{u}^i
429  lower_vec(uconrel,ptrgeom,ucovrel);
430  // normuconrel_fromui=sqrt(dot(uconrel,ucovrel));
431  normuconrel_fromui=dot(uconrel,ucovrel);
432 
433  }
434 
435  if(dopl[VSQ] && VSQ>0){
436  if(WHICHVEL==VEL3){
437  // pr2ucon(WHICHVEL,p2interp,ptrgeom,ucon);
438  // ucon2pr(VELREL4,ucon,ptrgeom,newpr); // now have relative 4-velocity
439  // uconrel[TT]=0.0;
440  // SLOOPA(j) uconrel[j]=newpr[UU+j];
441  // not done, but should be possible
442  myexit(1);
443  }
444  else{ // WHICHVEL==VEL4 or WHICHVEL==VELREL4 : acts same on both
445  // directly renormalize primitives
446  // SLOOPA(j) pr[UU+j] = p2interp[UU+j]*p2interp[VSQ]/absuconrel_fromui;
447  SLOOPA(j) pr[UU+j] = p2interp[UU+j]*fabs(p2interp[VSQ]/normuconrel_fromui);
448  // done!
449  }
450  }
451 
452 
453 
454 
455 
456  }
457  else{
458  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
459  myexit(100);
460  }
461 
462 
463 
464 #elif(VARTOINTERP==PRIMTOINTERP_3VELREL_GAMMAREL)
465 
466 
467 #if(DOEXTRAINTERP==0)
468  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
469  myexit(1);
470 #endif
471 
472 
473  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
474 
475  if(dopl[U1]){
476  // get relative 4-velocity
477  if(WHICHVEL!=VELREL4){
478  pr2ucon(WHICHVEL,pr, ptrgeom ,ucon);
479  ucon2pr(VELREL4,ucon,ptrgeom,newpr);
480  uconrel[TT]=0.0;
481  SLOOPA(j) uconrel[j]=newpr[UU+j];
482  }
483  else{
484  uconrel[TT]=0.0;
485  SLOOPA(j) uconrel[j]=pr[UU+j];
486  }
487 
488  // get Lorentz factor w.r.t. relative 4-velocity
489  gamma_calc_fromuconrel(uconrel,ptrgeom,&gamma,&qsq);
490  }
491 
492  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
493 
494  if(dopl[U1]){
495  // interpolate relative 3-velocity
496  for(pl=U1;pl<=U3;pl++) p2interp[pl]= uconrel[pl-U1+1]/gamma;
497  }
498 
499  if(dopl[VSQ] && VSQ>=0){
500  // interpolate \gamma separately
501  p2interp[VSQ]=gamma;
502  }
503 
504  }
505  else if(which==UNRESCALE){ // after interpolation
506 
507  // assign over everything, adjust velocity below
508  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
509 
510  if(dopl[U1] && dopl[VSQ] && VSQ>=0){
511  // get relative 4-velocity from \gamma and relative 3-velocity
512  uconrel[TT]=0;
513  SLOOPA(j) uconrel[j]=p2interp[UU+j]*p2interp[VSQ];
514 
515  // get WHICHVEL velocity
516  if(WHICHVEL!=VELREL4){
517  pr2ucon(VELREL4,p2interp, ptrgeom ,ucon);
518  ucon2pr(WHICHVEL,ucon,ptrgeom,newpr);
519  SLOOPA(j) pr[UU+j]=newpr[UU+j];
520  }
521  else{
522  SLOOPA(j) pr[UU+j]=uconrel[j];
523  }
524 
525  }
526 
527 
528  }
529  else{
530  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
531  myexit(100);
532  }
533 
534 
535 
536 #elif(VARTOINTERP==PRIMTOINTERP_3VELREL_GAMMAREL_DXDXP)
537 
538 
539 #if(DOEXTRAINTERP==0)
540  dualfprintf(fail_file,"Shouldn't be trying to do VARTOINTERP=%d if DOEXTRAINTERP=%d\n",VARTOINTERP,DOEXTRAINTERP);
541  myexit(1);
542 #endif
543 
544  //compute jacobian
545  //dxdxprim(X, V, dxdxp);
546  dxdxprim_ijk( ptrgeom->i, ptrgeom->j, ptrgeom->k, ptrgeom->p, dxdxp );
547 
548  if(which==DORESCALE){ // before interpolation, get quantities to interpolate
549 
550  if(dopl[U1]){
551  // get relative 4-velocity
552  if(WHICHVEL!=VELREL4){
553  pr2ucon(WHICHVEL,pr, ptrgeom ,ucon);
554  ucon2pr(VELREL4,ucon,ptrgeom,newpr);
555  uconrel[TT]=0.0;
556  SLOOPA(j) uconrel[j]=newpr[UU+j];
557  }
558  else{
559  uconrel[TT]=0.0;
560  SLOOPA(j) uconrel[j]=pr[UU+j];
561  }
562 
563  // get Lorentz factor w.r.t. relative 4-velocity
564  gamma_calc_fromuconrel(uconrel,ptrgeom,&gamma,&qsq);
565  }
566 
567  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
568 
569  if(dopl[U1]){
570  // interpolate relative 3-velocity
571  for(pl=U1;pl<=U3;pl++) p2interp[pl]= uconrel[pl-U1+1]/gamma;
572  }
573  if(dopl[U2]){
574  //Interpolate u^\theta and B^\theta instead of u^2 and B^2:
575  p2interp[U2] *= dxdxp[2][2];
576  }
577  if(dopl[B2]) p2interp[B2] *= dxdxp[2][2];
578 
579  if(dopl[VSQ]&&VSQ>=0){
580  // interpolate \gamma separately
581  p2interp[VSQ]=gamma;
582  }
583 
584  }
585  else if(which==UNRESCALE){ // after interpolation
586 
587  //return back to u^2 and B^2 from u^\theta and B^\theta
588  if(dopl[U2]) p2interp[U2] /= dxdxp[2][2];
589  if(dopl[B2]) p2interp[B2] /= dxdxp[2][2];
590 
591  // assign over everything, adjust velocity below
592  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
593 
594  if(dopl[U1]){
595  // get relative 4-velocity from \gamma and relative 3-velocity
596  uconrel[TT]=0;
597  SLOOPA(j) uconrel[j]=p2interp[UU+j]*p2interp[VSQ];
598 
599 
600  // get WHICHVEL velocity
601  if(WHICHVEL!=VELREL4){
602  pr2ucon(VELREL4,p2interp, ptrgeom ,ucon);
603  ucon2pr(WHICHVEL,ucon,ptrgeom,newpr);
604  SLOOPA(j) pr[UU+j]=newpr[UU+j];
605  }
606  else{
607  SLOOPA(j) pr[UU+j]=uconrel[j];
608  }
609  }
610 
611 
612 
613  }
614  else{
615  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
616  myexit(100);
617  }
618 
619 
620 
621 #elif(VARTOINTERP==PRIMTOINTERP_RAMESH1)
622 
623 
624  if(which==DORESCALE){ // rescale before interpolation
625  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
626 
627  // for the fields, do ramesh-like interpolation
628  if(dopl[B1]){ pl=B1; p2interp[pl] = pr[pl] * (sqrt(ptrgeom->gcov[GIND(1,1)])*pow(r,2.0-nu) ); }
629  if(dopl[B2]){ pl=B2; p2interp[pl] = pr[pl] * (sqrt(ptrgeom->gcov[GIND(2,2)])*pow(r,2.0-nu) ); }
630  if(dopl[B3]){ pl=B3; p2interp[pl] = pr[pl] * (sqrt(ptrgeom->gcov[GIND(3,3)])*pow(r,2.0-nu) ); }
631  }
632  else if(which==UNRESCALE){ // unrescale after interpolation
633  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
634 
635  // for the fields, do ramesh-like interpolation
636  if(dopl[B1]){ pl=B1; pr[pl] = p2interp[pl] / (sqrt(ptrgeom->gcov[GIND(1,1)])*pow(r,2.0-nu) ); }
637  if(dopl[B2]){ pl=B2; pr[pl] = p2interp[pl] / (sqrt(ptrgeom->gcov[GIND(2,2)])*pow(r,2.0-nu) ); }
638  if(dopl[B3]){ pl=B3; pr[pl] = p2interp[pl] / (sqrt(ptrgeom->gcov[GIND(3,3)])*pow(r,2.0-nu) ); }
639 
640  }
641  else{
642  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
643  myexit(100);
644  }
645 
646 
647 #elif(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION)
648 
649 
650  if(which==DORESCALE){ // rescale before interpolation
651 
652  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
653 
654  // if(dopl[U1]) p2interp[U1]=pr[U1]*(ptrgeom->gdet);
655  // if(dopl[U2]) p2interp[U2]=pr[U2]*(ptrgeom->gdet);
656  // if(dopl[U3]) p2interp[U3]=pr[U3]*(ptrgeom->gdet);
657 
658  if(dir==1||dir==3){
659  if(dopl[U1]) p2interp[U1]=pr[U1]*(ptrgeom->gdet);
660  if(dopl[U2]) p2interp[U2]=pr[U2]*(ptrgeom->gdet);
661  if(dopl[U3]) p2interp[U3]=pr[U3]*(ptrgeom->gdet);
662 
663  if(URAD1>=0){
664  if(dopl[URAD1]) p2interp[URAD1]=pr[URAD1]*(ptrgeom->gdet);
665  if(dopl[URAD2]) p2interp[URAD2]=pr[URAD2]*(ptrgeom->gdet);
666  if(dopl[URAD3]) p2interp[URAD3]=pr[URAD3]*(ptrgeom->gdet);
667  }
668 
669  if(dopl[B1]) p2interp[B1]=pr[B1]*(ptrgeom->gdet);
670  if(dopl[B2]) p2interp[B2]=pr[B2]*(ptrgeom->gdet);
671  if(dopl[B3]) p2interp[B3]=pr[B3]*(ptrgeom->gdet);
672  }
673  if(0&&dir==2){
674  //if(dopl[U1]) p2interp[U1]=pr[U1]*(ptrgeom->gdet);
675  //if(dopl[U2]) p2interp[U2]=pr[U2]*(ptrgeom->gdet);
676  //if(dopl[U3]) p2interp[U3]=pr[U3]*(ptrgeom->gdet);
677 
678  if(dopl[B1]) p2interp[B1]=pr[B1]*(ptrgeom->gdet); //sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3);
679  if(dopl[B2]) p2interp[B2]=pr[B2]*(ptrgeom->gdet);//*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*pow(V[1],3);
680  if(dopl[B3]) p2interp[B3]=pr[B3]*(ptrgeom->gdet);//*pow(V[1],3);
681  }
682  }
683  else if(which==UNRESCALE){ // unrescale after interpolation
684 
685  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
686 
687  //if(dopl[U1]) pr[U1]=p2interp[U1]/(ptrgeom->gdet);
688  //if(dopl[U2]) pr[U2]=p2interp[U2]/(ptrgeom->gdet);
689  //if(dopl[U3]) pr[U3]=p2interp[U3]/(ptrgeom->gdet);
690 
691  if(dir==1||dir==3){
692  if(dopl[U1]) pr[U1]=p2interp[U1]/(ptrgeom->gdet);
693  if(dopl[U2]) pr[U2]=p2interp[U2]/(ptrgeom->gdet);
694  if(dopl[U3]) pr[U3]=p2interp[U3]/(ptrgeom->gdet);
695 
696  if(URAD1>=0){
697  if(dopl[URAD1]) pr[URAD1]=p2interp[URAD1]/(ptrgeom->gdet);
698  if(dopl[URAD2]) pr[URAD2]=p2interp[URAD2]/(ptrgeom->gdet);
699  if(dopl[URAD3]) pr[URAD3]=p2interp[URAD3]/(ptrgeom->gdet);
700  }
701 
702  if(dopl[B1]) pr[B1]=p2interp[B1]/(ptrgeom->gdet);
703  if(dopl[B2]) pr[B2]=p2interp[B2]/(ptrgeom->gdet);
704  if(dopl[B3]) pr[B3]=p2interp[B3]/(ptrgeom->gdet);
705  }
706  if(0&&dir==2){
707  //if(dopl[U1]) pr[U1]=p2interp[U1]/(ptrgeom->gdet);
708  //if(dopl[U2]) pr[U2]=p2interp[U2]/(ptrgeom->gdet);
709  //if(dopl[U3]) pr[U3]=p2interp[U3]/(ptrgeom->gdet);
710 
711  if(dopl[B1]) pr[B1]=p2interp[B1]/(ptrgeom->gdet);
712  if(dopl[B2]) pr[B2]=p2interp[B2]/(ptrgeom->gdet);
713  if(dopl[B3]) pr[B3]=p2interp[B3]/(ptrgeom->gdet);
714  }
715  }
716  else{
717  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
718  myexit(100);
719  }
720 
721 
722 #elif(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION_WALD)
723 
724 
725  if(which==DORESCALE){ // rescale before interpolation
726 
727 
728  PINTERPLOOP(pliter,pl) p2interp[pl]=pr[pl];
729 
730  // if(dopl[U1]) p2interp[U1]=pr[U1]*(ptrgeom->gdet);
731  // if(dopl[U2]) p2interp[U2]=pr[U2]*(ptrgeom->gdet);
732  // if(dopl[U3]) p2interp[U3]=pr[U3]*(ptrgeom->gdet);
733 
734  if(dir==1||dir==3){
735  if(dopl[U1]) p2interp[U1]=pr[U1]*(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*(V[1]));
736  if(dopl[U2]) p2interp[U2]=pr[U2]*(sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*(V[1]));
737  if(dopl[U3]) p2interp[U3]=pr[U3]*(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
738 
739  if(URAD1>=0){
740  if(dopl[URAD1]) p2interp[URAD1]=pr[URAD1]*(ptrgeom->gdet);
741  if(dopl[URAD2]) p2interp[URAD2]=pr[URAD2]*(ptrgeom->gdet);
742  if(dopl[URAD3]) p2interp[URAD3]=pr[URAD3]*(ptrgeom->gdet);
743  }
744 
745  if(dopl[B1]) p2interp[B1]=pr[B1]*sqrt(fabs(ptrgeom->gcov[GIND(1,1)]));
746  if(dopl[B2]) p2interp[B2]=pr[B2]*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]));
747  if(dopl[B3]) p2interp[B3]=pr[B3]*(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
748  }
749  if(dir==2){
750  if(dopl[U1]) p2interp[U1]=pr[U1]*(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*(V[1]));
751  if(dopl[U2]) p2interp[U2]=pr[U2]*(sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*(V[1]));
752  if(dopl[U3]) p2interp[U3]=pr[U3]*(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
753 
754  if(dopl[B1]) p2interp[B1]=pr[B1]*sqrt(fabs(ptrgeom->gcov[GIND(1,1)]));
755  if(dopl[B2]) p2interp[B2]=pr[B2]*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]));
756  if(dopl[B3]) p2interp[B3]=pr[B3]*(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
757  }
758  }
759  else if(which==UNRESCALE){ // unrescale after interpolation
760 
761  PINTERPLOOP(pliter,pl) pr[pl]=p2interp[pl];
762 
763  //if(dopl[U1]) pr[U1]=p2interp[U1]/(ptrgeom->gdet);
764  //if(dopl[U2]) pr[U2]=p2interp[U2]/(ptrgeom->gdet);
765  //if(dopl[U3]) pr[U3]=p2interp[U3]/(ptrgeom->gdet);
766 
767  if(dir==1||dir==3){
768  if(dopl[U1]) pr[U1]=p2interp[U1]/(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*(V[1]));
769  if(dopl[U2]) pr[U2]=p2interp[U2]/(sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*(V[1]));
770  if(dopl[U3]) pr[U3]=p2interp[U3]/(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
771 
772  if(URAD1>=0){
773  if(dopl[URAD1]) pr[URAD1]=p2interp[URAD1]/(ptrgeom->gdet);
774  if(dopl[URAD2]) pr[URAD2]=p2interp[URAD2]/(ptrgeom->gdet);
775  if(dopl[URAD3]) pr[URAD3]=p2interp[URAD3]/(ptrgeom->gdet);
776  }
777 
778  if(dopl[B1]) pr[B1]=p2interp[B1]/sqrt(fabs(ptrgeom->gcov[GIND(1,1)]));//(ptrgeom->gdet);///(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3));
779  if(dopl[B2]) pr[B2]=p2interp[B2]/sqrt(fabs(ptrgeom->gcov[GIND(2,2)]));//(ptrgeom->gdet);///(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3));
780  if(dopl[B3]) pr[B3]=p2interp[B3]/(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));//(ptrgeom->gdet);///(pow(V[1],3));
781  }
782  if(dir==2){
783  if(dopl[U1]) pr[U1]=p2interp[U1]/(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*(V[1]));
784  if(dopl[U2]) pr[U2]=p2interp[U2]/(sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*(V[1]));
785  if(dopl[U3]) pr[U3]=p2interp[U3]/(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
786 
787  if(dopl[B1]) pr[B1]=p2interp[B1]/sqrt(fabs(ptrgeom->gcov[GIND(1,1)]));
788  if(dopl[B2]) pr[B2]=p2interp[B2]/sqrt(fabs(ptrgeom->gcov[GIND(2,2)]));
789  if(dopl[B3]) pr[B3]=p2interp[B3]/(sqrt(fabs(ptrgeom->gcov[GIND(3,3)]))*(V[1]));
790  }
791  }
792  else{
793  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
794  myexit(100);
795  }
796 
797 
798 #endif // end over choices for VARTOINTERP
799 
800 
801 
802 
803 
804 
805 
806 
808  //
809  // SEPARATE FIELD-ONLY TRANSFORMATIONS. Assumes VARTOINTERP=PRIMTOINTERP so nothing done otherwise for fields.
810  //
812 
813 
814 #if(VARTOINTERPFIELD==NOFIELDRESCALE)
815 
816  // DO NOTHING, assume VARTOINTERP takes care of any transformation or non-transformation
817 
818 #elif(VARTOINTERPFIELD==PULSARFIELD)
819 
820 
821  if(dir==1){
822  if(which==DORESCALE){ // rescale before interpolation
823 
824  Bconin[0]=0.0;
825  Bconin[1]=pr[B1];
826  Bconin[2]=pr[B2];
827  Bconin[3]=pr[B3];
828 
829  dxdxprim(X, V, dxdxp);
830 
831  getconsts(Bconin, V, ptrgeom, dxdxp,Bconout);
832 
833  if(dopl[B1]) p2interp[B1]=Bconout[1];
834  if(dopl[B2]) p2interp[B2]=Bconout[2];
835  if(dopl[B3]) p2interp[B3]=Bconout[3];
836  }
837  else if(which==UNRESCALE){ // unrescale after interpolation
838 
839  Bconin[0]=0.0;
840  Bconin[1]=p2interp[B1];
841  Bconin[2]=p2interp[B2];
842  Bconin[3]=p2interp[B3];
843 
844  dxdxprim(X, V, dxdxp);
845 
846  undoconsts(Bconin, V, ptrgeom, dxdxp, Bconout);
847 
848  if(dopl[B1]) pr[B1]=Bconout[1];
849  if(dopl[B2]) pr[B2]=Bconout[2];
850  if(dopl[B3]) pr[B3]=Bconout[3];
851 
852  }
853  else{
854  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
855  myexit(100);
856  }
857  }
858  // problem when used with dir=2 since near axis end up dividing by 0
859 
860 
861 #elif(VARTOINTERPFIELD==PULSARFIELD2)
862 
863 
864  if(dir==1){
865  if(which==DORESCALE){ // rescale before interpolation
866 
867  Bconin[0]=0.0;
868  Bconin[1]=pr[B1];
869  Bconin[2]=pr[B2];
870  Bconin[3]=pr[B3];
871 
872  dxdxprim(X, V, dxdxp);
873 
874  if(dopl[B1]) p2interp[B1]=Bconin[1]*dxdxp[1][1]*pow(V[1],3);
875  if(dopl[B2]) p2interp[B2]=Bconin[2]*dxdxp[2][2]*pow(V[1],4);
876  if(dopl[B3]) p2interp[B3]=Bconin[3]*pow(V[1],3);
877  }
878  else if(which==UNRESCALE){ // unrescale after interpolation
879 
880  dxdxprim(X, V, dxdxp);
881 
882  Bconin[0]=0.0;
883  Bconin[1]=p2interp[B1]/(dxdxp[1][1]*pow(V[1],3));
884  Bconin[2]=p2interp[B2]/(dxdxp[2][2]*pow(V[1],4));
885  Bconin[3]=p2interp[B3]/(pow(V[1],3));
886 
887  if(dopl[B1]) pr[B1]=Bconin[1];
888  if(dopl[B2]) pr[B2]=Bconin[2];
889  if(dopl[B3]) pr[B3]=Bconin[3];
890 
891  }
892  else{
893  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
894  myexit(100);
895  }
896  }
897  // problem when used with dir=2 since near axis end up dividing by 0
898 
899 
900 #elif(VARTOINTERPFIELD==PULSARFIELD3)
901 
902 
903  if(dir==1){
904  if(which==DORESCALE){ // rescale before interpolation
905 
906  Bconin[0]=0.0;
907  Bconin[1]=pr[B1];
908  Bconin[2]=pr[B2];
909  Bconin[3]=pr[B3];
910 
911  dxdxprim(X, V, dxdxp);
912 
913  if(dopl[B1]) p2interp[B1]=Bconin[1]*sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3);
914  if(dopl[B2]) p2interp[B2]=Bconin[2]*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*pow(V[1],3);
915  if(dopl[B3]) p2interp[B3]=Bconin[3]*pow(V[1],3);
916  }
917  else if(which==UNRESCALE){ // unrescale after interpolation
918 
919  dxdxprim(X, V, dxdxp);
920 
921  Bconin[0]=0.0;
922  Bconin[1]=p2interp[B1]/(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3));
923  Bconin[2]=p2interp[B2]/(sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3));
924  Bconin[3]=p2interp[B3]/(pow(V[1],3));
925 
926  if(dopl[B1]) pr[B1]=Bconin[1];
927  if(dopl[B2]) pr[B2]=Bconin[2];
928  if(dopl[B3]) pr[B3]=Bconin[3];
929 
930  }
931  else{
932  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
933  myexit(100);
934  }
935  }
936  // problem when used with dir=2 since near axis end up dividing by 0
937 
938 
939 #elif(VARTOINTERPFIELD==GDETVERSION)
940 
941  // dxdxprim_ijk( ptrgeom->i, ptrgeom->j, ptrgeom->k, ptrgeom->p, dxdxp );
942 
943 
944  if(which==DORESCALE){ // rescale before interpolation
945 
946  Bconin[0]=0.0;
947  Bconin[1]=pr[B1];
948  Bconin[2]=pr[B2];
949  Bconin[3]=pr[B3];
950 
951  if(dir==1){
952  if(dopl[B1]) p2interp[B1]=Bconin[1]*(ptrgeom->gdet); //sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3);
953  if(dopl[B2]) p2interp[B2]=Bconin[2]*(ptrgeom->gdet);//*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*pow(V[1],3);
954  if(dopl[B3]) p2interp[B3]=Bconin[3]*(ptrgeom->gdet);//*pow(V[1],3);
955  }
956  else if(dir==2){
957  if(dopl[B1]) p2interp[B1]=Bconin[1];
958  if(dopl[B2]) p2interp[B2]=Bconin[2];
959  if(dopl[B3]) p2interp[B3]=Bconin[3];
960  }
961  else{
962  if(dopl[B1]) p2interp[B1]=Bconin[1];
963  if(dopl[B2]) p2interp[B2]=Bconin[2];
964  if(dopl[B3]) p2interp[B3]=Bconin[3];
965  }
966  }
967  else if(which==UNRESCALE){ // unrescale after interpolation
968 
969  Bconin[0]=0.0;
970  Bconin[1]=p2interp[B1];
971  Bconin[2]=p2interp[B2];
972  Bconin[3]=p2interp[B3];
973 
974  if(dir==1){
975  if(dopl[B1]) pr[B1]=p2interp[B1]/(ptrgeom->gdet);
976  if(dopl[B2]) pr[B2]=p2interp[B2]/(ptrgeom->gdet);
977  if(dopl[B3]) pr[B3]=p2interp[B3]/(ptrgeom->gdet);
978  }
979  else if(dir==2){
980  if(dopl[B1]) pr[B1]=Bconin[1];
981  if(dopl[B2]) pr[B2]=Bconin[2];
982  if(dopl[B3]) pr[B3]=Bconin[3];
983  }
984  else{
985  if(dopl[B1]) pr[B1]=Bconin[1];
986  if(dopl[B2]) pr[B2]=Bconin[2];
987  if(dopl[B3]) pr[B3]=Bconin[3];
988  }
989 
990  }
991  else{
992  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
993  myexit(100);
994  }
995 
996 
997 #elif(VARTOINTERPFIELD==GDETFULLVERSION)
998 
999 
1000  if(which==DORESCALE){ // rescale before interpolation
1001 
1002  if(dir==1||dir==3){
1003  if(dopl[B1]) p2interp[B1]=pr[B1]*(ptrgeom->gdet); //sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3);
1004  if(dopl[B2]) p2interp[B2]=pr[B2]*(ptrgeom->gdet);//*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*pow(V[1],3);
1005  if(dopl[B3]) p2interp[B3]=pr[B3]*(ptrgeom->gdet);//*pow(V[1],3);
1006  }
1007  if(0&&dir==2){
1008  if(dopl[B1]) p2interp[B1]=pr[B1]*(ptrgeom->gdet); //sqrt(fabs(ptrgeom->gcov[GIND(1,1)]))*pow(V[1],3);
1009  if(dopl[B2]) p2interp[B2]=pr[B2]*(ptrgeom->gdet);//*sqrt(fabs(ptrgeom->gcov[GIND(2,2)]))*pow(V[1],3);
1010  if(dopl[B3]) p2interp[B3]=pr[B3]*(ptrgeom->gdet);//*pow(V[1],3);
1011  }
1012  }
1013  else if(which==UNRESCALE){ // unrescale after interpolation
1014 
1015  if(dir==1||dir==3){
1016  if(dopl[B1]) pr[B1]=p2interp[B1]/(ptrgeom->gdet);
1017  if(dopl[B2]) pr[B2]=p2interp[B2]/(ptrgeom->gdet);
1018  if(dopl[B3]) pr[B3]=p2interp[B3]/(ptrgeom->gdet);
1019  }
1020  if(0&&dir==2){
1021  if(dopl[B1]) pr[B1]=p2interp[B1]/(ptrgeom->gdet);
1022  if(dopl[B2]) pr[B2]=p2interp[B2]/(ptrgeom->gdet);
1023  if(dopl[B3]) pr[B3]=p2interp[B3]/(ptrgeom->gdet);
1024  }
1025 
1026  }
1027  else{
1028  dualfprintf(fail_file,"rescale(): no such rescale type! which=%d\n",which);
1029  myexit(100);
1030  }
1031 
1032 
1033 #endif // end over choices for VARTOINTERPFIELD
1034 
1035 
1036 
1037  return(0);
1038 }
1039 
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1