HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
metric.c
Go to the documentation of this file.
1 
17 #include "decs.h"
18 
19 
20 
21 
24 void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc)
25 {
26  void set_gcov_cylminkmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
27  void set_gcov_spcminkmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
28  void set_gcov_cartminkmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
29  void set_gcov_unigravity (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
30  void set_gcov_htmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
31  void set_gcov_htmetric_accurate(FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
32  void set_gcov_ksmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
33  void set_gcov_ks_jp1_metric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
34  void set_gcov_ks_bh_tov_metric (FTYPE *X, FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
35  extern void set_gcov_ks_tov_metric(FTYPE *X, FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
36  extern void set_gcov_bl_tov_metric(FTYPE *X, FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
37  void set_gcov_blmetric (FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc);
38  void gcov2gcovprim(struct of_geom *ptrgeom, FTYPE *X, FTYPE *V, FTYPE *gcovinfunc, FTYPE *gcovpertinfunc, FTYPE *gcovinfuncprim, FTYPE *gcovpertinfuncprim);
39  FTYPE gcovselfpert[NDIM];
40  extern int set_gcov_selfspcmetric(FTYPE *X, FTYPE *V, FTYPE *gcovselfpert);
41  FTYPE V[NDIM],Vmetric[NDIM],Xmetric[NDIM];
42  int j,k;
43  int presenttime;
44  FTYPE phi;
45  LOCALMETRICTEMPVARS;
46 
47 
48 
49 
50 
51 
53  //
54  // determine if asking for metric now or in past
55  //
57  if(DOEVOLVEMETRIC){
58 
59  // if(ptrgeom->i==7 && nstep==1084){
60  // dualfprintf(fail_file,"X[0]=%21.15g Xmetricold[0]=%21.15g Xmetricnew[0]=%21.15g t=%21.15g p=%d\n",X[0],Xmetricold[0],Xmetricnew[0],t,ptrgeom->p);
61  // }
62 
63  if(X[0]>=Xmetricnew[0]-SMALL){ // SMALL trying to account for numerical roundoff error. Assume never want metric at a time difference SMALL before
64  // by present we mean last time latest metric was computed
65  presenttime=1;
66  }
67  else{
68  // should only ever possibly get here if DOEVOLVEMETRIC==1
69  presenttime=0;
70 #if(DOEVOLVEMETRIC!=1)
71  dualfprintf(fail_file,"presenttime=0 and DOEVOLVEMETRIC==0 are incompatible: X[0]=%21.15g t=%21.15g\n",X[0],t);
72  myexit(5523);
73 #endif
74  if(ptrgeom->p==NOWHERE){
75  dualfprintf(fail_file,"Not capable of obtaining past metric that wasn't stored\n");
76  myexit(2356);
77  // for example, can't just interpolate metric since then won't necessarily satisfy divg=0
78  // although probably not big error for anything requiring this interpolation
79  // that is, Connection and most evolution things will use metric at known location
80  // so could infact interpolate if this is needed -- wait and see
81  }
82  if(getprim==0){
83  // not designed to feed back anything not in PRIMECOORDS when asking for past metric
84  dualfprintf(fail_file,"getprim==0 is not compatible with wanting past metric\n");
85  myexit(6246);
86  }
87  }
88  }
89  else presenttime=1; // no evolution, so always present
90 
91 
92 
93 
95  //
96  // Determine current full X, get V, then get g_{\mu\nu}
97  //
99  if(presenttime){
100 
101 
102  // if PRIMECOORDS, did store the data before, and choosing a gridded position, then can just grab from memory
103  if(whichcoord==PRIMECOORDS && didstoremetricdata==1 && ptrgeom->p!=NOWHERE){
104 
105  // get local metric quantities for this loc,i,j,k
106  GETLOCALMETRIC(ptrgeom->p,ptrgeom->i,ptrgeom->j,ptrgeom->k);
107 
108  DLOOP(j,k){
109  gcovinfunc[GIND(j,k)]=localgcov[GIND(j,k)];
110  }
111  DLOOPA(j){
112  gcovpertinfunc[j]=localgcovpert[j];
113  }
114 
115  }
116  else{ // then not defined yet or arbitrary position
117 
118 
119  // before here, X set by i,j,k or chosen directly
120  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
121  //bl_coord(X,V);
122 
123 
124  //dualfprintf(fail_file,"blcoordcalled: i=%d j=%d k=%d p=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p);
125 
126 
128  //
129  // Use V[X] and get Vmetric[V[X]] since all set_gcov's are in terms of Vmetric=r,h,ph (i.e. axisymmetric-type old/original versions)
130  //
132  rotate_VtoVmetric(whichcoord,THETAROT,V,Vmetric);
133  int jj; DLOOPA(jj) Xmetric[jj]=X[jj]; // no change? Depends upon how Xmetric used in set_gcov. GODMARK
134 
135 
136 
138  //
139  // Use Vmetric[] and get metric
140  //
142  if(whichcoord>=0){
143  if(whichcoord==BLCOORDS){
144  set_gcov_blmetric(Vmetric, gcovinfunc, gcovpertinfunc);
145  }
146  else if(whichcoord==KSCOORDS){
147  set_gcov_ksmetric(Vmetric, gcovinfunc, gcovpertinfunc);
148  }
149  else if(whichcoord==KSCARTCOORDS){
150  // set_gcov_kscartmetric(Vmetric, gcovinfunc, gcovpertinfunc); // can just call wrapper around set_gcov_ksmetric()
151  }
152  else if(whichcoord==KS_JP1_COORDS){
153  set_gcov_ks_jp1_metric(Vmetric, gcovinfunc, gcovpertinfunc);
154  }
155  else if(whichcoord==KS_BH_TOV_COORDS){
156  set_gcov_ks_bh_tov_metric(Xmetric, Vmetric, gcovinfunc, gcovpertinfunc);
157  }
158  else if(whichcoord==KS_TOV_COORDS){
159  set_gcov_ks_tov_metric(Xmetric, Vmetric, gcovinfunc, gcovpertinfunc);
160  }
161  else if(whichcoord==BL_TOV_COORDS){
162  set_gcov_bl_tov_metric(Xmetric, Vmetric, gcovinfunc, gcovpertinfunc);
163  }
164  else if(whichcoord==HTMETRIC){
165  set_gcov_htmetric(Vmetric, gcovinfunc, gcovpertinfunc);
166  }
167  else if(whichcoord==HTMETRICACCURATE){
168  set_gcov_htmetric_accurate(Vmetric, gcovinfunc, gcovpertinfunc);
169  }
170  else if(whichcoord==CARTMINKMETRIC || whichcoord==CARTMINKMETRIC2){
171  set_gcov_cartminkmetric(Vmetric, gcovinfunc, gcovpertinfunc);
172  }
173  else if(whichcoord==UNIGRAVITY){
174  set_gcov_unigravity(Vmetric, gcovinfunc, gcovpertinfunc);
175  }
176  else if(whichcoord==CYLMINKMETRIC){
177  set_gcov_cylminkmetric(Vmetric, gcovinfunc, gcovpertinfunc);
178  }
179  else if(whichcoord==SPCMINKMETRIC){
180  set_gcov_spcminkmetric(Vmetric, gcovinfunc, gcovpertinfunc);
181  }
182  else{
183  dualfprintf(fail_file,"gcov_func(): no such whichcoord=%d\n",whichcoord);
184  myexit(12626);
185  }
186 
187  // dualfprintf(fail_file,"blcoordcalled2: i=%d j=%d k=%d p=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p);
188 
189  }
190  else{
191  dualfprintf(fail_file,"your request makes no sense (i.e. can't get prim gcov from prim gcov): getprim=%d whichcoord=%d\n",getprim,whichcoord);
192  myexit(26632);
193  }
194 
196  //
197  // Account for self-gravity term. Add it to g_{\mu\nu}.
198  //
200 #if(DOSELFGRAVVSR)
201  if(
202  whichcoord==BLCOORDS ||
203  whichcoord==KSCOORDS ||
204  whichcoord==KS_JP1_COORDS ||
205  whichcoord==HTMETRIC ||
206  whichcoord==HTMETRICACCURATE ||
207  whichcoord==SPCMINKMETRIC
208  ){
209  // then using spherical polar coordinates and can do self-gravity as designed
210  set_gcov_selfspcmetric(Xmetric,Vmetric,gcovselfpert);
211 
212  // OLD DEBUG
213  // phi = -0.1/Vmetric[1];
214  //gcovinfunc[GIND(TT,TT)] += GINDASSIGNFACTOR(TT,TT)*(-2.0*phi);
215  // gcovpertinfunc[TT] += -2.0*phi;
216  // gcovinfunc[GIND(RR,RR)] += GINDASSIGNFACTOR(RR,RR)*(-2.0*phi);
217  // gcovpertinfunc[RR] += -2.0*phi;
218 
219  // gcovinfunc[GIND(TT,TT)] += GINDASSIGNFACTOR(TT,TT)*gcovselfpert[TT];
220  // gcovpertinfunc[TT] += gcovselfpert[RR];
221  // gcovinfunc[GIND(RR,RR)] += GINDASSIGNFACTOR(RR,RR)*gcovselfpert[TH];
222  // gcovpertinfunc[RR] += gcovselfpert[PH];
223 
224  // if(gcovselfpert[TT]+2.0*phi>0.01){
225  // dualfprintf(fail_file,"got TT difference: %21.15g %21.15g\n",gcovselfpert[TT],-2.0*phi);
226  // }
227 
228  // if(gcovselfpert[RR]+2.0*phi>0.01){
229  // dualfprintf(fail_file,"got RR difference: %21.15g %21.15g\n",gcovselfpert[RR],-2.0*phi);
230  // }
231 
232  // if(gcovselfpert[TH]>0.01){
233  // dualfprintf(fail_file,"got TH difference: %21.15g %21.15g\n",gcovselfpert[TH],0.0);
234  // }
235 
236  // if(gcovselfpert[PH]>0.01){
237  // dualfprintf(fail_file,"got PH difference: %21.15g %21.15g\n",gcovselfpert[PH],0.0);
238  // }
239 
240 
241  // DLOOPA(j){
242  // dualfprintf(fail_file,"t=%21.15g %ld %d :: X1=%21.15g :: postgcov[%d][%d]=%21.15g :: gcovselfpert[%d]=%21.15g mypert=%2.15g\n",t,steppart,nstep,Xmetric[1],j,j,gcovinfunc[GIND(j,j)],j,gcovselfpert[j],-2.0*phi);
243  // }
244 
245  // DLOOPA(j){
246  // dualfprintf(fail_file,"t=%21.15g %ld %d :: X1=%21.15g :: gcovselfpert[%d]=%21.15g\n",t,steppart,nstep,Xmetric[1],j,gcovselfpert[j]);
247  // }
248 
249 
251  //
252  // add self-gravity perturbation to metric
253  //
255  DLOOPA(j){
256  gcovinfunc[GIND(j,j)] += GINDASSIGNFACTOR(j,j)*gcovselfpert[j];
257  gcovpertinfunc[j]+=gcovselfpert[j]; // in this way the perturbation part is always resolved
258  // dualfprintf(fail_file,"gcovselfpert[%d]=%21.15g\n",j,gcovselfpert[j]);
259  }
260  if(whichcoord==KSCOORDS || whichcoord==KS_JP1_COORDS){
261  // KS-form has these terms
262  gcovinfunc[GIND(TT,RR)] += GINDASSIGNFACTOR(TT,RR)*gcovselfpert[TT];
263  gcovinfunc[GIND(RR,TT)] = gcovinfunc[GIND(TT,RR)];
264  }
265  }
266  else if(whichcoord==KS_BH_TOV_COORDS || whichcoord==KS_TOV_COORDS || whichcoord==BL_TOV_COORDS){
267  // then doing full TOV-like solution that doesn't just come in as a perturbation
268 
269  }
270  else{
271  // then not setup for these coordinates
272  dualfprintf(fail_file,"DOSELFGRAVVSR=1 will not work with this whichcoord=%d\n",whichcoord);
273  myexit(145);
274  }
275 #endif
276 
277 
278  // DLOOP(j,k) { stderrfprintf("1gcov[%d][%d]=%21.15g\n",j,k,gcovinfunc[GIND(j,k)]); fflush(stderr);}
279 
280 
281 
283  //
284  // If rotating metric, now we perform transformation on the metric differentials
285  //
286  // X[] inputted is really X, not Xmetric or anything else.
287  // i.e., X[] is associated with V, not Vmetric.
289  transVmetrictoV(whichcoord,ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, THETAROT, X, V, Xmetric, Vmetric, gcovinfunc,gcovpertinfunc);
290 
291 
293  //
294  // Convert to prim coords
295  //
297  if(getprim==1){
298  // all the above are analytic, so have to convert to prim coords.
299  gcov2gcovprim(ptrgeom, Xmetric, Vmetric, gcovinfunc,gcovpertinfunc, gcovinfunc, gcovpertinfunc);
300  }
301  // DLOOP(j,k) { stderrfprintf("2gcov[%d][%d]=%21.15g\n",j,k,gcovinfunc[GIND(j,k)]); fflush(stderr);}
302 
303  // if(ptrgeom->i==7 && nstep==1084){
304  // // DLOOP(j,k) dualfprintf(fail_file,"present time: gcov[%d][%d]=%21.15g\n",j,k,gcovinfunc[GIND(j,k)]);
305  // }
306  }// end if needing to recompute metric
307 
308 
309  }
310  else{
311 
313  //
314  // then not present time, and assume past time metric stored, so recover
315  // also assume never ask for past time at arbitrary location, checked previously
316  //
318 
319  if(ptrgeom->p!=NOWHERE){
320 
321  // get local metric quantities for this loc,i,j,k
322  GETLASTLOCALMETRIC(ptrgeom->p,ptrgeom->i,ptrgeom->j,ptrgeom->k);
323 
324  DLOOP(j,k){
325  gcovinfunc[GIND(j,k)]=localgcov[GIND(j,k)];
326  }
327  DLOOPA(j){
328  gcovpertinfunc[j]=localgcovpert[j];
329  }
330 
331  // if(ptrgeom->i==7 && nstep==1084){
332  // DLOOP(j,k) dualfprintf(fail_file,"past time: gcov[%d][%d]=%21.15g\n",j,k,gcovinfunc[GIND(j,k)]);
333  // }
334  }
335  else{
336  // then don't have grid at right location and need to interpolate
337 #if(NEWMETRICSTORAGE)
338  interpX_gcov(Xmetric, GLOBALPOINT(compgeomlast), NULL, NULL, gcovinfunc, gcovpertinfunc);
339 #else
340  interpX_gcov(Xmetric, NULL, GLOBALPOINT(gcovlast), GLOBALPOINT(gcovpertlast), gcovinfunc, gcovpertinfunc);
341 #endif
342  // above untested, so fail for now
343  dualfprintf(fail_file,"interpX_gcov() unexpectedly called\n");
344  myexit(13515);
345  }
346  }// end else if doing past time
347 
348  //dualfprintf(fail_file,"blcoordcalled3: i=%d j=%d k=%d p=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p);
349 
350 
351 
352 }
353 
354 
356 FTYPE arctanmath(FTYPE x, FTYPE y)
357 {
358  return(atan2(y,x));
359 }
360 
361 
378 int rotate_VtoVmetric(int whichcoord, FTYPE ROTANGLE, FTYPE *V, FTYPE *Vmetric)
379 {
380 
381  if(ROTANGLE!=0.0 && ALLOWMETRICROT && ISSPCMCOORD(whichcoord)){
382  // see metricrot.nb.
383  FTYPE tnew,rnew,hnew,phnew; // true V used in bl_coord() to map to X
384 
385  tnew=V[0];
386  rnew=V[1];
387  hnew=V[2];
388  phnew=V[3];
389 
390 
391  FTYPE told,rold,hold,phold; // what is inside metric functions like set_gcov()
392  told=tnew;
393  rold=rnew;
394 
395  FTYPE b0=ROTANGLE;
396 
397  hold=arctanmath (cos(b0)*cos(hnew) - 1.*cos(phnew)*sin(b0)*sin(hnew),pow(pow(cos(hnew)*sin(b0) + cos(b0)*cos(phnew)*sin(hnew),2) + pow(sin(hnew),2)*pow(sin(phnew),2),0.5));
398 
399  phold=arctanmath (cos(hnew)*sin(b0) + cos(b0)*cos(phnew)*sin(hnew),sin(hnew)*sin(phnew));
400 
401 
402  fix_hp(&hold,&phold);
403 
404  Vmetric[0]=told;
405  Vmetric[1]=rold;
406  Vmetric[2]=hold;
407  Vmetric[3]=phold;
408 
409  }
410  else{
411  int jj; DLOOPA(jj) Vmetric[jj]=V[jj];
412  }
413 
414 
415  return(0);
416 
417 }
418 
419 
421 int fix_hp(FTYPE *h, FTYPE *p)
422 {
423  FTYPE th=*h,ph=*p;
424 
425 #if(0)
426  // keep \theta between 0 and \pi
427  if(th<0.0) th+=M_PI;
428  if(th>=M_PI) th-=M_PI;
429 
430  // keep \phi between 0 and 2\pi
431  if(ph<0.0) ph+=2.0*M_PI;
432  if(ph>=2.0*M_PI) ph-=2.0*M_PI;
433 #else
434 
435  // keep \phi between 0 and 2\pi. Can always do that full rotation.
436  // assume never more out of phase that 1 full rotation
437  if(ph<0.0) ph+=2.0*M_PI;
438  if(ph>=2.0*M_PI) ph-=2.0*M_PI;
439 
440  // keep \theta between 0 and \pi and \phi between 0 and 2\pi
441  // but need to be at same physical SPC location, not arbitrary rotation
442  if(th>=0.0 && th<=M_PI){
443  // do nothing
444  }
445  else if(th<0.0){
446  th*=-1.0;
447  if(ph<=M_PI) ph+=M_PI;
448  else if(ph>M_PI) ph-=M_PI;
449  else{ dualfprintf(fail_file,"Shouldn't be here1 with th=%g ph=%g\n",th,ph); myexit(1);}
450  }
451  else if(th>=M_PI){
452  th=M_PI-th;
453  if(ph<=M_PI) ph+=M_PI;
454  else if(ph>M_PI) ph-=M_PI;
455  else{ dualfprintf(fail_file,"Shouldn't be here2 with th=%g ph=%g\n",th,ph); myexit(1);}
456  }
457  else{ dualfprintf(fail_file,"Shouldn't be here3 with th=%g ph=%g\n",th,ph); myexit(1);}
458 #endif
459 
460  *h=th;
461  *p=ph;
462 
463  return(0);
464 
465 }
466 
485 int rotate_VmetrictoV(int whichcoord, FTYPE ROTANGLE, FTYPE *Vmetric, FTYPE *V)
486 {
487 
488  if(ROTANGLE!=0.0 && ALLOWMETRICROT && ISSPCMCOORD(whichcoord)){
489  // see metricrot.nb. Note that mathematica's Arctan[x,y] = C's atan2(y,x) (i.e. args are flipped in order)
490  FTYPE t,r,h,ph;
491 
492  t=Vmetric[0];
493  r=Vmetric[1];
494  h=Vmetric[2];
495  ph=Vmetric[3];
496 
497  FTYPE xc,yc,zc,Rold;
498  xc=mysin(h)*mycos(ph);
499  yc=mysin(h)*mysin(ph);
500  zc=mycos(h);
501  Rold=sqrt(xc*xc + yc*yc);
502 
503  // rotation around y-axis using right-hand rule
504  FTYPE xcnew,ycnew,zcnew,Rnew;
505  xcnew=xc*mycos(ROTANGLE)-zc*mysin(ROTANGLE);
506  ycnew=yc;
507  zcnew=zc*mycos(ROTANGLE)+xc*mysin(ROTANGLE);
508  Rnew=sqrt(xcnew*xcnew + ycnew*ycnew);
509 
510  // Below uses atan2 so gets back correct quadrant
511  FTYPE tnew,rnew,hnew,phnew;
512  tnew=t;
513  rnew=r;
514  // these return in range [-\pi,\pi]
515  hnew=atan2(Rnew,zcnew);
516  phnew=atan2(ycnew,xcnew);
517 
518 
519  fix_hp(&hnew,&phnew);
520 
521 
522  V[0]=tnew;
523  V[1]=rnew;
524  V[2]=hnew;
525  V[3]=phnew;
526 
527  }
528  else{
529  int jj; DLOOPA(jj) V[jj]=Vmetric[jj];
530  }
531 
532 
533 
534  return(0);
535 }
536 
537 
538 
541 int interpX_gcov(FTYPE *X, struct of_compgeom PTRDEFMETMACP1A0(compgeom,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3), FTYPE PTRDEFMETMACP1A2(gcovgrid,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM,NDIM), FTYPE PTRDEFMETMACP1A1(gcovpertgrid,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM), FTYPE *gcov, FTYPE *gcovpert)
542 {
543  int i,j,k;
544  int jj,kk;
545  extern void icoord(FTYPE *X,int loc, int *i, int *j, int *k);
546  int ip,jp,kp;
547  FTYPE Xijk[NDIM],Xip[NDIM],Xjp[NDIM],Xkp[NDIM];
548  FTYPE gijk,gip,gjp,gkp;
549  FTYPE dist[4+1],totaldist;
550  int loc;
551 
552 
553  // CENT is chosen as reference location for all interpolations
554  loc=CENT;
555 
556  // use icoord() to get i,j,k from X so can get potential at any location from interpolation of existing potential at certain locations
557 
558  // get centered i,j,k
559  icoord(X,loc,&i,&j,&k); // truncates to lower i
560  coord_ijk(i,j,k, loc, Xijk);
561  // only want effective grid location (i.e. normalized X position) so that differences are 1 between grid cells -- don't care about absolute magnitude of X once normalized
562  DLOOPA(jj) Xijk[jj]/=dx[jj];
563 
565  // X1 dir
566  //
567  // limit lower value of i
568  if(i<=-N1BND-SHIFT1){
569  i=-N1BND;
570  }
571  // limit upper value of i
572  if(i>=N1+N1BND-SHIFT1){ // -1 is to offset so +1 is in range
573  i=N1+N1BND-SHIFT1-SHIFT1; // first -1 is to offset from out of range, and second -1 is so ip is located within array range
574  }
575 
576  // get X-position of coordinates i and i+1
577  ip=i+SHIFT1;
578  coord_ijk(ip,j,k, loc, Xip);
579  DLOOPA(jj) Xip[jj]/=dx[jj];
580 
581 
582 
584  // X2 dir
585  //
586  // limit lower value of j
587  if(j<=-N2BND-SHIFT2){
588  j=-N2BND;
589  }
590  // limit upper value of j
591  if(j>=N2+N2BND-SHIFT2){ // -1 is to offset so +1 is in range
592  j=N2+N2BND-SHIFT2-SHIFT2; // first -1 is to offset from out of range, and second -1 is so jp is located within array range
593  }
594 
595  // get X-position of coordinates j and j+1
596  jp=j+SHIFT2;
597  coord_ijk(i,jp,k, loc, Xjp);
598  DLOOPA(jj) Xjp[jj]/=dx[jj];
599 
600 
601 
603  // X3 dir
604  //
605  // limit lower value of k
606  if(k<=-N3BND-SHIFT3){
607  k=-N3BND;
608  }
609  // limit upper value of k
610  if(k>=N3+N3BND-SHIFT3){ // -1 is to offset so +1 is in range
611  k=N3+N3BND-SHIFT3-SHIFT3; // first -1 is to offset from out of range, and second -1 is so kp is located within array range
612  }
613 
614  // get X-position of coordinates k and k+1
615  kp=k+SHIFT3;
616  coord_ijk(i,j,kp, loc, Xkp);
617  DLOOPA(jj) Xkp[jj]/=dx[jj];
618 
619 
620  // now use some interpolation of 3 points
621 
622  // use this position to interpolate in X (which happens to be the same as interpolating in ijk)
623  // distance away from one point toward ijk
624  dist[1]=(1-(X[1]-Xijk[1]))*(1-(X[2]-Xijk[2]))*(1-(X[3]-Xijk[3]));
625  // away from one point toward ip
626  dist[2]=(1-(X[1]-Xip[1]))*(1-(X[2]-Xip[2]))*(1-(X[3]-Xip[3]))*(1.0-dist[1]);
627  // away from one point toward jp
628  dist[3]=(1-(X[1]-Xjp[1]))*(1-(X[2]-Xjp[2]))*(1-(X[3]-Xjp[3]))*(1.0-dist[2])*(1.0-dist[1]);
629  // away from one point toward kp
630  dist[4]=(1-(X[1]-Xkp[1]))*(1-(X[2]-Xkp[2]))*(1-(X[3]-Xkp[3]))*(1.0-dist[3])*(1.0-dist[2])*(1.0-dist[1]);
631  // normalization of total distance
632  totaldist=dist[1]+dist[2]+dist[3]+dist[4];
633 
634 
635 
636  DLOOP(jj,kk){
637 #if(NEWMETRICSTORAGE)
638  gijk=METMACP1A0(compgeom,loc,i,j,k).gcov[GIND(jj,kk)];
639  gip=METMACP1A0(compgeom,loc,ip,j,k).gcov[GIND(jj,kk)];
640  gjp=METMACP1A0(compgeom,loc,i,jp,k).gcov[GIND(jj,kk)];
641  gkp=METMACP1A0(compgeom,loc,i,j,kp).gcov[GIND(jj,kk)];
642 #else
643  gijk=METMACP1A2(gcovgrid,loc,i,j,k,jj,kk);
644  gip=METMACP1A2(gcovgrid,loc,ip,j,k,jj,kk);
645  gjp=METMACP1A2(gcovgrid,loc,i,jp,k,jj,kk);
646  gkp=METMACP1A2(gcovgrid,loc,i,j,kp,jj,kk);
647 #endif
648 
649  // X's are per unit dX's
650  gcov[GIND(jj,kk)] = (gijk*dist[1]+gip*dist[2]+gjp*dist[3]+gkp*dist[4])/totaldist;
651 
652  }
653 
654  DLOOPA(jj){
655 #if(NEWMETRICSTORAGE)
656  gijk=METMACP1A0(compgeom,loc,i,j,k).gcovpert[jj];
657  gip=METMACP1A0(compgeom,loc,ip,j,k).gcovpert[jj];
658  gjp=METMACP1A0(compgeom,loc,i,jp,k).gcovpert[jj];
659  gkp=METMACP1A0(compgeom,loc,i,j,kp).gcovpert[jj];
660 #else
661  gijk=METMACP1A1(gcovpertgrid,loc,i,j,k,jj);
662  gip=METMACP1A1(gcovpertgrid,loc,ip,j,k,jj);
663  gjp=METMACP1A1(gcovpertgrid,loc,i,jp,k,jj);
664  gkp=METMACP1A1(gcovpertgrid,loc,i,j,kp,jj);
665 #endif
666  // X's are per unit dX's
667  gcovpert[jj] = (gijk*dist[1]+gip*dist[2]+gjp*dist[3]+gkp*dist[4])/totaldist;
668 
669  }
670 
671  return(0);
672 
673 }
674 
675 
676 
678 void gcon_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon)
679 {
680  void set_gcon_blmetric(FTYPE *V, FTYPE *gcon);
681  void set_gcon_ksmetric(FTYPE *V, FTYPE *gcon);
682  void gcon2gconprim(struct of_geom *ptrgeom, FTYPE *X, FTYPE *V, FTYPE *gcon,FTYPE *gconprim);
683  void matrix_inverse_metric(int whichcoord, FTYPE *gcov, FTYPE *gcon);
684  int j,k;
685  FTYPE V[NDIM];
686 
687 
688  if(whichcoord>=0){
689  if(whichcoord==BLCOORDS){
690  //dualfprintf(fail_file,"mi in BLCOORDS\n");
691  if(ANALYTICGCON){
692  bl_coord(X, V);
693  set_gcon_blmetric(V, gcon);
694  if(getprim) gcon2gconprim(ptrgeom, X, V, gcon,gcon);
695  }
696  // since don't have gcon and want to keep things simple by only having to specify gcov
697  else matrix_inverse_metric(whichcoord,gcov,gcon);
698 
699  }
700  else if(whichcoord==KSCOORDS){
701  //dualfprintf(fail_file,"mi in KSCOORDS\n");
702  //DLOOP(j,k) dualfprintf(fail_file,"ks gcov[%d][%d]=%g\n",j,k,gcov[GIND(j,k)]);
703 
704  // DEBUG:
705  //bl_coord(X, V);
706  //dualfprintf(fail_file,"i=%d j=%d k=%d loc=%d :: %21.15g\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p,V[1]);
707 
708 
709  if(ANALYTICGCON){
710  bl_coord(X, V);
711  set_gcon_ksmetric(V,gcon);
712  if(getprim) gcon2gconprim(ptrgeom, X, V,gcon,gcon);
713  }
714  // since don't have gcon and want to keep things simple by only having to specify gcov
715  else matrix_inverse_metric(whichcoord,gcov,gcon);
716  }
717  else if(whichcoord==KS_JP1_COORDS || whichcoord==HTMETRIC || whichcoord==HTMETRICACCURATE || whichcoord==CARTMINKMETRIC || whichcoord==CARTMINKMETRIC2 || whichcoord==UNIGRAVITY || whichcoord==CYLMINKMETRIC || whichcoord==SPCMINKMETRIC || whichcoord==KS_BH_TOV_COORDS || whichcoord==KS_TOV_COORDS || whichcoord==BL_TOV_COORDS || whichcoord==KSCARTCOORDS){
718  // do not have analytic gcon, so invert numerically
719  matrix_inverse_metric(whichcoord,gcov,gcon);
720  }
721  else{
722  dualfprintf(fail_file,"gcon_func(): no such whichcoord=%d\n",whichcoord);
723  myexit(2687);
724  }
725  }
726  else{
727  dualfprintf(fail_file,"your request makes no sense (i.e. can't get prim gcon from prim gcon\n");
728  myexit(3783);
729  }
730 
731 
732 
733 }
734 
735 
736 
737 
738 
739 
740 
741 
742 
744 void conn_func(int whichcoord, FTYPE *X, struct of_geom *geom,
745  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
746 {
747  void set_conn_general(FTYPE *X, struct of_geom *geom,
748  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
749  void set_conn_ksmetric(FTYPE *X, struct of_geom *geom,
750  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
751  void set_conn_cylminkmetric(FTYPE *X, struct of_geom *geom,
752  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
753  void set_conn_cartminkmetric(FTYPE *X, struct of_geom *geom,
754  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
755 
756 
757  if(whichcoord==KS_JP1_COORDS || whichcoord==BLCOORDS || whichcoord==KS_BH_TOV_COORDS || whichcoord==KS_TOV_COORDS || whichcoord==BL_TOV_COORDS || whichcoord==HTMETRIC || whichcoord==HTMETRICACCURATE || whichcoord==UNIGRAVITY || whichcoord==SPCMINKMETRIC || whichcoord==KSCARTCOORDS){
758  set_conn_general(X,geom,conn,conn2);
759  }
760  else if(whichcoord==KSCOORDS){
761  set_conn_ksmetric(X,geom,conn,conn2);
762  }
763  else if(whichcoord==CYLMINKMETRIC){
764  set_conn_cylminkmetric(X,geom,conn,conn2);
765  }
766  else if(whichcoord==CARTMINKMETRIC || whichcoord==CARTMINKMETRIC2){
767  set_conn_cartminkmetric(X,geom,conn,conn2);
768  }
769  else{
770  dualfprintf(fail_file,"conn_func(): no such whichcoord=%d\n",whichcoord);
771  myexit(7367);
772  }
773 }
774 
775 
776 
781 void eomfunc_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *EOMFUNCNAME)
782 {
783 
784  FTYPE gcovmcoord[SYMMATRIXNDIM];
785  FTYPE gcovpertcoord[NDIM];
786  FTYPE r,th;
787  int j,k;
788  int pl,pliter;
789  FTYPE V[NDIM];
790  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
791  FTYPE ftemp;
792 
793 
794  if(WHICHEOM!=WITHGDET || NEWMETRICSTORAGE==1){
795  // then do eomfunc_func() calculation
796  }
797  else{
798  // don't
799  // ensures avoid memory leak even if users mistakenly calls this function
800  dualfprintf(fail_file,"Don't call eomfunc_func() whichcoord=%d\n",whichcoord);
801  return;
802  }
803 
804 
805 
806 
807  if(WHICHEOM==WITHGDET){
808  gcov_func(ptrgeom, getprim, whichcoord,X,gcovmcoord,gcovpertcoord); // actually returns primcoords version of whichcoord
809  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
810  if(gdet_func_metric(whichcoord,V,gcovmcoord,&ftemp)!=0){
811  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() problem in eomfunc_func()\n");
812  }
813  PLOOP(pliter,pl) EOMFUNCASSIGN(pl)=ftemp;
814  }
815  else if(WHICHEOM==WITHNOGDET){
816  PLOOP(pliter,pl) EOMFUNCASSIGN(pl)=1.0;
817  }
818  else if(WHICHEOM==WITHSINSQ){ // obviously coordinate dependent
819  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
820 
821  // rotate consistently with rotation in metric.c
822  // that is, anytime take i,j,k -> V for metric stuff, need to first rotate V.
823  r=V[1]; th=V[2]; // NOTEMARK :No rotation_V here because th and sin(th)^2 below are for polar grid itself, not the metric.
824  ftemp=sin(th)*sin(th);
825  PLOOP(pliter,pl) EOMFUNCASSIGN(pl)=ftemp;
826  }
827  else{
828  dualfprintf(fail_file,"eomfunc_func(): no such WHICHEOM=%d\n",WHICHEOM);
829  myexit(798436);
830  }
831 
832  if(FORCEGDETPOSITIVE==1){
833  PLOOP(pliter,pl) EOMFUNCASSIGN(pl)=fabs(EOMFUNCASSIGN(pl));
834  }
835 
836 }
837 
838 
839 
854 void set_gcov_htmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
855 {
856  FTYPE P2,Q12, Q22;
857  FTYPE z,OO,RO,SS,AA,alphasq,betaphi;
858  FTYPE SSpert,OOpert,AApert;
859  FTYPE M,J,Q;
860  FTYPE r,th;
861  FTYPE ftemp;
862  // size of r relative to M is controlled by M, not r.
863  FTYPE r2small;
864 #if(SMOOTHSING)
865  FTYPE signr,rsmooth;
866 #endif
867  FTYPE rsharp;
868 
869  r = rsharp = V[1];
870 #if(SMOOTHSING)
871  signr=mysign(r);
872  rsmooth = signr*(fabs(r)+SMALL+drsing);
873  r = rsmooth;
874  r2small = r*r;
875 #else
876  r2small = r*r + SMALL;
877 #endif
878 
879 
880  th=V[2];
881 
882 
883 #define FLATSPACE 0
884 #define NOSPIN 1
885 #define FULLHT 2
886 
887 #define HTMETRICTYPE FULLHT
888 
889  // here G=c=1 and M is either 0 or 1 usually, such that
890  // J=a M^2 , Q~J^2/M =a^2 M
891  // NOW: uses MBH as mass of NS (or any compact object within the boundary)
892  M=MBH;
893 
894  if(HTMETRICTYPE==FULLHT){
895  J=a*MBH; // I\Omega = J
896  Q=a*a*MBH; // approximately, perhaps, like Kerr geometry, if hard enough EOS
897  }
898  else if(HTMETRICTYPE==NOSPIN){
899  a = 0.0;
900  J=a*MBH; // no spin of metric
901  Q=a*a*MBH;
902  }
903  else if(HTMETRICTYPE==FLATSPACE){
904  M=MBH=1E-100; // no mass of star, just flat space (GJ model)
905  a=0.0;
906  J=a*MBH;
907  Q=a*a*MBH;
908  }
909 
910 
911  // surface radius is not constant radius if Q!=0
912  // if Q=0, then surface radius is 3M? depends on EOS.
913 
914  // let's set M=1, such that J is up to M^2 and Q is up to M^3 such that radius is in units of GM/c^2
915 
916  // For Sun: M/r<M/R~2E-6, J/r^2<J/R^2~1E-12 (J/M^2=4E-24), Q/r^3<Q/R^3<~1E-10 (Q/R^3=8E-28) and surface roughly spherical
917 
918  // Fastest pulsar PSR 1937+21 \tau=1.56ms , 640times/sec
919  // Deepto Chakrabarty of MIT
920  // Vela: PSR 0833-45 \tau=89.3ms
921  // strongest: PSR 0329+54 \tau=0.715s
922  // peak: 760times/sec (limit by GR?), 2-3X is theoretical limit but would mass shed.
923  // formation spin rate: 30/sec
924 
925 
926  // see Hartle & Thorne 1968 or Kim Lee Lee Lee 2005
927 
928  // Legendre polynomial
929  P2=0.5*(3.0*cos(th)*cos(th)-1.0);
930 
931  z=rsharp/M-1.0;
932 
933  // Mathematica's LegendreQ (associated Legendre of 2nd kind) is such that:
934  // Q12=Im[Legendre[2,1,z]] and Q22=-Re[Legendre[2,2,z]]
935  Q12=sqrt(z*z-1.0)*( (3.0*z*z-2.0)/(z*z-1.0)-3.0/2.0*z*log((z+1)/(z-1)) );
936 
937  Q22=(3.0/2.0*(z*z-1.0)*log((z+1)/(z-1))-(3.0*z*z*z-5.0*z)/(z*z-1.0) );
938 
939 
940 
941  // metric stuff
942  OOpert=-2.0*M/r+2.0*J*J/(r*r*r*r);
943  OO=1.0+OOpert;
944 
945  RO=1.0+2.0*(J*J/(M*r*r*r)*(1.0+M/r)+5.0/8.0*(Q-J*J/M)/(M*M*M)*Q22)*P2;
946 
947  SSpert=-2.0*(J*J/(M*r*r*r)*(1.0-5.0*M/r)+5.0/8.0*(Q-J*J/M)/(M*M*M)*Q22)*P2;
948  SS=1.0+SSpert;
949 
950  AApert=+2.0*(-J*J/(M*r*r*r)*(1.0+2.0*M/r)+5.0/8.0*(Q-J*J/M)/(M*M*M)*( (2.0*M/(sqrt(r*r*(1.0-2.0*M/r))))*Q12 - Q22 ))*P2;
951  AA=1.0+AApert;
952 
953  alphasq=OO*RO;
954 
955  betaphi=-2.0*J/(r*r*r); // note that omega_{zamo}=-betaphi
956 
957  gcov[GIND(PH,PH)] = rsharp*rsharp*AA*sin(th)*sin(th) ;
958  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0;
959 
960  gcov[GIND(TT,TT)] = -(alphasq-betaphi*betaphi*gcov[GIND(PH,PH)]) ;
961  ftemp = (1.0-alphasq);
962  // g_{tt} + 1
963  gcovpert[TT] = ftemp + betaphi*betaphi*gcov[GIND(PH,PH)] ;
964 
965  gcov[GIND(TT,RR)] = 0.0 ;
966  gcov[GIND(TT,TH)] = 0.0 ;
967  gcov[GIND(TT,PH)] = betaphi*gcov[GIND(PH,PH)] ;
968 
969  gcov[GIND(RR,TT)] = 0.0 ;
970  gcov[GIND(RR,RR)] = SS/OO ;
971  // SS/OO = (SSpert+1.0)/OO = SSpert/OO + 1.0/OO , so SS/OO-1.0 = SSpert/OO + (1.0/OO-1.0)
972  ftemp=(1.0/OO - 1.0);
973  gcovpert[RR] = SSpert/OO + ftemp ; // might help with radial gravity terms
974 
975  gcov[GIND(RR,TH)] = 0.0 ;
976  gcov[GIND(RR,PH)] = 0.0 ;
977 
978  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
979  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
980  gcov[GIND(TH,TH)] = rsharp*rsharp*AA ;
981  // r^2 AA - 1 = r^2(AApert+1.0) -1 = r^2 AApert + r^2 -1, so not a big issue with catastrophic cancellation in non-rel limit
982  // probably have to define background as being r^2 to have gravity work at large radii (GODMARK)
983  gcovpert[TH] = gcov[GIND(TH,TH)] - 1.0;
984 
985  gcov[GIND(TH,PH)] = 0.0 ;
986 
987  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
988  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
989  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
990 
991 
992 
993 
994 
995 
996 }
997 
998 #undef HTMETRICTYPE
999 
1001 void set_gcov_htmetric_accurate(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1002 {
1003  FTYPE M,J,Q;
1004  FTYPE u,p,A1,A2,W,F1,F2,H1,H2,L,G1;
1005  FTYPE r,th;
1006  FTYPE ftemp;
1007  FTYPE r2small;
1008 #if(SMOOTHSING)
1009  FTYPE signr,rsmooth;
1010 #endif
1011  FTYPE rsharp;
1012 
1013  r = rsharp = V[1];
1014 #if(SMOOTHSING)
1015  signr=mysign(r);
1016  rsmooth = signr*(fabs(r)+SMALL+drsing);
1017  r = rsmooth;
1018  r2small = r*r;
1019 #else
1020  r2small = r*r + SMALL;
1021 #endif
1022 
1023 
1024  th=V[2];
1025 
1026 
1027 
1028 #define FLATSPACE 0
1029 #define NOSPIN 1
1030 #define FULLHT 2
1031 
1032 #define HTMETRICTYPE FULLHT
1033 
1034  // here G=c=1 and M is either 0 or 1 usually, such that
1035  // J=a M^2 , Q~J^2/M =a^2 M
1036  // NOW: uses MBH as mass of NS (or any compact object within the boundary)
1037  M=MBH;
1038 
1039  if(HTMETRICTYPE==FULLHT){
1040  J=a*MBH; // I\Omega = J
1041  Q=a*a*MBH; // approximately, perhaps, like Kerr geometry, if hard enough EOS
1042  }
1043  else if(HTMETRICTYPE==NOSPIN){
1044  a = 0.0;
1045  J=a*MBH; // no spin of metric
1046  Q=a*a*MBH;
1047  }
1048  else if(HTMETRICTYPE==FLATSPACE){
1049  M=MBH=1E-100; // no mass of star, just flat space (GJ model)
1050  a=0.0;
1051  J=a*MBH;
1052  Q=a*a*MBH;
1053  }
1054 
1055 
1056 
1057 
1058 
1059  L=(80.0*pow(M,6)+8.0*pow(M,4)*r*r+10.0*M*M*M*r*r*r+20.0*M*M*pow(r,4)-45.0*M*pow(r,5)+15.0*pow(r,6));
1060 
1061  u=cos(th);
1062  p=1.0/(8.0*M*pow(r,4)*(r-2.0*M)); // typo in paper for parenthesis
1063 
1064 
1065  A1=(15.0*r*(r-2.0*M)*(1.0-3.0*u*u))/(16.0*M*M)*log(r/(r-2.0*M));
1066  A2=(15.0*(r*r-2.0*M*M)*(3.0*u*u-1.0))/(16.0*M*M)*log(r/(r-2.0*M));
1067 
1068  W=(r-M)*(16.0*pow(M,5)+8.0*pow(M,4)*r-10*M*M*r*r*r-30.0*M*pow(r,4)+15.0*pow(r,5))+u*u*(48.0*pow(M,6)-8.0*pow(M,5)*r-24.0*pow(M,4)*r*r-30.0*M*M*M*r*r*r-60.0*M*M*pow(r,4)+135.0*M*pow(r,5)-45.0*pow(r,6));
1069 
1070  F1=-p*W+A1;
1071  F2=5.0*r*r*r*p*(3.0*u*u-1.0)*(r-M)*(2.0*M*M+6.0*M*r-3.0*r*r)-A1;
1072 
1073  H1=A2+(1.0/(8.0*M*r*r*r*r))*(1.0-3.0*u*u) * (16.0*pow(M,5)+8.0*pow(M,4)*r-10.0*M*M*r*r*r+15.0*M*pow(r,4)+15.0*pow(r,5)) ;
1074  H2=-A2 + (1.0/(8.0*M*r))*5.0*(1.0-3.0*u*u)*(2.0*M*M-3.0*M*r-3.0*r*r);
1075 
1076  G1=p*((L-72.0*M*M*M*M*M*r)-3.0*u*u*(L-56.0*M*M*M*M*M*r))-A1;
1077 
1078 
1079 
1080 
1081 
1082  ftemp = J*J*F1-Q*F2;
1083  gcov[GIND(TT,TT)] = -(1.0-2.0*M/r)*(1.0+ftemp);
1084  // g_{tt} + 1
1085  // gcovpert[TT] = -(J*J*F1-Q*F2) + (2.0*M/r)*(1.0+J*J*F1-Q*F2);
1086  gcovpert[TT] = -ftemp +(1.0+ftemp)*(2.0*M/r);
1087 
1088 
1089 
1090  gcov[GIND(TT,RR)] = 0.0 ;
1091  gcov[GIND(TT,TH)] = 0.0 ;
1092  gcov[GIND(TT,PH)] = (2.0*J*M*M/r)*sin(th)*sin(th);
1093 
1094  gcov[GIND(RR,TT)] = 0.0 ;
1095 
1096  ftemp=J*J*G1+Q*F2;
1097  gcov[GIND(RR,RR)] = (1.0/(1.0-2.0*M/r))*(1.0+ftemp);
1098  // g_{rr} - 1
1099  gcovpert[RR] = (ftemp+2.0*M/r)/(1.0-2.0*M/r); // should help with radial gravity in non-rel regime
1100 
1101  gcov[GIND(RR,TH)] = 0.0 ;
1102  gcov[GIND(RR,PH)] = 0.0 ;
1103 
1104  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
1105  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
1106  gcov[GIND(TH,TH)] = rsharp*rsharp*(1.0+J*J*H1-Q*H2) ;
1107  gcovpert[TH] = gcov[GIND(TH,TH)]-1.0;
1108 
1109  gcov[GIND(TH,PH)] = 0.0 ;
1110 
1111  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
1112  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
1113  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
1114  gcov[GIND(PH,PH)] = gcov[GIND(TH,TH)]*sin(th)*sin(th);
1115  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0;
1116 
1117 
1118 
1119 
1120 
1121 
1122 }
1123 #undef HTMETRICTYPE
1124 
1125 
1126 
1129 void set_gcov_cartminkmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1130 {
1131 
1132  gcov[GIND(TT,TT)] = -1.0;
1133  gcovpert[TT] = 0.0;
1134 
1135  gcov[GIND(TT,RR)] = 0.0 ;
1136  gcov[GIND(TT,TH)] = 0.0 ;
1137  gcov[GIND(TT,PH)] = 0.0 ;
1138 
1139  gcov[GIND(RR,TT)] = 0.0 ;
1140  gcov[GIND(RR,RR)] = 1.0 ;
1141  gcovpert[RR] = 0.0;
1142 
1143  gcov[GIND(RR,TH)] = 0.0 ;
1144  gcov[GIND(RR,PH)] = 0.0 ;
1145 
1146  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
1147  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
1148  gcov[GIND(TH,TH)] = 1.0 ;
1149  gcovpert[TH] = 0.0;
1150 
1151  gcov[GIND(TH,PH)] = 0.0 ;
1152 
1153  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
1154  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
1155  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
1156  gcov[GIND(PH,PH)] = 1.0 ;
1157  gcovpert[PH] = 0.0;
1158 
1159 }
1160 
1161 
1163 void set_gcov_unigravity(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1164 {
1165 
1166  FTYPE phi;
1167  FTYPE x,y,z;
1168  FTYPE FORCEMAG;
1169 
1170  x=V[1];
1171  y=V[2];
1172  z=V[3];
1173 
1174  //this is the hard-cored value for test 28, right? SASMARK
1175  FORCEMAG=0.1/(coordparams.timescalefactor*coordparams.timescalefactor); // strength of force, scales like v^2
1176 
1177  // uniform force field in +y direction (i.e. dphi/dy = FORCEMAG, so Force = -FORCEMAG)
1178  phi = FORCEMAG*y;
1179 
1180  gcov[GIND(TT,TT)] = -1.0-2.0*phi;
1181  gcovpert[TT] = -2.0*phi;
1182 
1183  gcov[GIND(TT,RR)] = 0.0 ;
1184  gcov[GIND(TT,TH)] = 0.0 ;
1185  gcov[GIND(TT,PH)] = 0.0 ;
1186 
1187  gcov[GIND(RR,TT)] = 0.0 ;
1188  gcov[GIND(RR,RR)] = 1.0-2.0*phi ;
1189  gcovpert[RR] = -2.0*phi;
1190 
1191  gcov[GIND(RR,TH)] = 0.0 ;
1192  gcov[GIND(RR,PH)] = 0.0 ;
1193 
1194  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
1195  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
1196  gcov[GIND(TH,TH)] = 1.0-2.0*phi ;
1197  gcovpert[TH] = -2.0*phi;
1198 
1199  gcov[GIND(TH,PH)] = 0.0 ;
1200 
1201  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
1202  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
1203  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
1204 
1205  gcov[GIND(PH,PH)] = 1.0-2.0*phi ;
1206  gcovpert[PH] = -2.0*phi;
1207 
1208  // dualfprintf(fail_file,"tsf=%21.15g :: %21.15g %21.15g :: %21.15g %21.15g %21.15g %21.15g\n",coordparams.timescalefactor,FORCEMAG,phi,gcovpert[TT],gcovpert[RR],gcovpert[TH],gcovpert[PH]);
1209 
1210 
1211 
1212 }
1213 
1214 
1215 
1217 void set_gcov_cylminkmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1218 {
1219 
1220  FTYPE r,Mass,RSTAR;
1221  FTYPE phi;
1222  FTYPE R,z;
1223 
1224 
1225 
1226  R = V[1];
1227  z = V[2];
1228 
1229 
1230 
1231 
1232  Mass=MBH; // NOW: use MBH as any compact object mass (in length units)
1233  RSTAR=0.1; // could use dxdxp[RR][RR]*dx[RR] to get size
1234  r=sqrt(R*R+z*z);
1235  // These SMALL's assume not further squaring or otherwise making smaller to yield numerical 0
1236  // These SMALL's also assume eventually apply gdet as 0 at axis.
1237  FTYPE rsq=SMALL+r*r;
1238  FTYPE Rsq=SMALL+R*R;
1239  FTYPE zsq=SMALL+z*z;
1240 
1241  if(r<RSTAR){
1242  phi = (Mass/(2.0*RSTAR))*((r/RSTAR)*(r/RSTAR)-3.0);
1243  }
1244  else phi = -Mass/r;
1245 
1246 
1247  gcov[GIND(TT,TT)] = -1.0-2.0*phi;
1248  gcovpert[TT] = -2.0*phi;
1249 
1250  gcov[GIND(TT,RR)] = 0.0 ;
1251  gcov[GIND(TT,TH)] = 0.0 ;
1252  gcov[GIND(TT,PH)] = 0.0 ;
1253 
1254  gcov[GIND(RR,TT)] = 0.0 ;
1255  gcov[GIND(RR,RR)] = 1.0-2.0*phi*Rsq/(rsq);
1256  gcovpert[RR] = -2.0*phi*Rsq/(rsq);
1257 
1258  gcov[GIND(RR,TH)] = -2.0*phi*R*z/(rsq) ; // order v^4
1259  gcov[GIND(RR,PH)] = 0.0 ;
1260 
1261  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
1262  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
1263  gcov[GIND(TH,TH)] = 1.0 -2.0*phi*zsq/(rsq) ;
1264  gcovpert[TH] =-2.0*phi*zsq/(rsq) ; // order v^4
1265 
1266  gcov[GIND(TH,PH)] = 0.0 ;
1267 
1268  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
1269  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
1270  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
1271  gcov[GIND(PH,PH)] = Rsq;
1272  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0; // doesn't matter that -1.0 subtracted
1273 
1274  // dualfprintf(fail_file,"Mass=%g r=%g R=%g z=%g phi=%g gtt=%g\n",Mass,r,R,z,phi,gcov[GIND(TT,TT)]);
1275  // int jj,kk;
1276  // DLOOP(jj,kk) dualfprintf(fail_file,"jj=%d kk=%d g=%g\n",jj,kk,gcov[GIND(jj,kk)]);
1277 
1278 }
1279 
1280 
1281 
1282 
1283 
1284 
1286 void set_gcov_spcminkmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1287 {
1288  FTYPE r,th;
1289  FTYPE r2small;
1290 #if(SMOOTHSING)
1291  FTYPE signr,rsmooth;
1292 #endif
1293  FTYPE rsharp;
1294 
1295  r = rsharp = V[1];
1296 #if(SMOOTHSING)
1297  signr=mysign(r);
1298  rsmooth = signr*(fabs(r)+SMALL+drsing);
1299  r = rsmooth;
1300  r2small = r*r;
1301 #else
1302  r2small = r*r + SMALL;
1303 #endif
1304 
1305 
1306 
1307 
1308  th=V[2];
1309 
1310  if(POSDEFMETRIC){
1311  if(th<0.0){ th=-th;}
1312  if(th>M_PI) { th=M_PI-th; }
1313  }
1314  else{
1315  }
1316 
1317  // avoid singularity at polar axis
1318 #if(COORDSINGFIX) // just for metric components // hack
1319  if(fabs(th)<SINGSMALL){
1320  if(th>=0) th=SINGSMALL;
1321  if(th<0) th=-SINGSMALL;
1322  }
1323  if(fabs(M_PI-th)<SINGSMALL){
1324  if(th>=M_PI) th=M_PI+SINGSMALL;
1325  if(th<M_PI) th=M_PI-SINGSMALL;
1326  }
1327 #endif
1328 
1329 
1330 
1331  gcov[GIND(TT,TT)] = -1.0;
1332  gcovpert[TT] = 0.0;
1333 
1334  gcov[GIND(TT,RR)] = 0.0 ;
1335  gcov[GIND(TT,TH)] = 0.0 ;
1336  gcov[GIND(TT,PH)] = 0.0 ;
1337 
1338  gcov[GIND(RR,TT)] = 0.0 ;
1339  gcov[GIND(RR,RR)] = 1.0;
1340  gcovpert[RR] = 0.0;
1341 
1342  gcov[GIND(RR,TH)] = 0.0 ;
1343  gcov[GIND(RR,PH)] = 0.0 ;
1344 
1345  gcov[GIND(TH,TT)] = gcov[GIND(TT,TH)] ;
1346  gcov[GIND(TH,RR)] = gcov[GIND(RR,TH)] ;
1347  gcov[GIND(TH,TH)] = rsharp*rsharp ;
1348  gcovpert[TH] = gcov[GIND(TH,TH)]-1.0;
1349 
1350  gcov[GIND(TH,PH)] = 0.0 ;
1351 
1352  gcov[GIND(PH,TT)] = gcov[GIND(TT,PH)] ;
1353  gcov[GIND(PH,RR)] = gcov[GIND(RR,PH)] ;
1354  gcov[GIND(PH,TH)] = gcov[GIND(TH,PH)] ;
1355 
1356  gcov[GIND(PH,PH)] = (rsharp*sin(th))*(rsharp*sin(th));
1357  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0;
1358 
1359 }
1360 
1361 
1365 void set_gcov_ks_bh_tov_metric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1366 {
1367  FTYPE r,th;
1368  FTYPE mysin(FTYPE th);
1369  FTYPE mycos(FTYPE th);
1370  FTYPE MSQ;
1371  FTYPE phi;
1372  int j,k;
1373  void set_gcov_ksmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert);
1374  extern int set_gcov_ks_tov_spcmetric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert, SFTYPE *MOrself, SFTYPE *phiself, SFTYPE *vrsqself);
1375  FTYPE gcovtovks[SYMMATRIXNDIM],gcovtovkspert[NDIM];
1376  FTYPE gcovbhks[SYMMATRIXNDIM],gcovbhkspert[NDIM];
1377  FTYPE MtotOr,MBHprime,MBHprimeOr;
1378  SFTYPE MOrself,phiself,vrsqself;
1379  FTYPE starpot, totalpot;
1380  FTYPE r2small;
1381  FTYPE disc;
1382 #if(SMOOTHSING)
1383  FTYPE signr,rsmooth;
1384 #endif
1385  FTYPE rsharp;
1386 
1387 
1388 
1389 #if(0)
1390  // debug test
1391  DLOOP(j,k) gcov[GIND(j,k)] = 0.0;
1392  DLOOPA(j,j) gcov[GIND(j,j)] = 1.0;
1393  gcovpert[TT] gcov[GIND(TT,TT)]=-1.0;
1394  DLOOPA(j) gcovpert[j]= 0.0;
1395 #endif
1396 
1397 
1398 
1399 
1400 
1401  r = rsharp = V[1];
1402 #if(SMOOTHSING)
1403  signr=mysign(r);
1404  rsmooth = signr*(fabs(r)+SMALL+drsing);
1405  r = rsmooth;
1406  r2small = r*r;
1407 #else
1408  r2small = r*r + SMALL;
1409 #endif
1410 
1411 
1412  th=V[2];
1413 
1414 
1415  // get TOV ks metric assuming pre-existing phivsr_tot, MvsrOr_tot, vrsqvsr_tot
1416  // needs X for interpolation
1417  set_gcov_ks_tov_spcmetric(X, V, gcovtovks, gcovtovkspert, &MOrself, &phiself, &vrsqself);
1418 
1419  // get black hole KS metric
1420  set_gcov_ksmetric(V, gcovbhks, gcovbhkspert);
1421 
1422  // default is KS BH metric
1423  DLOOP(j,k) gcov[GIND(j,k)] = gcovbhks[GIND(j,k)];
1424  DLOOPA(j) gcovpert[j] = gcovbhkspert[j];
1425 
1426  // now compute TOV modifications
1427  // MBHprime = 2.0*MBH*r*r/(a*a+2.0*r*r+a*a*cos(2.0*th));
1428  MBHprime = MBH/(1.0 + (a*a + a*a*cos(2.0*th))/(2.0*r2small));
1429  MBHprimeOr = MBHprime/r;
1430  MtotOr = MBHprimeOr + MOrself; // assume MOrself has sign built-in
1431  starpot=vrsqself - 2.0*MOrself;
1432  totalpot = -2.0*MtotOr;
1433 
1434  // assign mixed metric components
1435  gcov[GIND(TT,TT)] = gcovbhks[GIND(TT,TT)]*(-gcovtovks[GIND(TT,TT)]);
1436  gcovpert[TT] = gcov[GIND(TT,TT)] + 1.0; // no obvious way to remove perturbed part
1437 
1438  // gcov[GIND(RR,RR)] = 1.0 + 2.0*MtotOr-vrsqself;
1439  gcov[GIND(RR,RR)] = 1.0 + fabs(- totalpot);
1440  gcovpert[RR] = fabs(-totalpot);
1441 
1442  disc=-gcovtovks[GIND(TT,TT)]/(1.0-fabs(starpot));
1443  if(disc<0.0){
1444  disc=0.0; // assume this value never used
1445  // and can't trust that signature of metric will work out now, so artificially force signature to be negative
1446  gcov[GIND(TT,TT)]=-fabs(gcov[GIND(TT,TT)]); // should never be used
1447  }
1448 
1449  gcov[GIND(TT,RR)] = gcov[GIND(RR,TT)] = (-totalpot)*sqrt(disc);
1450 
1451  // DEBUG (overwrite)
1452  // DLOOP(j,k) gcov[GIND(j,k)] = gcovtovks[GIND(j,k)];
1453  // DLOOPA(j) gcovpert[j] = gcovtovkspert[j];
1454 
1455  // DLOOP(j,k) dualfprintf(fail_file,"gcovtovks[%d][%d]=%21.15g : %21.15g %21.15g %21.15g\n",j,k,gcovtovks[GIND(j,k)],MOrself,phiself,vrsqself);
1456 
1457 
1458 }
1459 
1460 
1462 void set_gcov_ks_tov_metric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1463 {
1464  extern int set_gcov_ks_tov_spcmetric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert, SFTYPE *MOrself, SFTYPE *phiself, SFTYPE *vrsqself);
1465  SFTYPE MOrself,phiself,vrsqself;
1466 
1467 
1468  // get TOV ks metric assuming pre-existing phivsr_tot, Mvsr_tot, vrsqvsr_tot
1469  // needs X for interpolation
1470  set_gcov_ks_tov_spcmetric(X, V, gcov, gcovpert, &MOrself, &phiself, &vrsqself);
1471 
1472 }
1473 
1475 void set_gcov_bl_tov_metric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1476 {
1477  extern int set_gcov_bl_tov_spcmetric(FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert, SFTYPE *MOrself, SFTYPE *phiself, SFTYPE *vrsqself);
1478  SFTYPE MOrself,phiself,vrsqself;
1479 
1480 
1481  // get TOV BL metric assuming pre-existing phivsr_tot, Mvsr_tot, vrsqvsr_tot
1482  // needs X for interpolation
1483  set_gcov_bl_tov_spcmetric(X, V, gcov, gcovpert, &MOrself, &phiself, &vrsqself);
1484 
1485 }
1486 
1487 
1488 
1489 
1490 
1492 void set_gcov_ksmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1493 {
1494  FTYPE sth, cth, s2, a2, r2, r2small, r3;
1495  FTYPE rho2,rho2small;
1496  FTYPE r,th;
1497  FTYPE mysin(FTYPE th);
1498  FTYPE mycos(FTYPE th);
1499  FTYPE MSQ;
1500  FTYPE phi;
1501  int j,k;
1502 #if(SMOOTHSING)
1503  FTYPE signr,rsmooth;
1504 #endif
1505  FTYPE rsharp;
1506 
1507 
1508 
1509  r = rsharp = V[1];
1510 #if(SMOOTHSING)
1511  signr=mysign(r);
1512  rsmooth = signr*(fabs(r)+SMALL+drsing);
1513  r = rsmooth;
1514  r2small = r*r;
1515  // dualfprintf(fail_file,"got here r=%21.15g drsing=%21.15g\n",r,drsing);
1516 #else
1517  r2small = r*r + SMALL;
1518 #endif
1519 
1520 
1521  // dualfprintf(fail_file,"r=%21.15g\n",r);
1522 
1523  // theta gives well-defined metric for all th
1524  th=V[2];
1525 
1526 
1527  sth=mysin(th);
1528  cth=mycos(th);
1529 
1530 
1531  s2 = sth * sth;
1532  a2 = a * a;
1533  r2 = rsharp * rsharp;
1534  r3 = r2 * r;
1535  rho2 = r2 + a * a * cth * cth;
1536  rho2small = r2small + a * a * cth * cth;
1537  MSQ=(2.*MBH*r-QBH*QBH);
1538 
1539  // really gcov00 diverges at r=0, but force non-divergence to avoid nan's only
1540 #define ks_gcov00 (-1. + MSQ/rho2small)
1541 #define ks_gcov01 (MSQ/rho2small)
1542 #define ks_gcov02 (0)
1543 #define ks_gcov03 (-MSQ*a*s2/rho2small)
1544 #define ks_gcov10 (ks_gcov01)
1545 #define ks_gcov11 (1. + MSQ/rho2small)
1546 #define ks_gcov12 (0)
1547 #define ks_gcov13 (-a*s2*(1. + MSQ/rho2small))
1548 #define ks_gcov20 (0)
1549 #define ks_gcov21 (0)
1550 #define ks_gcov22 (rho2)
1551 #define ks_gcov23 (0)
1552 #define ks_gcov30 (ks_gcov03)
1553 #define ks_gcov31 (ks_gcov13)
1554 #define ks_gcov32 (0)
1555  // old
1556  //#define ks_gcov33 (s2*(rho2 + a*a*s2*(1. + 2.*r/rho2small)))
1557  // new
1558 #define ks_gcov33 (s2*(rho2 + a*a*s2*(1. + MSQ/rho2small)))
1559  // Living review format
1560  //#define ks_gcov33 (s2*(r*r+a*a + a*a*s2*(MSQ/rho2small)))
1561 
1562 
1563  gcov[GIND(TT,TT)] = ks_gcov00 ;
1564  gcovpert[TT] = MSQ/rho2small;
1565 
1566 
1567  gcov[GIND(TT,RR)] = ks_gcov01 ;
1568  gcov[GIND(TT,TH)] = ks_gcov02 ;
1569  gcov[GIND(TT,PH)] = ks_gcov03 ;
1570 
1571  gcov[GIND(RR,TT)] = ks_gcov10 ;
1572 
1573  gcov[GIND(RR,RR)] = ks_gcov11 ;
1574  gcovpert[RR] = MSQ/rho2small ;
1575 
1576  gcov[GIND(RR,TH)] = ks_gcov12 ;
1577  gcov[GIND(RR,PH)] = ks_gcov13 ;
1578 
1579  gcov[GIND(TH,TT)] = ks_gcov20 ;
1580  gcov[GIND(TH,RR)] = ks_gcov21 ;
1581 
1582  gcov[GIND(TH,TH)] = ks_gcov22 ;
1583  gcovpert[TH] = gcov[GIND(TH,TH)]-1.0;
1584 
1585  gcov[GIND(TH,PH)] = ks_gcov23 ;
1586 
1587  gcov[GIND(PH,TT)] = ks_gcov30 ;
1588  gcov[GIND(PH,RR)] = ks_gcov31 ;
1589  gcov[GIND(PH,TH)] = ks_gcov32 ;
1590 
1591  gcov[GIND(PH,PH)] = ks_gcov33 ;
1592  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0;
1593 
1594 
1595  // if(nstep>1051){
1596  // dualfprintf(fail_file,"MBH=%21.15g a=%21.15g gcov[GIND(TT,TT)]=%21.15g\n",MBH,a,gcov[GIND(TT,TT)]);
1597 
1598  // DLOOP(j,k) dualfprintf(fail_file,"ks gcov[%d][%d]=%g\n",j,k,gcov[GIND(j,k)]);
1599  //}
1600 
1601 }
1602 
1603 
1605 void set_gcov_ks_jp1_metric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1606 {
1607  FTYPE r,th;
1608  FTYPE mysin(FTYPE th);
1609  FTYPE mycos(FTYPE th);
1610  FTYPE MSQ;
1611  FTYPE phi;
1612  int j,k;
1613  void set_gcov_ksmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert);
1614  FTYPE gcovksjp1pert[SYMMATRIXNDIM];
1615  FTYPE gcovbhks[SYMMATRIXNDIM],gcovbhkspert[NDIM];
1616  FTYPE r2small;
1617  FTYPE disc;
1618 #if(SMOOTHSING)
1619  FTYPE signr,rsmooth;
1620 #endif
1621  FTYPE rsharp;
1622 
1623 
1624 #if(0)
1625  // debug test
1626  DLOOP(j,k) gcov[GIND(j,k)] = 0.0;
1627  DLOOP(j,k) gcovksjp1pert[GIND(j,k)] = 0.0;
1628  DLOOPA(j,j) gcov[GIND(j,j)] = 1.0;
1629  gcovpert[TT] gcov[GIND(TT,TT)]=-1.0;
1630  DLOOPA(j) gcovpert[j]= 0.0;
1631 #endif
1632 
1633 
1634  r = rsharp = V[1];
1635 #if(SMOOTHSING)
1636  signr=mysign(r);
1637  rsmooth = signr*(fabs(r)+SMALL+drsing);
1638  r = rsmooth;
1639  r2small = r*r;
1640 #else
1641  r2small = r*r + SMALL;
1642 #endif
1643 
1644 
1645  th=V[2];
1646 
1647 
1648  // Difference between Johannsen & Psaltis 2011 metric (/data/jon/mathematica/mathematica_merged/timj_kerrschild_form2.nb)
1649  // See also checkinghowtimjmetricisdifferent.nb for printing of the below
1650  // only non-zero and non-repeated (due to symmetry) terms:
1651  DLOOP(j,k) gcovksjp1pert[GIND(j,k)] = 0.0;
1652 
1653  gcovksjp1pert[GIND(0,0)]=-4.*EP3*r*(2.*r*(-2.*MBH + r) + pow(a,2) + cos(2.*th)*pow(a,2))*pow(MBH,3)*pow(pow(a,2) + cos(2.*th)*pow(a,2) + 2.*pow(r,2),-3);
1654 
1655  gcovksjp1pert[GIND(0,1)]=2.*EP3*pow(MBH,4)*pow(r,2)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-3);
1656 
1657  gcovksjp1pert[GIND(1,1)]=-1. - 4.*pow(MBH,2)*pow(r,2)*(EP3*r*pow(MBH,3) + pow(r,4) + 2.*pow(a,2)*pow(r,2)*pow(cos(th),2) + pow(a,4)*pow(cos(th),4))*pow(r*(-2.*MBH + r) + pow(a,2),-1)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-3) - 2.*MBH*r*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-1) + (EP3*r*pow(MBH,3)*(pow(r,2) + pow(a,2)*pow(cos(th),2)) + pow(pow(r,2) + pow(a,2)*pow(cos(th),2),3))*pow((r*(-2.*MBH + r) + pow(a,2))*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),2) + EP3*r*pow(a,2)*pow(MBH,3)*pow(sin(th),2),-1) + pow(a,2)*pow(r*(-2.*MBH + r) + pow(a,2),-2)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-3)*pow(sin(th),2)*(-4.*EP3*pow(MBH,5)*pow(r,3) - 4.*pow(MBH,2)*pow(r,2)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),2) + (pow(a,2) + pow(r,2))*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),3) + MBH*r*pow(a,2)*(2.*EP3*r*pow(MBH,3) + EP3*pow(MBH,2)*(pow(r,2) + pow(a,2)*pow(cos(th),2)) + 2.*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),2))*pow(sin(th),2));
1658 
1659  gcovksjp1pert[GIND(0,3)]=-2.*a*EP3*pow(MBH,4)*pow(r,2)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-3)*pow(sin(th),2);
1660 
1661  gcovksjp1pert[GIND(1,3)]=0.125*a*EP3*r*pow(MBH,3)*(-4.*r*(2.*MBH + r)*pow(a,2) + 4.*r*(2.*MBH + r)*cos(2.*th)*pow(a,2) - 1.*pow(a,4) + cos(4.*th)*pow(a,4) + 32.*pow(MBH,2)*pow(r,2))*pow(r*(-2.*MBH + r) + pow(a,2),-1)*pow(pow(r,2) + pow(a,2)*pow(cos(th),2),-3)*pow(sin(th),2);
1662 
1663  gcovksjp1pert[GIND(3,3)]=4.*EP3*r*pow(a,2)*(2.*r*(2.*MBH + r) + pow(a,2) + cos(2.*th)*pow(a,2))*pow(MBH,3)*pow(pow(a,2) + cos(2.*th)*pow(a,2) + 2.*pow(r,2),-3)*pow(sin(th),4);
1664 
1665 
1666 
1667  // get black hole KS metric
1668  set_gcov_ksmetric(V, gcovbhks, gcovbhkspert);
1669 
1670 
1671  // default is KS BH metric
1672  // and add "perturbation" that is JP1 metric
1673  // NOTEMARK: when using DLOOP(j,k) over GIND(j,k) for metric, *cannot* use += since that would go over same value twice!
1674  DLOOP(j,k) gcov[GIND(j,k)] = gcovbhks[GIND(j,k)] + gcovksjp1pert[GIND(j,k)];
1675  DLOOPA(j) gcovpert[j] = gcovbhkspert[j] + gcovksjp1pert[GIND(j,j)];
1676  // i.e. JP1 terms are deviation from Kerr. NOTEMARK: For g_{rr}, this adds and subtracts unity, so not acurate machine accurate. But that's ok since applications with this metric are always for near the BH.
1677 
1678 
1679 
1680 }
1681 
1682 
1683 
1684 
1685 
1687 void set_gcov_blmetric(FTYPE *V, FTYPE *gcov, FTYPE *gcovpert)
1688 {
1689  FTYPE sth, cth, s2, a2, r2, r2small, r3, DD, mu;
1690  FTYPE mupert,DDpert;
1691  FTYPE r,th;
1692  FTYPE MSQ;
1693  int j,k;
1694 #if(SMOOTHSING)
1695  FTYPE signr,rsmooth;
1696 #endif
1697  FTYPE rsharp;
1698 
1699  r = rsharp = V[1];
1700 #if(SMOOTHSING)
1701  signr=mysign(r);
1702  rsmooth = signr*(fabs(r)+SMALL+drsing);
1703  r = rsmooth;
1704  r2small = r*r;
1705 #else
1706  r2small = r*r + SMALL;
1707 #endif
1708 
1709 
1710 
1711  th=V[2];
1712 
1713 
1714 
1715 
1716  cth = cos(th);
1717  sth=sin(th);
1718 
1719 
1720  MSQ=(2.*MBH*rsharp-QBH*QBH);
1721 
1722  s2 = sth * sth;
1723  a2 = a * a;
1724  r2 = rsharp * rsharp;
1725  r3 = r2 * r;
1726 
1727  DDpert =- MSQ / r2small + a2 / r2small;
1728  DD = 1. + DDpert;
1729  mupert=+ a2 * cth * cth / r2small;
1730  mu = 1. + mupert;
1731 
1732 
1733 
1734 #define bl_gcov00 (-(1. - MSQ/(r2small*mu)))
1735 #define bl_gcov01 (0)
1736 #define bl_gcov02 (0)
1737 #define bl_gcov03 (-MSQ*a*s2/(r2small*mu))
1738 #define bl_gcov10 (0)
1739 #define bl_gcov11 (mu/DD)
1740 #define bl_gcov12 (0)
1741 #define bl_gcov13 (0)
1742 #define bl_gcov20 (0)
1743 #define bl_gcov21 (0)
1744 #define bl_gcov22 (r2*mu)
1745 #define bl_gcov23 (0)
1746 #define bl_gcov30 (bl_gcov03)
1747 #define bl_gcov31 (0)
1748 #define bl_gcov32 (0)
1749 #define bl_gcov33 (r2*sth*sth*(1. + a2/r2small + MSQ*a2*s2/(r2small*r2small*mu)))
1750 
1751  gcov[GIND(TT,TT)] = bl_gcov00 ;
1752  gcovpert[TT] = MSQ/(r2small*mu);
1753 
1754  gcov[GIND(TT,RR)] = bl_gcov01 ;
1755  gcov[GIND(TT,TH)] = bl_gcov02 ;
1756  gcov[GIND(TT,PH)] = bl_gcov03 ;
1757 
1758  gcov[GIND(RR,TT)] = bl_gcov10 ;
1759 
1760  gcov[GIND(RR,RR)] = bl_gcov11 ;
1761  // mu/DD = (1.0+mupert)/(1.0+DDpert) =
1762  gcovpert[RR] = (DDpert - mupert) / (1.0 + DDpert);
1763 
1764  gcov[GIND(RR,TH)] = bl_gcov12 ;
1765  gcov[GIND(RR,PH)] = bl_gcov13 ;
1766 
1767  gcov[GIND(TH,TT)] = bl_gcov20 ;
1768  gcov[GIND(TH,RR)] = bl_gcov21 ;
1769 
1770  gcov[GIND(TH,TH)] = bl_gcov22 ;
1771  gcovpert[TH] = gcov[GIND(TH,TH)]-1.0;
1772 
1773  gcov[GIND(TH,PH)] = bl_gcov23 ;
1774 
1775  gcov[GIND(PH,TT)] = bl_gcov30 ;
1776  gcov[GIND(PH,RR)] = bl_gcov31 ;
1777  gcov[GIND(PH,TH)] = bl_gcov32 ;
1778 
1779  gcov[GIND(PH,PH)] = bl_gcov33 ;
1780  gcovpert[PH] = gcov[GIND(PH,PH)]-1.0;
1781 
1782 
1783 
1784  // DLOOP(j,k) dualfprintf(fail_file,"bl gcov[%d][%d]=%g\n",j,k,gcov[GIND(j,k)]);
1785 
1786 
1787 
1788 }
1789 
1790 
1797 void set_gcon_ksmetric(FTYPE *V, FTYPE *gcon)
1798 {
1799  FTYPE sth, cth, s2, a2, r2, r2small,r3;
1800  FTYPE rho2;
1801  FTYPE r,th;
1802  FTYPE MSQ;
1803 #if(SMOOTHSING)
1804  FTYPE signr,rsmooth;
1805 #endif
1806  FTYPE rsharp;
1807 
1808  r = rsharp = V[1];
1809 #if(SMOOTHSING)
1810  signr=mysign(r);
1811  rsmooth = signr*(fabs(r)+SMALL+drsing);
1812  r = rsmooth;
1813  r2small = r*r;
1814 #else
1815  r2small = r*r + SMALL;
1816 #endif
1817  // size of r and M controlled by changing M, not rescaling r
1818 
1819 
1820 
1821 
1822  th=V[2];
1823 
1824  cth = cos(th);
1825  sth=sin(th);
1826 
1827 
1828  s2 = sth * sth;
1829  a2 = a * a;
1830  r2 = r * r;
1831  r3 = r2 * r;
1832  rho2 = r2small + a * a * cth * cth; // always divide by this
1833  MSQ=(2.*MBH*r-QBH*QBH);
1834 
1835 
1836 #define ks_gcon00 (-(1.+ MSQ/rho2))
1837 #define ks_gcon01 (MSQ/rho2)
1838 #define ks_gcon02 (0)
1839 #define ks_gcon03 (0)
1840 #define ks_gcon10 (ks_gcon01)
1841  //#define ks_gcon11 ((r*(r-2.)+a*a)/rho2)
1842 #define ks_gcon11 ((r*r+a*a-MSQ)/rho2)
1843 #define ks_gcon12 (0)
1844 #define ks_gcon13 (a/rho2)
1845 #define ks_gcon20 (ks_gcon02)
1846 #define ks_gcon21 (ks_gcon12)
1847 #define ks_gcon22 (1./rho2)
1848 #define ks_gcon23 (0)
1849 #define ks_gcon30 (ks_gcon03)
1850 #define ks_gcon31 (ks_gcon13)
1851 #define ks_gcon32 (ks_gcon23)
1852 #define ks_gcon33 (1./(rho2*s2))
1853 
1854 
1855  gcon[GIND(TT,TT)] = ks_gcon00 ;
1856  gcon[GIND(TT,RR)] = ks_gcon01 ;
1857  gcon[GIND(TT,TH)] = ks_gcon02 ;
1858  gcon[GIND(TT,PH)] = ks_gcon03 ;
1859 
1860  gcon[GIND(RR,TT)] = ks_gcon10 ;
1861  gcon[GIND(RR,RR)] = ks_gcon11 ;
1862  gcon[GIND(RR,TH)] = ks_gcon12 ;
1863  gcon[GIND(RR,PH)] = ks_gcon13 ;
1864 
1865  gcon[GIND(TH,TT)] = ks_gcon20 ;
1866  gcon[GIND(TH,RR)] = ks_gcon21 ;
1867  gcon[GIND(TH,TH)] = ks_gcon22 ;
1868  gcon[GIND(TH,PH)] = ks_gcon23 ;
1869 
1870  gcon[GIND(PH,TT)] = ks_gcon30 ;
1871  gcon[GIND(PH,RR)] = ks_gcon31 ;
1872  gcon[GIND(PH,TH)] = ks_gcon32 ;
1873  if(s2!=0.0){
1874  gcon[GIND(PH,PH)] = ks_gcon33 ;
1875  }
1876  else gcon[GIND(PH,PH)]=1.0; // avoid coordinate singularity -- although should never use this value
1877 
1878 }
1879 
1881 void set_gcon_blmetric(FTYPE *V, FTYPE *gcon)
1882 {
1883  FTYPE sth, cth, s2, a2, r2, r2small, r3, DD, mu;
1884  FTYPE r,th;
1885  FTYPE MSQ;
1886 #if(SMOOTHSING)
1887  FTYPE signr,rsmooth;
1888 #endif
1889  FTYPE rsharp;
1890 
1891  r = rsharp = V[1];
1892 #if(SMOOTHSING)
1893  signr=mysign(r);
1894  rsmooth = signr*(fabs(r)+SMALL+drsing);
1895  r = rsmooth;
1896  r2small = r*r;
1897 #else
1898  r2small = r*r + SMALL;
1899 #endif
1900 
1901 
1902  // size of r relative to M is controlled by M, not r.
1903 
1904  th=V[2];
1905 
1906  cth = cos(th);
1907  sth=sin(th);
1908 
1909  MSQ=(2.*MBH*r-QBH*QBH);
1910 
1911  s2 = sth * sth;
1912  a2 = a * a;
1913  r2 = r * r;
1914  r3 = r2 * r;
1915  DD = 1. - MSQ / r2 + a2 / r2;
1916  mu = 1. + a2 * cth * cth / r2;
1917 
1918 
1919 #define bl_gcon00 (-1. - MSQ*(1. + a2/r2small)/(r2small*DD*mu))
1920 #define bl_gcon01 (0)
1921 #define bl_gcon02 (0)
1922 #define bl_gcon03 (-MSQ*a/(r2small*r2small*DD*mu))
1923 #define bl_gcon10 (0)
1924 #define bl_gcon11 (DD/mu)
1925 #define bl_gcon12 (0)
1926 #define bl_gcon13 (0)
1927 #define bl_gcon20 (0)
1928 #define bl_gcon21 (0)
1929 #define bl_gcon22 (1./(r2small*mu))
1930 #define bl_gcon23 (0)
1931 #define bl_gcon30 (bl_gcon03)
1932 #define bl_gcon31 (0)
1933 #define bl_gcon32 (0)
1934 #define bl_gcon33 ((1. - MSQ/(r2small*mu))/(r2small*sth*sth*DD))
1935 
1936 
1937  gcon[GIND(TT,TT)] = bl_gcon00 ;
1938  gcon[GIND(TT,RR)] = bl_gcon01 ;
1939  gcon[GIND(TT,TH)] = bl_gcon02 ;
1940  gcon[GIND(TT,PH)] = bl_gcon03 ;
1941 
1942  gcon[GIND(RR,TT)] = bl_gcon10 ;
1943  gcon[GIND(RR,RR)] = bl_gcon11 ;
1944  gcon[GIND(RR,TH)] = bl_gcon12 ;
1945  gcon[GIND(RR,PH)] = bl_gcon13 ;
1946 
1947  gcon[GIND(TH,TT)] = bl_gcon20 ;
1948  gcon[GIND(TH,RR)] = bl_gcon21 ;
1949  gcon[GIND(TH,TH)] = bl_gcon22 ;
1950  gcon[GIND(TH,PH)] = bl_gcon23 ;
1951 
1952  gcon[GIND(PH,TT)] = bl_gcon30 ;
1953  gcon[GIND(PH,RR)] = bl_gcon31 ;
1954  gcon[GIND(PH,TH)] = bl_gcon32 ;
1955  gcon[GIND(PH,PH)] = bl_gcon33 ;
1956 
1957 }
1958 
1964 void set_conn_general(FTYPE *X, struct of_geom *geom,
1965  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
1966 {
1967  void conn_func_numerical1(FTYPE DELTA,FTYPE *X, struct of_geom *geom,
1968  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
1969 
1970  conn_func_numerical1(CONNDELTA,X, geom, conn, conn2);
1971 }
1972 
1973 
1974 
1975 void set_conn_cylminkmetric(FTYPE *X, struct of_geom *geom,
1976  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
1977 {
1978  void conn_func_numerical1(FTYPE DELTA,FTYPE *X, struct of_geom *geom,
1979  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
1980  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
1981  int i,j,k;
1982  FTYPE gcovmid[SYMMATRIXNDIM];
1983  FTYPE gcovpertmid[NDIM];
1984  FTYPE gdetmid;
1985 
1986  FTYPE dxdxp[NDIM][NDIM]; //atch
1987  FTYPE V[NDIM]; //atch
1988 
1989  if(ANALYTICCONNECTION&&defcoord==UNIFORMCOORDS){ // uniform grid SUPERSASMARK
1990  // could directly use gdet in global memory
1991  // only works for X1=R and X2=z
1992  gcov_func(geom, 1,CYLMINKMETRIC,X, gcovmid,gcovpertmid);
1993 
1994  bl_coord( X, V ); //actually, dxdxprim() does not use V or X since metric is uniform
1995  if(gdet_func_metric(MCOORD,V,gcovmid,&gdetmid)!=0){
1996  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() problem in set_conn_cylminkmetric()\n");
1997  }
1998 
1999 
2000 
2001  // see transforms.c and mettometp() and see gcov2gcovprim() //atch
2002  dxdxprim(X, V, dxdxp); //atch
2003 
2004  for (k = 0; k < NDIM; k++) conn2[k]= 0.0;
2005  //conn2[RR]=-1.0/gdetmid; //wrong as well... shouldn't it be 0 since there is no 2nd connection when WITHGDET? SUPERSASMARK
2006 
2007 
2008  for (i = 0; i < NDIM; i++)
2009  for (j = 0; j < NDIM; j++)
2010  for (k = 0; k < NDIM; k++) {
2011  conn[i][j][k] = 0.;
2012  }
2013  conn[PH][RR][PH]=1.0 / X[1]; //1.0/gdetmid; //apparently, wrong because should not care about the 2nd dimension SUPERSASMARK
2014  conn[PH][PH][RR]=1.0 / X[1]; //1.0/gdetmid; //apparently, wrong because should not care about the 2nd dimension SUPERSASMARK
2015  conn[RR][PH][PH]= - X[1] * dxdxp[3][3] * dxdxp[3][3]; //-gdetmid;
2016  }
2017  else{
2018  conn_func_numerical1(CONNDELTA,X, geom, conn, conn2);
2019  }
2020 
2021 
2022 
2023 }
2024 
2026 void set_conn_cartminkmetric(FTYPE *X, struct of_geom *geom,
2027  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
2028 {
2029  void conn_func_numerical1(FTYPE DELTA,FTYPE *X, struct of_geom *geom,
2030  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
2031 
2032  int i,j,k;
2033 
2034 
2035  if(ANALYTICCONNECTION&&defcoord==UNIFORMCOORDS){// uniform grid
2036  for (k = 0; k < NDIM; k++) conn2[k]= 0.0;
2037 
2038  for (i = 0; i < NDIM; i++)
2039  for (j = 0; j < NDIM; j++)
2040  for (k = 0; k < NDIM; k++) {
2041  conn[i][j][k] = 0.;
2042  }
2043  }
2044  else{
2045  conn_func_numerical1(CONNDELTA,X, geom, conn, conn2);
2046  }
2047 
2048 }
2049 
2050 
2051 
2052 void set_conn_ksmetric(FTYPE *X, struct of_geom *geom,
2053  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
2054 {
2055  void conn_func_numerical1(FTYPE DELTA,FTYPE *X, struct of_geom *geom,
2056  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2);
2057  void mks_conn_func(FTYPE *X, struct of_geom *geom,
2058  FTYPE (*lconn)[NDIM][NDIM],FTYPE *conn2);
2059 
2060  // FTYPE DELTA; // for debug below
2061 
2062  // the analytic form can be used with WHICHEOM=WITHNOGDET
2063  // determine for which we have analytic expressions
2064  // this currently competes with mks_source_conn in phys.c
2065  if((!ANALYTICCONNECTION)||(defcoord!=LOGRSINTH)) conn_func_numerical1(CONNDELTA,X, geom, conn, conn2);
2066  else mks_conn_func(X,geom,conn,conn2);
2067 
2068 
2069  /*
2070  // some debug stuff
2071  // if((geom->i==62)&&(geom->j==32)){ // where c002 is bad
2072  if(0&&(geom->i==32)&&(geom->j==0)){ // where c023 is bad
2073  for(DELTA=1E-15;DELTA<1E5;DELTA*=1.01){
2074  conn_func_numerical1(DELTA,X, geom, conn, conn2);
2075  dualfprintf(fail_file,"%30.20Lg %30.20Lg\n",DELTA,conn[0][2][3]); fflush(fail_file);
2076  // dualfprintf(fail_file,"%21.15g %21.15g\n",DELTA,conn[0][2][3]); fflush(fail_file);
2077  }
2078  exit(0);
2079  }
2080  else{
2081  conn_func_numerical1(CONNDELTA,X, geom, conn, conn2);
2082  }
2083  //mks_conn_func(X,geom,conn,conn2);
2084  */
2085 
2086 }
2087 
2088 
2089 
2090 
2091 
2092 
2093 
2094 
2095 
2096 
2097 
2098 
2099 
2100 
2101 
2102 
2103 
2104 
2105 
2106 
2107 
2108 
2109 
2110 
2111 
2112 
2113 
2114 
2115 
2116 
2117 
2118 
2119 
2120 
2121 
2122 
2123 
2124 
2125 
2126 
2127 
2128 
2130 //
2131 // below very specific to defcoord==LOGRSINTH and MCOORD=KSCOORDS:
2132 // gives: analytic connection and source
2133 //
2135 
2136 
2137 
2140 void mks_conn_func(FTYPE *X, struct of_geom *ptrgeom,
2141  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
2142 {
2143  int i, j, k, l;
2144  FTYPE V[NDIM];
2145  FTYPE r,th,sigma,dxdxptrue[NDIM][NDIM];
2146 #ifdef WIN32
2147  extern FTYPE cot(FTYPE arg);
2148 #endif
2149  FTYPE dxdxp[NDIM];
2150 
2151 
2152  if(MBH!=1.0){
2153  dualfprintf(fail_file,"mks_conn_func not setup for MBH!=1.0\n");
2154  myexit(10);
2155  }
2156 
2157 
2158  // get bl coordinates
2159  bl_coord_ijk(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, V);
2160  r=V[1];
2161  th=V[2];
2162 
2163 
2164  // the connection
2165 
2166  // this is not exactly right, since derivative of metric is derivative of absolute values, but shouldn't/doesn't seem to matter much
2167  // follows gcov_func()
2168  if(POSDEFMETRIC){
2169  if(th<0.0){ th=-th;}
2170  if(th>M_PI) { th=M_PI-th; }
2171  }
2172  else{
2173  }
2174 
2175  // avoid singularity at polar axis
2176 #if(COORDSINGFIX)
2177  if(fabs(th)<SINGSMALL){
2178  if(th>=0) th=SINGSMALL;
2179  if(th<0) th=-SINGSMALL;
2180  }
2181  if(fabs(M_PI-th)<SINGSMALL){
2182  if(th>=M_PI) th=M_PI+SINGSMALL;
2183  if(th<M_PI) th=M_PI-SINGSMALL;
2184  }
2185 #endif
2186 
2187  // set aux vars
2188  dxdxprim_ijk(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, dxdxptrue);
2189  DLOOPA(j) dxdxp[j]=dxdxptrue[j][j]; // defcoord==LOGRSINTH assumes transformation is diagonal
2190  sigma=r*r+a*a*cos(th)*cos(th);
2191 
2192 
2193  conn[0][0][0]=(-2.*r*sigma + 4.*pow(r,3.))*pow(sigma,-3.);
2194  conn[0][0][1]=dxdxp[1]*(2.*r + sigma)*(-1.*sigma + 2.*pow(r,2.))*
2195  pow(sigma,-3.);
2196  conn[0][0][2]=-1.*dxdxp[2]*r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th);
2197  conn[0][0][3]=-2.*a*r*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*
2198  pow(sin(th),2.);
2199  conn[0][1][0]=dxdxp[1]*(2.*r + sigma)*(-1.*sigma + 2.*pow(r,2.))*
2200  pow(sigma,-3.);
2201  conn[0][1][1]=2.*(r + sigma)*pow(dxdxp[1],2.)*
2202  (-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.);
2203  conn[0][1][2]=-1.*dxdxp[1]*dxdxp[2]*r*pow(a,2.)*pow(sigma,-2.)*
2204  sin(2.*th);
2205  conn[0][1][3]=dxdxp[1]*a*(2.*r + sigma)*(sigma - 2.*pow(r,2.))*
2206  pow(sigma,-3.)*pow(sin(th),2.);
2207  conn[0][2][0]=-1.*dxdxp[2]*r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th);
2208  conn[0][2][1]=-1.*dxdxp[1]*dxdxp[2]*r*pow(a,2.)*pow(sigma,-2.)*
2209  sin(2.*th);
2210  conn[0][2][2]=-2.*pow(dxdxp[2],2.)*pow(r,2.)*pow(sigma,-1.);
2211  conn[0][2][3]=2.*dxdxp[2]*r*cos(th)*pow(a,3.)*pow(sigma,-2.)*
2212  pow(sin(th),3.);
2213  conn[0][3][0]=-2.*a*r*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*
2214  pow(sin(th),2.);
2215  conn[0][3][1]=dxdxp[1]*a*(2.*r + sigma)*(sigma - 2.*pow(r,2.))*
2216  pow(sigma,-3.)*pow(sin(th),2.);
2217  conn[0][3][2]=2.*dxdxp[2]*r*cos(th)*pow(a,3.)*pow(sigma,-2.)*
2218  pow(sin(th),3.);
2219  conn[0][3][3]=2.*r*pow(sigma,-3.)*pow(sin(th),2.)*
2220  (-1.*r*pow(sigma,2.) + pow(a,2.)*(-1.*sigma + 2.*pow(r,2.))*
2221  pow(sin(th),2.));
2222  conn[1][0][0]=pow(dxdxp[1],-1.)*(-1.*sigma + 2.*pow(r,2.))*
2223  pow(sigma,-3.)*(-2.*r + sigma + pow(a,2.)*pow(sin(th),2.));
2224  conn[1][0][1]=0.5*(4.*r - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
2225  (sigma - 2.*pow(r,2.))*pow(sigma,-3.);
2226  conn[1][0][2]=0.;
2227  conn[1][0][3]=0.5*a*pow(dxdxp[1],-1.)*
2228  (4.*r - 2.*sigma - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
2229  (-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*pow(sin(th),2.);
2230  conn[1][1][0]=0.5*(4.*r - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
2231  (sigma - 2.*pow(r,2.))*pow(sigma,-3.);
2232  conn[1][1][1]=pow(sigma,-3.)*
2233  (-1.*dxdxp[1]*(2.*r + sigma)*(-1.*sigma + 2.*pow(r,2.)) +
2234  pow(sigma,3.) + dxdxp[1]*pow(a,2.)*(-1.*sigma + 2.*pow(r,2.))*
2235  pow(sin(th),2.));
2236  conn[1][1][2]=-1.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*sin(th)
2237  ;
2238  conn[1][1][3]=0.5*a*(pow(a,2.)*(sigma - 2.*pow(r,2.)) +
2239  cos(2.*th)*pow(a,2.)*(-1.*sigma + 2.*pow(r,2.)) +
2240  2.*r*((-2. + sigma)*sigma + 4.*pow(r,2.)))*pow(sigma,-3.)*
2241  pow(sin(th),2.);
2242  conn[1][2][0]=0.;
2243  conn[1][2][1]=-1.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*sin(th)
2244  ;
2245  conn[1][2][2]=-1.*r*pow(dxdxp[1],-1.)*pow(dxdxp[2],2.)*
2246  pow(sigma,-1.)*(-2.*r + sigma + pow(a,2.)*pow(sin(th),2.));
2247  conn[1][2][3]=0.;
2248  conn[1][3][0]=0.5*a*pow(dxdxp[1],-1.)*
2249  (4.*r - 2.*sigma - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
2250  (-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*pow(sin(th),2.);
2251  conn[1][3][1]=0.5*a*(pow(a,2.)*(sigma - 2.*pow(r,2.)) +
2252  cos(2.*th)*pow(a,2.)*(-1.*sigma + 2.*pow(r,2.)) +
2253  2.*r*((-2. + sigma)*sigma + 4.*pow(r,2.)))*pow(sigma,-3.)*
2254  pow(sin(th),2.);
2255  conn[1][3][2]=0.;
2256  conn[1][3][3]=-1.*pow(dxdxp[1],-1.)*pow(sigma,-3.)*pow(sin(th),2.)*
2257  (-2.*r + sigma + pow(a,2.)*pow(sin(th),2.))*
2258  (r*pow(sigma,2.) + pow(a,2.)*(sigma - 2.*pow(r,2.))*pow(sin(th),2.));
2259  conn[2][0][0]=-1.*r*pow(dxdxp[2],-1.)*pow(a,2.)*pow(sigma,-3.)*
2260  sin(2.*th);
2261  conn[2][0][1]=-1.*dxdxp[1]*r*pow(dxdxp[2],-1.)*pow(a,2.)*
2262  pow(sigma,-3.)*sin(2.*th);
2263  conn[2][0][2]=0.;
2264  conn[2][0][3]=2.*a*r*cos(th)*pow(dxdxp[2],-1.)*pow(sigma,-3.)*
2265  (sigma + pow(a,2.)*pow(sin(th),2.))*sin(th);
2266  conn[2][1][0]=-1.*dxdxp[1]*r*pow(dxdxp[2],-1.)*pow(a,2.)*
2267  pow(sigma,-3.)*sin(2.*th);
2268  conn[2][1][1]=-1.*r*pow(dxdxp[1],2.)*pow(dxdxp[2],-1.)*pow(a,2.)*
2269  pow(sigma,-3.)*sin(2.*th);
2270  conn[2][1][2]=dxdxp[1]*r*pow(sigma,-1.);
2271  conn[2][1][3]=dxdxp[1]*a*pow(dxdxp[2],-1.)*pow(sigma,-3.)*sin(th)*
2272  (sigma*(2.*r + sigma)*cos(th) + r*pow(a,2.)*sin(th)*sin(2.*th));
2273  conn[2][2][0]=0.;
2274  conn[2][2][1]=dxdxp[1]*r*pow(sigma,-1.);
2275  conn[2][2][2]=4.*(M_PI*X[2] - 1.*th)*pow(dxdxp[2],-1.)*
2276  pow(M_PI,2.) - 1.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*sin(th)\
2277  ;
2278  conn[2][2][3]=0.;
2279  conn[2][3][0]=2.*a*r*cos(th)*pow(dxdxp[2],-1.)*pow(sigma,-3.)*
2280  (sigma + pow(a,2.)*pow(sin(th),2.))*sin(th);
2281  conn[2][3][1]=dxdxp[1]*a*pow(dxdxp[2],-1.)*pow(sigma,-3.)*sin(th)*
2282  (sigma*(2.*r + sigma)*cos(th) + r*pow(a,2.)*sin(th)*sin(2.*th));
2283  conn[2][3][2]=0.;
2284  conn[2][3][3]=-1.*cos(th)*pow(dxdxp[2],-1.)*pow(sigma,-3.)*
2285  (pow(sigma,3.) + sigma*(4.*r + sigma)*pow(a,2.)*pow(sin(th),2.) +
2286  2.*r*pow(a,4.)*pow(sin(th),4.))*sin(th);
2287  conn[3][0][0]=a*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.);
2288  conn[3][0][1]=dxdxp[1]*a*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)
2289  ;
2290  conn[3][0][2]=-2.*dxdxp[2]*a*r*cot(th)*pow(sigma,-2.);
2291  conn[3][0][3]=-1.*pow(a,2.)*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*
2292  pow(sin(th),2.);
2293  conn[3][1][0]=dxdxp[1]*a*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)
2294  ;
2295  conn[3][1][1]=a*pow(dxdxp[1],2.)*(-1.*sigma + 2.*pow(r,2.))*
2296  pow(sigma,-3.);
2297  conn[3][1][2]=-1.*dxdxp[1]*dxdxp[2]*a*(2.*r + sigma)*cot(th)*
2298  pow(sigma,-2.);
2299  conn[3][1][3]=dxdxp[1]*pow(sigma,-3.)*
2300  (r*pow(sigma,2.) + pow(a,2.)*(sigma - 2.*pow(r,2.))*pow(sin(th),2.));
2301  conn[3][2][0]=-2.*dxdxp[2]*a*r*cot(th)*pow(sigma,-2.);
2302  conn[3][2][1]=-1.*dxdxp[1]*dxdxp[2]*a*(2.*r + sigma)*cot(th)*
2303  pow(sigma,-2.);
2304  conn[3][2][2]=-1.*a*r*pow(dxdxp[2],2.)*pow(sigma,-1.);
2305  conn[3][2][3]=dxdxp[2]*
2306  (cot(th) + r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th));
2307  conn[3][3][0]=-1.*pow(a,2.)*(-1.*sigma + 2.*pow(r,2.))*pow(sigma,-3.)*
2308  pow(sin(th),2.);
2309  conn[3][3][1]=dxdxp[1]*pow(sigma,-3.)*
2310  (r*pow(sigma,2.) + pow(a,2.)*(sigma - 2.*pow(r,2.))*pow(sin(th),2.));
2311  conn[3][3][2]=dxdxp[2]*
2312  (cot(th) + r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th));
2313  conn[3][3][3]=pow(sigma,-3.)*
2314  (-1.*a*r*pow(sigma,2.)*pow(sin(th),2.) +
2315  pow(a,3.)*(-1.*sigma + 2.*pow(r,2.))*pow(sin(th),4.));
2316  conn2[0]=0.;
2317  conn2[1]=-1.*pow(sigma,-1.)*(2.*dxdxp[1]*r + pow(r,2.) +
2318  pow(a,2.)*pow(cos(th),2.));
2319  conn2[2]=-1.*dxdxp[2]*cot(th) +
2320  4.*(-1.*M_PI*X[2] + th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) +
2321  dxdxp[2]*pow(a,2.)*pow(sigma,-1.)*sin(2.*th);
2322  conn2[3]=0.;
2323 
2324 
2325 }
2326 
2327 
2328 
2329 /*********************************************************************************************
2330  Scott's s MKS connection that can be used for any transformation between r,th <-> X1,X2 :
2331 *********************************************************************************************/
2332 #define M (1.)
2333 void mks_conn_func_general(FTYPE *X, struct of_geom *geom, FTYPE (*conn)[NDIM][NDIM] )
2334 {
2335  int i, j, k, l;
2336  FTYPE V[NDIM];
2337  FTYPE r,th,sigma,dxdxp[NDIM][NDIM],dxdxp_dxp[NDIM][NDIM][NDIM];
2338 
2339  FTYPE t1, t10, t102, t1024, t1035, t1037, t104, t11, t114, t116, t119;
2340  FTYPE t12, t121, t123, t126, t129, t130, t132, t14, t148, t149, t15, t152;
2341  FTYPE t154, t156, t157, t158, t159, t161, t169, t17, t171, t172, t175, t177;
2342  FTYPE t185, t2, t203, t204, t208, t209, t21, t210, t212, t214, t22, t221;
2343  FTYPE t222, t224, t227, t23, t236, t24, t240, t241, t242, t243, t245, t246;
2344  FTYPE t247, t248, t25, t250, t251, t258, t26, t260, t261, t263, t264, t271;
2345  FTYPE t273, t275, t276, t278, t28, t280, t281, t283, t284, t285, t286, t288;
2346  FTYPE t289, t29, t297, t299, t3, t30, t300, t303, t305, t306, t308, t309;
2347  FTYPE t31, t310, t313, t314, t320, t325, t327, t328, t329, t330, t333, t336;
2348  FTYPE t338, t34, t340, t342, t344, t346, t35, t358, t361, t363, t366, t367;
2349  FTYPE t368, t370, t372, t375, t38, t380, t381, t384, t385, t387, t39, t399;
2350  FTYPE t4, t40, t400, t402, t404, t405, t406, t408, t409, t41, t411, t412;
2351  FTYPE t418, t42, t421, t425, t428, t431, t432, t433, t434, t437, t440, t442;
2352  FTYPE t448, t451, t453, t454, t459, t462, t467, t469, t480, t481, t486, t487;
2353  FTYPE t488, t491, t492, t498, t501, t504, t507, t508, t510, t512, t52, t521;
2354  FTYPE t528, t53, t530, t553, t556, t56, t57, t588, t60, t607, t627, t628;
2355  FTYPE t63, t630, t631, t632, t634, t636, t637, t64, t651, t652, t654, t656;
2356  FTYPE t657, t659, t661, t662, t670, t673, t675, t677, t686, t689, t7, t712;
2357  FTYPE t74, t748, t75, t78, t793, t794, t795, t799, t8, t800, t801, t803;
2358  FTYPE t806, t807, t813, t816, t822, t83, t831, t84, t845, t86, t89, t891;
2359  FTYPE t90, t91, t916, t917, t920, t924, t928, t940, t946, t968, t97, t970,t991;
2360 
2361 
2362  if(MBH!=1.0){
2363  dualfprintf(fail_file,"mks_conn_func_general not setup for MBH!=1.0\n");
2364  myexit(11);
2365  }
2366 
2367 
2368  // get bl coordinates
2369  bl_coord(X,V);
2370  r=V[1];
2371  th=V[2];
2372 
2373  // the connection
2374 
2375  // this is not exactly right, since derivative of metric is derivative of absolute values, but shouldn't/doesn't seem to matter much
2376  // follows gcov_func()
2377  if(POSDEFMETRIC){
2378  if(th<0.0){ th=-th;}
2379  if(th>M_PI) { th=M_PI-th; }
2380  }
2381  else{
2382  }
2383  // avoid singularity at polar axis
2384 #if(COORDSINGFIX)
2385  if(fabs(th)<SINGSMALL){
2386  if(th>=0) th=SINGSMALL;
2387  if(th<0) th=-SINGSMALL;
2388  }
2389  if(fabs(M_PI-th)<SINGSMALL){
2390  if(th>=M_PI) th=M_PI+SINGSMALL;
2391  if(th<M_PI) th=M_PI-SINGSMALL;
2392  }
2393 #endif
2394 
2395  // set aux vars
2396  dxdxprim(X,V,dxdxp);
2397  // DLOOPA(j) dxdxp[j]=dxdxptrue[j][j]; // defcoord==LOGRSINTH assumes transformation is diagonal
2398  // dx_dxp_calc(X,r,th,dx_dxp);
2399  //dx_dxp_dxp_calc(X,r,th,dx_dxp_dxp);
2400 
2401 
2402 
2403  // GODMARK
2404  // need to set second derivative analytically
2405  for(i=0;i<NDIM;i++) for(j=0;j<NDIM;j++) for(k=0;k<NDIM;k++){
2406  dxdxp_dxp[i][j][k]=0.0;
2407  }
2408 
2409 
2410  // t1 = rf(X1,X2);
2411  t1 = r;
2412  // t2 = thf(X1,X2);
2413  t2 = th;
2414  t3 = cos(t2);
2415  t4 = a*t3;
2416  t7 = (-t1+t4)*(t1+t4);
2417  t8 = M*M;
2418  t10 = t1*t1;
2419  t11 = a*a;
2420  t12 = t3*t3;
2421  t14 = t10+t11*t12;
2422  t15 = t14*t14;
2423  t17 = 1/t15/t14;
2424  conn[0][0][0] = -2.0*t7*t8*t1*t17;
2425  t21 = t10*t10;
2426  // t22 = diff(rf(X1,X2),X1);
2427  t22 = dxdxp[RR][1];
2428  t23 = t21*t22;
2429  t24 = t10*t1;
2430  t25 = M*t24;
2431  t26 = t25*t22;
2432  t28 = t24*t3;
2433  t29 = sin(t2);
2434  // t30 = diff(thf(X1,X2),X1);
2435  t30 = dxdxp[TH][1];
2436  t31 = t29*t30;
2437  t34 = M*t1;
2438  t35 = t22*t12;
2439  t38 = t12*t12;
2440  t39 = t38*t22;
2441  t40 = t29*t1;
2442  t41 = t12*t3;
2443  t42 = t41*t30;
2444  conn[0][0][1] = -M*(-t23-2.0*t26+(2.0*t28*t31+2.0*t34*t35+(t39+2.0*t40*t42)*t11)*t11)*t17;
2445  // t52 = diff(rf(X1,X2),X2);
2446  t52 = dxdxp[RR][2];
2447  t53 = t21*t52;
2448  // t56 = diff(thf(X1,X2),X2);
2449  t56 = dxdxp[TH][2];
2450  t57 = t29*t56;
2451  t60 = t52*t12;
2452  t63 = t38*t52;
2453  t64 = t41*t1;
2454  conn[0][0][2] = -M*(-t53-2.0*t25*t52+(2.0*t28*t57+2.0*t34*t60+(t63+2.0*t64*t57)*t11)*t11)*t17;
2455  t74 = -1.0+t3;
2456  t75 = 1.0+t3;
2457  t78 = t7*t74*t75;
2458  conn[0][0][3] = -2.0*t78*a*t1*t8*t17;
2459  conn[0][1][0] = conn[0][0][1];
2460  t83 = t30*t30;
2461  t84 = t21*t10;
2462  t86 = t22*t22;
2463  t89 = t24*t22;
2464  t90 = t3*t29;
2465  t91 = t90*t30;
2466  t97 = t86*t12;
2467  t102 = t22*t29;
2468  t104 = t64*t102;
2469  conn[0][1][1] = -2.0*(t83*t84-t21*t86-t25*t86+(2.0*t89*t91+2.0*t83*t21*t12+
2470  t34*t97+(t83*t10*t38+t38*t86+2.0*t104*t30)*t11)*t11)*M*t17;
2471  t114 = t22*t52;
2472  t116 = t30*t56;
2473  t119 = t24*t52;
2474  t121 = t114*t12;
2475  t123 = t21*t12;
2476  t126 = t90*t56;
2477  t129 = t52*t29;
2478  t130 = t129*t30;
2479  t132 = t10*t38;
2480  conn[0][1][2] = -2.0*(-t25*t114+t116*t84-t23*t52+(t119*t91+t34*t121+2.0*
2481  t116*t123+t89*t126
2482  +(t39*t52+t64*t130+t116*t132+t104*t56)*t11)*t11)*M*t17;
2483  t148 = 2.0*t28*t30;
2484  t149 = t102*t12;
2485  t152 = t41*t24;
2486  t154 = 2.0*t152*t30;
2487  t156 = 2.0*t64*t30;
2488  t157 = t39*t29;
2489  t158 = t30*t1;
2490  t159 = t38*t3;
2491  t161 = 2.0*t158*t159;
2492  t169 = a*t29*t17;
2493  conn[0][1][3] = -(2.0*t25*t102+t23*t29+(-t148-2.0*t34*t149+t154+(-t156-t157+t161)*t11)*t11)*M*t169;
2494  conn[0][2][0] = conn[0][0][2];
2495  conn[0][2][1] = conn[0][1][2];
2496  t171 = t52*t52;
2497  t172 = t171*M;
2498  t175 = t56*t56;
2499  t177 = t1*t12;
2500  t185 = t52*t41;
2501  conn[0][2][2] = -2.0*(-t172*t24-t171*t21+t175*t84+(t172*t177+2.0*t119*t126+
2502  2.0*t175*t21*t12
2503  +(t171*t38+2.0*t185*t40*t56+t175*t10*t38)*t11)*t11)*M*t17;
2504  t203 = 2.0*t152*t56;
2505  t204 = t129*t12;
2506  t208 = 2.0*t28*t56;
2507  t209 = t63*t29;
2508  t210 = t56*t1;
2509  t212 = 2.0*t210*t159;
2510  t214 = 2.0*t64*t56;
2511  conn[0][2][3] = (-2.0*t25*t129-t53*t29+(-t203+2.0*t34*t204+t208+(t209-t212+t214)*t11)*t11)*M*t169;
2512  conn[0][3][0] = conn[0][0][3];
2513  conn[0][3][1] = conn[0][1][3];
2514  conn[0][3][2] = conn[0][2][3];
2515  t221 = t21*t1;
2516  t222 = t24*t12;
2517  t224 = t10*t12;
2518  t227 = t1*t38;
2519  t236 = (-t221+(-2.0*t222+(-t224+t10)*M+(-t227+(t38-t12)*M)*t11)*t11)*t74*t75;
2520  conn[0][3][3] = -2.0*t236*t34*t17;
2521  t240 = t56*M;
2522  t241 = t240*t24;
2523  t242 = 2.0*t241;
2524  t243 = t56*t21;
2525  t245 = 2.0*t240*t177;
2526  t246 = t56*t10;
2527  t247 = t246*t12;
2528  t248 = t1*t3;
2529  t250 = 2.0*t248*t129;
2530  t251 = t56*t12;
2531  t258 = t30*t52;
2532  t260 = 1/(-t22*t56+t258);
2533  t261 = t260*t17;
2534  conn[1][0][0] = -(-t242+t243+(t245-t247+t250+t246-t251*t11)*t11)*M*t261;
2535  t263 = M*t22;
2536  t264 = t56*t38;
2537  t271 = (-t242+(t245-t247+t250+t246+(-t251+t264)*t11)*t11)*t260*t17;
2538  conn[1][0][1] = -t263*t271;
2539  t273 = M*t52;
2540  conn[1][0][2] = -t273*t271;
2541  t275 = t24*t29;
2542  t276 = t240*t275;
2543  t278 = t119*t3;
2544  t280 = t243*t29;
2545  t281 = t40*t12;
2546  t283 = 2.0*t240*t281;
2547  t284 = t29*t12;
2548  t285 = t246*t284;
2549  t286 = t52*t3;
2550  t288 = 2.0*t286*t1;
2551  t289 = t57*t10;
2552  t297 = a*t29*t260*t17;
2553  conn[1][0][3] = (-2.0*t276+2.0*t278+t280+(t283-t285+t288+t289-t56*t11*t284)*t11)*M*t297;
2554  conn[1][1][0] = conn[1][0][1];
2555  // t299 = diff(diff(thf(X1,X2),X1),X1);
2556  t299 = dxdxp_dxp[TH][1][1];
2557  t300 = t52*t299;
2558  t303 = t258*t221*t22;
2559  t305 = t56*t84;
2560  // t306 = diff(diff(rf(X1,X2),X1),X1);
2561  t306 = dxdxp_dxp[RR][1][1] ;
2562  t308 = t83*t56;
2563  t309 = t21*t24;
2564  t310 = t308*t309;
2565  t313 = 2.0*t308*t84;
2566  t314 = t56*t24;
2567  t320 = t30*t22;
2568  t325 = t258*t89*t12;
2569  t327 = t308*t221;
2570  t328 = t52*t83;
2571  t329 = t90*t21;
2572  t330 = t328*t329;
2573  t333 = t306*t12;
2574  t336 = t221*t12;
2575  t338 = 2.0*t308*t336;
2576  t340 = 4.0*t308*t123;
2577  t342 = t248*t29;
2578  t344 = 2.0*t52*t86*t342;
2579  t346 = t56*t86;
2580  t358 = t306*t38;
2581  t361 = t1*t22;
2582  t363 = t258*t361*t38;
2583  t366 = 2.0*t308*t222;
2584  t367 = t24*t38;
2585  t368 = t308*t367;
2586  t370 = t41*t29*t10;
2587  t372 = 2.0*t328*t370;
2588  t375 = 2.0*t308*t132;
2589  t380 = t159*t29;
2590  t381 = t380*t56;
2591  t384 = t308*t227;
2592  t385 = t38*t12;
2593  t387 = t328*t380;
2594  conn[1][1][1] = -(-t300*t84-2.0*t303+t305*t306-t310+(-t243*t86+t313-2.0*t314*t86*M)*M
2595  +(-2.0*t320*t3*t280-4.0*t325-t327+t330-3.0*t300*t123+3.0*t243*t333
2596  -t338+(t340+t344-t246*t97+t346*t10+2.0*t210*t97*M)*M
2597  +(-4.0*t320*t41*t289-3.0*t300*t132+3.0*t246*t358-2.0*t363-t366-t368+t372
2598  +(-t346*t12+t375+2.0*t346*t38)*M
2599  +(-2.0*t320*t381-t384-t300*t385+t387+t56*t306*t385)*t11)*t11)*t11)*t260*t17;
2600  // t399 = diff(diff(thf(X1,X2),X1),X2);
2601  t399 = dxdxp_dxp[TH][1][2];
2602  t400 = t52*t399;
2603  // t402 = diff(diff(rf(X1,X2),X1),X2);
2604  t402 = dxdxp_dxp[RR][1][2];
2605  t404 = t30*t175;
2606  t405 = t404*t309;
2607  t406 = t171*t30;
2608  t408 = t56*t221;
2609  t409 = t114*t408;
2610  t411 = 2.0*t404*t84;
2611  t412 = t52*t56;
2612  t418 = t402*t12;
2613  t421 = t404*t221;
2614  t425 = t114*t314*t12;
2615  t428 = 2.0*t404*t336;
2616  t431 = t22*t175;
2617  t432 = t431*t329;
2618  t433 = t10*t22;
2619  t434 = t433*t12;
2620  t437 = 4.0*t404*t123;
2621  t440 = 2.0*t171*t22*t342;
2622  t442 = t210*M;
2623  t448 = 2.0*t370*t431;
2624  t451 = t114*t210*t38;
2625  t453 = 2.0*t404*t222;
2626  t454 = t402*t38;
2627  t459 = t404*t367;
2628  t462 = 2.0*t404*t132;
2629  t467 = t404*t227;
2630  t469 = t431*t380;
2631  conn[1][1][2] = (t400*t84-t305*t402+t405+t406*t221+t409+(-t411+t412*t23+2.0*t114*t241)*M
2632  +(-3.0*t243*t418+t421+3.0*t400*t123+2.0*t425+t428+2.0*t406*t222+
2633  t432+(t412*t434-t437-t440-t412*t433-2.0*t121*t442)*M
2634  +(t448+t406*t227+t451+t453-3.0*t246*t454+3.0*t400*t132+t459
2635  +(t114*t251-t462-2.0*t114*t264)*M+(t467+t400*t385+t469
2636  -t56*t402*t385)*t11)*t11)*t11)*t260*t17;
2637  t480 = t286*t21;
2638  t481 = t57*t221;
2639  t486 = 2.0*t185*t10;
2640  t487 = t57*t222;
2641  t488 = 2.0*t487;
2642  t491 = t57*t227;
2643  t492 = t52*t159;
2644  t498 = (-t491+t492+(-t57*t12+t57*t38)*M)*t11;
2645  t501 = t480-t481+(2.0*t278-2.0*t276)*M+(t486-t488+(-t285+t289+t288+t283)*M+t498)*t11;
2646  conn[1][1][3] = t501*t22*t297;
2647  conn[1][2][0] = conn[1][0][2];
2648  conn[1][2][1] = conn[1][1][2];
2649  t504 = t171*t56;
2650  // t507 = diff(diff(thf(X1,X2),X2),X2);
2651  t507 = dxdxp_dxp[TH][2][2];
2652  t508 = t52*t507;
2653  // t510 = diff(diff(rf(X1,X2),X2),X2);
2654  t510 = dxdxp_dxp[RR][2][2];
2655  t512 = t175*t56;
2656  t521 = t512*t221;
2657  t528 = t52*t175;
2658  t530 = t510*t12;
2659  t553 = t510*t38;
2660  t556 = t512*t24;
2661  conn[1][2][2] = -(-2.0*t504*t221-t508*t84+t305*t510-t512*t309
2662  +(2.0*t512*t84-t504*t21-2.0*t504*t25)*M
2663  +(-2.0*t521*t12-3.0*t508*t123-4.0*t504*t222-t521-t528*t329+3.0*t243*t530
2664  +(-t504*t224+t504*t10+2.0*t342*t171*t52+4.0*t512*t21*t12+2.0*t12*t171*t442)*M
2665  +(-3.0*t508*t132-2.0*t504*t227-2.0*t528*t370+3.0*t246*t553-2.0*t556*t12-t556*t38
2666  +(2.0*t512*t10*t38-t504*t12+2.0*t504*t38)*M
2667  +(-t528*t380-t512*t1*t38-t508*t385+t56*t510*t385)*t11)*t11)*t11)*t260*t17;
2668  conn[1][2][3] = t501*t52*t297;
2669  conn[1][3][0] = conn[1][0][3];
2670  conn[1][3][1] = conn[1][1][3];
2671  conn[1][3][2] = conn[1][2][3];
2672  t588 = t84*t29;
2673  t607 = t29*t38;
2674  conn[1][3][3] = -(-t56*t309*t29+t286*t84+2.0*t240*t588
2675  +(-t481-2.0*t408*t284+t480+2.0*t185*t21+(-4.0*t185*t24+t280+3.0*t243*t284+4.0*t278
2676  +(-2.0*t314*t29+2.0*t487)*M)*M
2677  +(t492*t10+t486-t314*t607-t488+(-2.0*t492*t1+t289+t288-2.0*t285+3.0*t246*t607
2678  +(-2.0*t491+2.0*t210*t284)*M)*M+t498)*t11)*t11)*t29*t261;
2679  t627 = t30*t21;
2680  t628 = t30*M;
2681  t630 = 2.0*t628*t24;
2682  t631 = t30*t10;
2683  t632 = t631*t12;
2684  t634 = 2.0*t628*t177;
2685  t636 = 2.0*t361*t90;
2686  t637 = t30*t12;
2687  conn[2][0][0] = -(-t627+t630+(t632-t634-t631-t636+t637*t11)*t11)*M*t261;
2688  t651 = (-t630+(t634+t636+t631-t632+(-t637+t30*t38)*t11)*t11)*t260*t17;
2689  conn[2][0][1] = t263*t651;
2690  conn[2][0][2] = t273*t651;
2691  t652 = t628*t275;
2692  t654 = t89*t3;
2693  t656 = t627*t29;
2694  t657 = t631*t284;
2695  t659 = 2.0*t628*t281;
2696  t661 = 2.0*t361*t3;
2697  t662 = t31*t10;
2698  conn[2][0][3] = (2.0*t652-2.0*t654-t656+(t657-t659-t661-t662+t30*t11*t284)*t11)*M*t297;
2699  conn[2][1][0] = conn[2][0][1];
2700  t670 = t30*t86;
2701  t673 = t30*t84;
2702  t675 = t22*t299;
2703  t677 = t83*t30;
2704  t686 = t677*t221;
2705  t689 = t83*t22;
2706  t712 = t677*t24;
2707  conn[2][1][1] = -(2.0*t670*t221-t673*t306+t675*t84+t677*t309+(-2.0*t677*t84+t627*t86+2.0*t670*t25)*M
2708  +(2.0*t686*t12+t689*t329-3.0*t627*t333+t686+4.0*t670*t222+3.0*t675*t123
2709  +(-2.0*t86*t22*t1*t90+t631*t97-t670*t10-4.0*t677*t21*t12-2.0*t637*t86*t34)*M
2710  +(2.0*t712*t12+t712*t38+2.0*t670*t227+3.0*t675*t132-3.0*t631*t358+2.0*t689*t370
2711  +(t670*t12-2.0*t677*t10*t38-2.0*t670*t38)*M
2712  +(t675*t385-t30*t306*t385+t689*t380+t677*t1*t38)*t11)*t11)*t11)*t260*t17;
2713  t748 = t22*t399;
2714  conn[2][1][2] = -(t303+t310-t673*t402+t748*t84+t346*t221
2715  +(-t313+t258*t23+2.0*t258*t26)*M
2716  +(t327+t330+2.0*t325-3.0*t627*t418+2.0*t346*t222+3.0*t748*t123+t338
2717  +(-t340-t258*t433-t344+t258*t434-2.0*t60*t30*t361*M)*M
2718  +(t346*t227+t372+t363+t366+t368-3.0*t631*t454+3.0*t748*t132
2719  +(t258*t35-t375-2.0*t258*t39)*M+(t387+t748*t385+t384
2720  -t30*t402*t385)*t11)*t11)*t11)*t260*t17;
2721  t793 = t31*t221;
2722  t794 = t3*t22;
2723  t795 = t794*t21;
2724  t799 = t31*t222;
2725  t800 = 2.0*t799;
2726  t801 = t22*t41;
2727  t803 = 2.0*t801*t10;
2728  t806 = t31*t227;
2729  t807 = t22*t159;
2730  t813 = (-t806+t807+(t31*t38-t31*t12)*M)*t11;
2731  t816 = -t793+t795+(2.0*t654-2.0*t652)*M+(-t800+t803+(t661-t657+t662+t659)*M+t813)*t11;
2732  conn[2][1][3] = -t816*t22*t297;
2733  conn[2][2][0] = conn[2][0][2];
2734  conn[2][2][1] = conn[2][1][2];
2735  t822 = t22*t507;
2736  t831 = t3*t56;
2737  t845 = t41*t56;
2738  conn[2][2][2] = -(-t673*t510+2.0*t409+t405+t822*t84+(t406*t21-t411+2.0*t406*t25)*M
2739  +(-3.0*t627*t530+t421-t432+2.0*t130*t831*t21+t428+3.0*t822*t123+4.0*t425
2740  +(t406*t224-t406*t10-t440-t437-2.0*t406*t177*M)*M
2741  +(4.0*t130*t845*t10+2.0*t451+3.0*t822*t132+t453+t459-t448-3.0*t631*t553
2742  +(t406*t12-t462-2.0*t406*t38)*M+(2.0*t258*t381+t822*t385-t30*t510*t385
2743  +t467-t469)*t11)*t11)*t11)*t260*t17;
2744  conn[2][2][3] = -t816*t52*t297;
2745  conn[2][3][0] = conn[2][0][3];
2746  conn[2][3][1] = conn[2][1][3];
2747  conn[2][3][2] = conn[2][2][3];
2748  t891 = t30*t24;
2749  conn[2][3][3] = (t794*t84+2.0*t628*t588-t30*t309*t29
2750  +(t795-t793-2.0*t30*t221*t284+2.0*t801*t21
2751  +(3.0*t627*t284+t656-4.0*t801*t24+4.0*t654+(-2.0*t891*t29+2.0*t799)*M)*M
2752  +(t807*t10+t803-t891*t607-t800+(3.0*t631*t607-2.0*t807*t1+t661-2.0*t657+t662
2753  +(2.0*t158*t284-2.0*t806)*M)*M+t813)*t11)*t11)*t29*t261;
2754  t917 = a*M;
2755  conn[3][0][0] = -t7*t917*t17;
2756  t920 = t102*t10;
2757  t924 = 1/t29;
2758  t916 = t924*t17;
2759  conn[3][0][1] = -t917*(-t920+t148+(t156+t149)*t11)*t916;
2760  t928 = t129*t10;
2761  conn[3][0][2] = -t917*(t208-t928+(t214+t204)*t11)*t916;
2762  conn[3][0][3] = -t78*M*t11*t17;
2763  conn[3][1][0] = conn[3][0][1];
2764  t940 = t83*t29;
2765  t946 = t86*t29;
2766  t968 = t916;
2767  conn[3][1][1] = -(t940*t221+2.0*t794*t627+(4.0*t794*t891-t946*t10)*M
2768  +(4.0*t801*t631+2.0*t940*t222+(4.0*t361*t42+t946*t12)*M+(t940*t227+2.0*t807*t30)*t11)
2769  *t11)*a*t968;
2770  t970 = t29*t221;
2771  t991 = t52*t1;
2772  conn[3][1][2] = -(t116*t970+t286*t627+t794*t243
2773  +(2.0*t286*t891-t102*t52*t10+2.0*t794*t314)*M
2774  +(2.0*t185*t631+2.0*t801*t246+2.0*t116*t275*t12+(2.0*t361*t845+2.0*t991*t42+t102*t60)*M
2775  +(t116*t40*t38+t492*t30+t807*t56)*t11)*t11)*a*t968;
2776  t1024 = t38*t41;
2777  conn[3][1][3] = (t3*t30*t84+t970*t22+(3.0*t42*t21+2.0*t275*t35
2778  +(t102*t224+t148-t154-t920)*M
2779  +(t40*t39+3.0*t159*t30*t10+(t156-t161+t149-t157)*M
2780  +t1024*t30*t11)*t11)*t11)*t924*t17;
2781  conn[3][2][0] = conn[3][0][2];
2782  conn[3][2][1] = conn[3][1][2];
2783  t1035 = t175*t29;
2784  t1037 = t171*t29;
2785  conn[3][2][2] = -(2.0*t286*t243+t1035*t221+(-t1037*t10+4.0*t286*t314)*M
2786  +(4.0*t185*t246+2.0*t1035*t222+(t1037*t12+4.0*t991*t845)*M
2787  +(t1035*t227+2.0*t492*t56)*t11)*t11)*a*t968;
2788  conn[3][2][3] = -(-t970*t52-t831*t84+(-2.0*t275*t60-3.0*t845*t21
2789  +(t203-t208-t129*t224+t928)*M
2790  +(-t40*t63-3.0*t159*t56*t10+(t209+t212-t204-t214)*M
2791  -t1024*t56*t11)*t11)*t11)*t924*t17;
2792  conn[3][3][0] = conn[3][0][3];
2793  conn[3][3][1] = conn[3][1][3];
2794  conn[3][3][2] = conn[3][2][3];
2795  conn[3][3][3] = -t236*a*t17;
2796 
2797  return;
2798 
2799 }
2800 
2801 #undef M
2802 
2803 
2804 
2805 
2806 
2807 
2810 void mks_source_conn(FTYPE *pr, struct of_geom *ptrgeom, struct of_state *q,FTYPE *dU)
2811 {
2812  int ii,jj,kk;
2813  int i=0, j=0, k=0, l=0;
2814  FTYPE r,th,X[NDIM],V[NDIM],sigma,dxdxptrue[NDIM][NDIM];
2815 #ifdef WIN32
2816  extern FTYPE cot(FTYPE arg);
2817 #endif
2818  extern FTYPE csc(FTYPE arg);
2819  FTYPE b[NDIM],u[NDIM],bsq,en,rho;
2820  FTYPE dxdxp[NDIM];
2821 
2822 
2823 
2824  if(MBH!=1.0){
2825  dualfprintf(fail_file,"mks_source_conn not setup for MBH!=1.0\n");
2826  myexit(12);
2827  }
2828 
2829 
2830  ii=ptrgeom->i;
2831  jj=ptrgeom->j;
2832  kk=ptrgeom->k;
2833 
2834 
2835  bsq = dot(q->bcon, q->bcov);
2836  u[TT]=q->ucon[TT];
2837  u[RR]=q->ucon[RR];
2838  u[TH]=q->ucon[TH];
2839  u[PH]=q->ucon[PH];
2840 
2841  b[TT]=q->bcon[TT];
2842  b[RR]=q->bcon[RR];
2843  b[TH]=q->bcon[TH];
2844  b[PH]=q->bcon[PH];
2845 
2846  rho=pr[RHO];
2847  en=pr[UU];
2848 
2849  coord(ptrgeom->i, ptrgeom->j, ptrgeom->k,ptrgeom->p, X);
2850  // get bl coordinates
2851  bl_coord(X,V);
2852  r=V[1];
2853  th=V[2];
2854 
2855 
2856  // this is not exactly right, since derivative of metric is derivative of absolute values, but shouldn't/doesn't seem to matter much
2857  // follows gcov_func()
2858  if(POSDEFMETRIC){
2859  if(th<0.0){ th=-th;}
2860  if(th>M_PI) { th=M_PI-th; }
2861  }
2862  else{
2863  }
2864  // avoid singularity at polar axis
2865 #if(COORDSINGFIX)
2866  if(fabs(th)<SINGSMALL){
2867  if(th>=0) th=SINGSMALL;
2868  if(th<0) th=-SINGSMALL;
2869  }
2870  if(fabs(M_PI-th)<SINGSMALL){
2871  if(th>=M_PI) th=M_PI+SINGSMALL;
2872  if(th<M_PI) th=M_PI-SINGSMALL;
2873  }
2874 #endif
2875  // set aux vars
2876  dxdxprim(X,V,dxdxptrue);
2877  DLOOPA(j) dxdxp[j]=dxdxptrue[j][j]; // defcoord==LOGRSINTH assumes transformation is diagonal
2878 
2879 
2880  sigma=r*r+a*a*cos(th)*cos(th);
2881 
2882 
2883 
2884  if((WHICHEOM==WITHNOGDET)&&(NOGDETU0==1)){
2885  // see grmhd-fullsource-simplify.nb
2886 
2887  dU[UU]+=pow(sigma,-1.)*(-1.*(-1. - 2.*dxdxp[1]*r*pow(sigma,-1.))*
2888  (b[RR]*(-1.*b[TT]*pow(a,2.)*pow(cos(th),2.) +
2889  r*(2.*b[TT] + 2.*b[RR]*dxdxp[1] - 1.*b[TT]*r -
2890  2.*b[PH]*a*pow(sin(th),2.))) +
2891  u[RR]*(bsq + en*gam + rho)*
2892  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
2893  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
2894  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.)))) +
2895  (b[TH]*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
2896  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
2897  2.*b[PH]*a*pow(sin(th),2.))) -
2898  1.*u[TH]*(bsq + en*gam + rho)*
2899  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
2900  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
2901  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.))))*
2902  (-1.*dxdxp[2]*cot(th) +
2903  4.*(-1.*M_PI*X[2] + th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) +
2904  dxdxp[2]*pow(a,2.)*pow(sigma,-1.)*sin(2.*th)));
2905 
2906  }
2907  // else nothing to add then
2908 
2909 
2910  if((WHICHEOM==WITHNOGDET)&&(NOGDETU1==1)){
2911 
2912  // see grmhd-ksksp-mhd.nb and the other source*simplify*.nb files
2913 
2914  dU[U1]+=0.0625*(16.*dxdxp[1]*r*
2915  (0.5*bsq + en*(-1. + gam) -
2916  1.*sigma*pow(b[TH],2.)*pow(dxdxp[2],2.) +
2917  (bsq + en*gam + rho)*sigma*pow(dxdxp[2],2.)*pow(u[TH],2.))*
2918  pow(sigma,-1.) + 2.*(-1. - 2.*dxdxp[1]*r*pow(sigma,-1.))*
2919  (4.*bsq + 8.*en*(-1. + gam) -
2920  1.*b[RR]*dxdxp[1]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
2921  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
2922  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
2923  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
2924  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.))*
2925  pow(sigma,-1.) + dxdxp[1]*u[RR]*(bsq + en*gam + rho)*
2926  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
2927  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
2928  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
2929  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
2930  4.*u[PH]*a*pow(r,2.))*pow(sigma,-1.)) -
2931  1.*dxdxp[1]*(4.*r - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
2932  (-1.*b[TT]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
2933  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
2934  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
2935  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
2936  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.)) +
2937  u[TT]*(bsq + en*gam + rho)*
2938  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
2939  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
2940  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
2941  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
2942  4.*u[PH]*a*pow(r,2.)))*pow(sigma,-4.)*
2943  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.)) +
2944  8.*a*pow(dxdxp[1],2.)*(-1.*b[RR]*
2945  (-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
2946  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
2947  cos(2.*th)*pow(a,2.)*
2948  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
2949  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
2950  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.)) +
2951  u[RR]*(bsq + en*gam + rho)*
2952  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
2953  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
2954  cos(2.*th)*pow(a,2.)*
2955  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
2956  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
2957  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.)))*pow(sigma,-4.)*
2958  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.) +
2959  8.*dxdxp[1]*a*(-1.*b[TT]*
2960  (-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
2961  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
2962  cos(2.*th)*pow(a,2.)*
2963  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
2964  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
2965  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.)) +
2966  u[TT]*(bsq + en*gam + rho)*
2967  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
2968  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
2969  cos(2.*th)*pow(a,2.)*
2970  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
2971  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
2972  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.)))*pow(sigma,-4.)*
2973  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.) +
2974  dxdxp[1]*a*(-1.*b[PH]*
2975  (16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r - 8.*b[PH]*a*r +
2976  4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*cos(2.*th) +
2977  4.*b[RR]*dxdxp[1]*pow(a,2.) - 1.*b[PH]*pow(a,3.) +
2978  b[PH]*cos(4.*th)*pow(a,3.) + 8.*b[RR]*dxdxp[1]*pow(r,2.) -
2979  4.*b[PH]*a*pow(r,2.)) +
2980  u[PH]*(bsq + en*gam + rho)*
2981  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
2982  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
2983  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
2984  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
2985  4.*u[PH]*a*pow(r,2.)))*pow(sigma,-4.)*
2986  (pow(a,2.)*(sigma - 2.*pow(r,2.)) +
2987  2.*r*(pow(1.,2.)*(-2.*sigma + 4.*pow(r,2.)) + pow(sigma,2.)) +
2988  cos(2.*th)*pow(a,2.)*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.)))*
2989  pow(sin(th),2.) + 8.*dxdxp[1]*pow(sigma,-3.)*
2990  (bsq + 2.*en*(-1. + gam) -
2991  1.*b[PH]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
2992  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
2993  cos(2.*th)*pow(a,2.)*
2994  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
2995  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
2996  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.))*pow(sigma,-1.)*
2997  pow(sin(th),2.) + u[PH]*(bsq + en*gam + rho)*
2998  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
2999  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3000  cos(2.*th)*pow(a,2.)*
3001  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
3002  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
3003  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.))*pow(sigma,-1.)*
3004  pow(sin(th),2.))*(r*pow(sigma,2.) +
3005  pow(a,2.)*(-1.*pow(r,2.) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)
3006  ) + 2.*pow(sigma,-3.)*(4.*bsq + 8.*en*(-1. + gam) -
3007  1.*b[RR]*dxdxp[1]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
3008  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
3009  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
3010  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
3011  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.))*
3012  pow(sigma,-1.) + dxdxp[1]*u[RR]*(bsq + en*gam + rho)*
3013  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
3014  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
3015  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
3016  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
3017  4.*u[PH]*a*pow(r,2.))*pow(sigma,-1.))*
3018  (-1.*dxdxp[1]*(2. + r)*pow(r,3.) + pow(sigma,3.) +
3019  dxdxp[1]*pow(a,2.)*(2.*r*pow(cos(th),2.) +
3020  pow(a,2.)*pow(cos(th),4.) +
3021  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.))) +
3022  16.*dxdxp[1]*a*pow(sigma,-4.)*
3023  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3024  (-1.*pow(r,2.) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)*
3025  (b[PH]*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
3026  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3027  2.*b[PH]*a*pow(sin(th),2.))) -
3028  1.*u[PH]*(bsq + en*gam + rho)*
3029  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3030  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3031  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.)))) -
3032  32.*pow(dxdxp[1],2.)*pow(sigma,-4.)*
3033  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3034  (r*(1. + r) + pow(a,2.)*pow(cos(th),2.))*
3035  (b[RR]*(-1.*b[TT]*pow(a,2.)*pow(cos(th),2.) +
3036  r*(2.*b[TT] + 2.*b[RR]*dxdxp[1] - 1.*b[TT]*r -
3037  2.*b[PH]*a*pow(sin(th),2.))) +
3038  u[RR]*(bsq + en*gam + rho)*
3039  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3040  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3041  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.)))) +
3042  16.*dxdxp[1]*pow(sigma,-3.)*
3043  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3044  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3045  (0.5*bsq + en*(-1. + gam) +
3046  b[TT]*pow(sigma,-1.)*
3047  (b[TT]*pow(a,2.)*pow(cos(th),2.) +
3048  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3049  2.*b[PH]*a*pow(sin(th),2.))) -
3050  1.*u[TT]*(bsq + en*gam + rho)*pow(sigma,-1.)*
3051  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3052  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3053  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.)))) -
3054  2.*dxdxp[1]*dxdxp[2]*cos(th)*pow(a,2.)*
3055  (-1.*b[TH]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
3056  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
3057  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
3058  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
3059  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.)) +
3060  u[TH]*(bsq + en*gam + rho)*
3061  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
3062  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
3063  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
3064  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
3065  4.*u[PH]*a*pow(r,2.)))*pow(sigma,-2.)*sin(th) -
3066  8.*dxdxp[1]*dxdxp[2]*a*cos(th)*
3067  (-1.*b[TH]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3068  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3069  cos(2.*th)*pow(a,2.)*
3070  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
3071  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
3072  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.)) +
3073  u[TH]*(bsq + en*gam + rho)*
3074  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
3075  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3076  cos(2.*th)*pow(a,2.)*
3077  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
3078  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
3079  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.)))*pow(sigma,-3.)*
3080  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*sin(th) +
3081  16.*dxdxp[1]*dxdxp[2]*r*
3082  (b[TH]*b[TT] - 1.*u[TH]*u[TT]*(bsq + en*gam + rho))*pow(a,2.)*
3083  pow(sigma,-2.)*sin(2.*th) +
3084  16.*dxdxp[2]*r*(b[RR]*b[TH] -
3085  1.*u[RR]*u[TH]*(bsq + en*gam + rho))*pow(dxdxp[1],2.)*pow(a,2.)*
3086  pow(sigma,-2.)*sin(2.*th) -
3087  16.*dxdxp[1]*dxdxp[2]*r*pow(a,2.)*pow(sigma,-3.)*
3088  (b[TH]*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
3089  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3090  2.*b[PH]*a*pow(sin(th),2.))) -
3091  1.*u[TH]*(bsq + en*gam + rho)*
3092  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3093  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3094  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.))))*sin(2.*th) +
3095  2.*dxdxp[1]*(-1.*b[TH]*
3096  (16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r - 8.*b[PH]*a*r +
3097  4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*cos(2.*th) +
3098  4.*b[RR]*dxdxp[1]*pow(a,2.) - 1.*b[PH]*pow(a,3.) +
3099  b[PH]*cos(4.*th)*pow(a,3.) + 8.*b[RR]*dxdxp[1]*pow(r,2.) -
3100  4.*b[PH]*a*pow(r,2.)) +
3101  u[TH]*(bsq + en*gam + rho)*
3102  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
3103  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
3104  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
3105  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
3106  4.*u[PH]*a*pow(r,2.)))*pow(sigma,-1.)*
3107  (-1.*dxdxp[2]*cot(th) +
3108  4.*(-1.*M_PI*X[2] + th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) +
3109  dxdxp[2]*pow(a,2.)*pow(sigma,-1.)*sin(2.*th)) -
3110  16.*dxdxp[1]*dxdxp[2]*a*
3111  (b[PH]*b[TH] - 1.*u[PH]*u[TH]*(bsq + en*gam + rho))*
3112  pow(sigma,-2.)*sin(th)*(r*(2. + r)*sigma*cos(th) +
3113  sigma*pow(a,2.)*pow(cos(th),3.) + r*pow(a,2.)*sin(th)*sin(2.*th)));
3114 
3115  }
3116  else{
3117 
3118  dU[U1]+=0.5*dxdxp[1]*pow(sigma,-5.)*
3119  (r*pow(sigma,4.)*(bsq - 2.*en + 2.*en*gam -
3120  2.*sigma*pow(dxdxp[2],2.)*pow(b[TH],2.) +
3121  (bsq + en*gam + rho)*pow(dxdxp[2],2.)*
3122  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.))*pow(u[TH],2.)) +
3123  a*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3124  (-1.*sigma*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH])*
3125  (r*(2. + 3.*r)*pow(a,2.) +
3126  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3127  2.*pow(r,4.)) - 2.*a*
3128  ((bsq + 2.*en*(-1. + gam))*r - 1.*dxdxp[1]*sigma*b[TT]*b[RR] +
3129  0.5*dxdxp[1]*(bsq + en*gam + rho)*u[TT]*u[RR]*
3130  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)))*
3131  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) -
3132  4.*a*r*sigma*(-0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*
3133  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) - 1.*pow(b[TT],2.) +
3134  (bsq + en*gam + rho)*pow(u[TT],2.)))*pow(sin(th),2.) +
3135  dxdxp[1]*a*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3136  (-4.*a*r*((bsq + 2.*en*(-1. + gam))*r - 1.*dxdxp[1]*sigma*b[TT]*b[RR] +
3137  dxdxp[1]*(bsq + en*gam + rho)*sigma*u[TT]*u[RR])*
3138  pow(dxdxp[1],-1.) +
3139  sigma*(r*(2. + 3.*r)*pow(a,2.) +
3140  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3141  2.*pow(r,4.))*(-1.*b[RR]*b[PH] + (bsq + en*gam + rho)*u[RR]*u[PH] +
3142  0.5*a*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-1.)*pow(sigma,-1.))
3143  - 1.*a*pow(dxdxp[1],-1.)*
3144  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3145  ((bsq + 2.*en*(-1. + gam))*((-2. + r)*r + pow(a,2.)) -
3146  2.*sigma*pow(dxdxp[1],2.)*pow(b[RR],2.) +
3147  (bsq + en*gam + rho)*pow(dxdxp[1],2.)*
3148  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.))*pow(u[RR],2.)))*
3149  pow(sin(th),2.) - 1.*sigma*
3150  (4.*r - 1.*pow(a,2.) + cos(2.*th)*pow(a,2.))*
3151  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3152  (((bsq + 2.*en*(-1. + gam))*r - 1.*dxdxp[1]*sigma*b[TT]*b[RR] +
3153  dxdxp[1]*(bsq + en*gam + rho)*sigma*u[TT]*u[RR])*pow(sigma,-1.)*
3154  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) +
3155  2.*r*(-0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*
3156  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) - 1.*pow(b[TT],2.) +
3157  (bsq + en*gam + rho)*pow(u[TT],2.)) -
3158  1.*a*(-1.*b[TT]*b[PH] + (bsq + en*gam + rho)*u[TT]*u[PH])*
3159  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)) +
3160  (4.*a*r*sigma*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH]) -
3161  1.*a*(a*(bsq + 2.*en*(-1. + gam)) -
3162  1.*dxdxp[1]*b[RR]*b[PH]*
3163  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)) +
3164  dxdxp[1]*(bsq + en*gam + rho)*u[RR]*u[PH]*
3165  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)))*
3166  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) +
3167  0.5*(r*(2. + 3.*r)*pow(a,2.) +
3168  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3169  2.*pow(r,4.))*((bsq + 2.*en*(-1. + gam))*pow(csc(th),2.) -
3170  2.*sigma*pow(b[PH],2.) + 2.*(bsq + en*gam + rho)*sigma*pow(u[PH],2.))
3171  )*pow(sin(th),2.)*(r*pow(sigma,2.) +
3172  pow(a,2.)*(-1.*pow(r,2.) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)
3173  ) + 2.*a*sigma*(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3174  (-1.*pow(r,2.) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)*
3175  (r*(a*(bsq + 2.*en*(-1. + gam)) - 2.*dxdxp[1]*sigma*b[RR]*b[PH] +
3176  2.*dxdxp[1]*(bsq + en*gam + rho)*sigma*u[RR]*u[PH])*pow(sigma,-1.)\
3177  - 1.*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH])*
3178  ((2. - 1.*r)*r - 1.*pow(a,2.)*pow(cos(th),2.)) -
3179  2.*a*r*(0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*
3180  pow(csc(th),2.) - 1.*pow(b[PH],2.) +
3181  (bsq + en*gam + rho)*pow(u[PH],2.))*pow(sin(th),2.)) +
3182  a*sigma*(pow(a,2.)*(sigma - 2.*pow(r,2.)) +
3183  2.*r*(pow(1.,2.)*(-2.*sigma + 4.*pow(r,2.)) + pow(sigma,2.)) +
3184  cos(2.*th)*pow(a,2.)*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.)))*
3185  pow(sin(th),2.)*(2.*r*(-1.*b[TT]*b[PH] +
3186  (bsq + en*gam + rho)*u[TT]*u[PH]) +
3187  dxdxp[1]*(-1.*b[RR]*b[PH] + (bsq + en*gam + rho)*u[RR]*u[PH] +
3188  0.5*a*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-1.)*pow(sigma,-1.))
3189  *(r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) -
3190  1.*a*(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3191  (0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*pow(csc(th),2.) -
3192  1.*pow(b[PH],2.) + (bsq + en*gam + rho)*pow(u[PH],2.))*
3193  pow(sin(th),2.)) + sigma*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3194  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3195  ((bsq + 2.*en*(-1. + gam))*sigma +
3196  2.*((-2. + r)*r + pow(a,2.)*pow(cos(th),2.))*pow(b[TT],2.) -
3197  1.*(bsq + en*gam + rho)*
3198  (2.*(-2. + r)*r + pow(a,2.) + cos(2.*th)*pow(a,2.))*pow(u[TT],2.) +
3199  4.*r*b[TT]*(-1.*dxdxp[1]*b[RR] + a*b[PH]*pow(sin(th),2.)) +
3200  4.*r*(bsq + en*gam + rho)*u[TT]*
3201  (dxdxp[1]*u[RR] - 1.*a*u[PH]*pow(sin(th),2.))) +
3202  2.*sigma*(2.*r*(-1.*b[TT]*b[RR] + (bsq + en*gam + rho)*u[TT]*u[RR] +
3203  (bsq + 2.*en*(-1. + gam))*r*pow(dxdxp[1],-1.)*pow(sigma,-1.)) +
3204  dxdxp[1]*(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3205  (0.5*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-2.)*
3206  ((-2. + r)*r + pow(a,2.))*pow(sigma,-1.) - 1.*pow(b[RR],2.) +
3207  (bsq + en*gam + rho)*pow(u[RR],2.)) -
3208  1.*a*(-1.*b[RR]*b[PH] + (bsq + en*gam + rho)*u[RR]*u[PH] +
3209  0.5*a*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-1.)*pow(sigma,-1.))
3210  *(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.))*
3211  (-1.*dxdxp[1]*(2. + r)*pow(r,3.) + pow(sigma,3.) +
3212  dxdxp[1]*pow(a,2.)*(2.*r*pow(cos(th),2.) +
3213  pow(a,2.)*pow(cos(th),4.) +
3214  (pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.))) +
3215  4.*dxdxp[1]*sigma*(pow(r,2.) - 1.*pow(a,2.)*pow(cos(th),2.))*
3216  (r*(1. + r) + pow(a,2.)*pow(cos(th),2.))*
3217  (b[TT]*b[RR]*((-2. + r)*r + pow(a,2.)*pow(cos(th),2.)) -
3218  2.*dxdxp[1]*r*pow(b[RR],2.) + 2.*a*r*b[RR]*b[PH]*pow(sin(th),2.) +
3219  0.5*(bsq + en*gam + rho)*u[RR]*
3220  (-1.*u[TT]*(2.*(-2. + r)*r + pow(a,2.) + cos(2.*th)*pow(a,2.)) +
3221  4.*r*(dxdxp[1]*u[RR] - 1.*a*u[PH]*pow(sin(th),2.)))) -
3222  1.*dxdxp[2]*a*cos(th)*pow(sigma,2.)*
3223  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3224  (4.*a*r*(b[TT]*b[TH] - 1.*(bsq + en*gam + rho)*u[TT]*u[TH]) -
3225  1.*(b[TH]*b[PH] - 1.*(bsq + en*gam + rho)*u[TH]*u[PH])*
3226  (r*(2. + 3.*r)*pow(a,2.) +
3227  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3228  2.*pow(r,4.)) + 2.*dxdxp[1]*a*
3229  (b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH])*
3230  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)))*sin(th) -
3231  2.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,3.)*
3232  (2.*r*(-1.*b[TT]*b[TH] + (bsq + en*gam + rho)*u[TT]*u[TH]) +
3233  dxdxp[1]*(-1.*b[RR]*b[TH] + (bsq + en*gam + rho)*u[RR]*u[TH])*
3234  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) -
3235  1.*a*(-1.*b[TH]*b[PH] + (bsq + en*gam + rho)*u[TH]*u[PH])*
3236  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.))*sin(th) +
3237  2.*dxdxp[2]*r*(b[TT]*b[TH] - 1.*(bsq + en*gam + rho)*u[TT]*u[TH])*
3238  pow(a,2.)*pow(sigma,3.)*sin(2.*th) +
3239  2.*dxdxp[1]*dxdxp[2]*r*
3240  (b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH])*pow(a,2.)*pow(sigma,3.)*
3241  sin(2.*th) - 2.*dxdxp[2]*r*pow(a,2.)*pow(sigma,2.)*
3242  (-2.*dxdxp[1]*r*(b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH]) -
3243  1.*(b[TT]*b[TH] - 1.*(bsq + en*gam + rho)*u[TT]*u[TH])*
3244  ((2. - 1.*r)*r - 1.*pow(a,2.)*pow(cos(th),2.)) +
3245  2.*a*r*(b[TH]*b[PH] - 1.*(bsq + en*gam + rho)*u[TH]*u[PH])*pow(sin(th),2.)
3246  )*sin(2.*th) - 2.*dxdxp[2]*a*
3247  (b[TH]*b[PH] - 1.*(bsq + en*gam + rho)*u[TH]*u[PH])*pow(sigma,3.)*sin(th)*
3248  (r*(2. + r)*sigma*cos(th) + sigma*pow(a,2.)*pow(cos(th),3.) +
3249  r*pow(a,2.)*sin(th)*sin(2.*th)));
3250  }
3251 
3252 
3253 
3254  if((WHICHEOM==WITHNOGDET)&&(NOGDETU2==1)){
3255 
3256  dU[U2]+=dxdxp[1]*r*(-1.*b[RR]*b[TH] +
3257  u[RR]*u[TH]*(bsq + en*gam + rho))*pow(dxdxp[2],2.) -
3258  1.*(-1.*b[RR]*b[TH] + u[RR]*u[TH]*(bsq + en*gam + rho))*
3259  (2.*dxdxp[1]*r + sigma)*pow(dxdxp[2],2.) -
3260  0.125*r*pow(dxdxp[2],2.)*((-2. + r)*r + pow(a,2.))*
3261  (-1.*b[TH]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
3262  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
3263  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
3264  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
3265  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.)) +
3266  u[TH]*(bsq + en*gam + rho)*
3267  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
3268  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
3269  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
3270  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
3271  4.*u[PH]*a*pow(r,2.)))*pow(sigma,-2.) -
3272  0.5*a*r*pow(dxdxp[2],2.)*(-1.*b[TH]*
3273  (-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3274  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3275  cos(2.*th)*pow(a,2.)*(-1.*b[RR]*dxdxp[1]*a +
3276  b[PH]*(-2. + r)*r + b[PH]*pow(a,2.)) -
3277  1.*b[RR]*dxdxp[1]*pow(a,3.) + b[PH]*pow(a,4.) +
3278  2.*b[PH]*pow(r,4.)) +
3279  u[TH]*(bsq + en*gam + rho)*
3280  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) + dxdxp[1]*u[RR]*r) +
3281  u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3282  cos(2.*th)*pow(a,2.)*(-1.*dxdxp[1]*u[RR]*a +
3283  u[PH]*(-2. + r)*r + u[PH]*pow(a,2.)) -
3284  1.*dxdxp[1]*u[RR]*pow(a,3.) + u[PH]*pow(a,4.) +
3285  2.*u[PH]*pow(r,4.)))*pow(sigma,-2.)*pow(sin(th),2.) -
3286  2.*pow(dxdxp[2],2.)*pow(r,2.)*pow(sigma,-2.)*
3287  (b[TH]*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
3288  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3289  2.*b[PH]*a*pow(sin(th),2.))) -
3290  1.*u[TH]*(bsq + en*gam + rho)*
3291  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3292  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3293  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.)))) +
3294  2.*dxdxp[2]*r*cos(th)*pow(a,3.)*pow(sigma,-3.)*
3295  (b[PH]*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
3296  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3297  2.*b[PH]*a*pow(sin(th),2.))) -
3298  1.*u[PH]*(bsq + en*gam + rho)*
3299  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3300  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3301  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.))))*pow(sin(th),3.) -
3302  1.*dxdxp[2]*a*r*cos(th)*(-1.*b[TT]*
3303  (-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3304  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3305  cos(2.*th)*pow(a,2.)*(-1.*b[RR]*dxdxp[1]*a +
3306  b[PH]*(-2. + r)*r + b[PH]*pow(a,2.)) -
3307  1.*b[RR]*dxdxp[1]*pow(a,3.) + b[PH]*pow(a,4.) +
3308  2.*b[PH]*pow(r,4.)) +
3309  u[TT]*(bsq + en*gam + rho)*
3310  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) + dxdxp[1]*u[RR]*r) +
3311  u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3312  cos(2.*th)*pow(a,2.)*(-1.*dxdxp[1]*u[RR]*a +
3313  u[PH]*(-2. + r)*r + u[PH]*pow(a,2.)) -
3314  1.*dxdxp[1]*u[RR]*pow(a,3.) + u[PH]*pow(a,4.) +
3315  2.*u[PH]*pow(r,4.)))*pow(sigma,-3.)*sin(th) -
3316  0.125*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*
3317  (4.*bsq + 8.*en*(-1. + gam) -
3318  1.*b[RR]*dxdxp[1]*(16.*b[TT]*r + 16.*b[RR]*dxdxp[1]*r -
3319  8.*b[PH]*a*r + 4.*a*(b[RR]*dxdxp[1]*a + b[PH]*r*(2. + r))*
3320  cos(2.*th) + 4.*b[RR]*dxdxp[1]*pow(a,2.) -
3321  1.*b[PH]*pow(a,3.) + b[PH]*cos(4.*th)*pow(a,3.) +
3322  8.*b[RR]*dxdxp[1]*pow(r,2.) - 4.*b[PH]*a*pow(r,2.))*
3323  pow(sigma,-1.) + dxdxp[1]*u[RR]*(bsq + en*gam + rho)*
3324  (16.*dxdxp[1]*u[RR]*r + 16.*u[TT]*r - 8.*u[PH]*a*r +
3325  4.*a*(dxdxp[1]*u[RR]*a + u[PH]*r*(2. + r))*cos(2.*th) +
3326  4.*dxdxp[1]*u[RR]*pow(a,2.) - 1.*u[PH]*pow(a,3.) +
3327  u[PH]*cos(4.*th)*pow(a,3.) + 8.*dxdxp[1]*u[RR]*pow(r,2.) -
3328  4.*u[PH]*a*pow(r,2.))*pow(sigma,-1.))*sin(th) -
3329  0.5*dxdxp[1]*dxdxp[2]*a*cos(th)*
3330  (-1.*b[RR]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3331  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3332  cos(2.*th)*pow(a,2.)*(-1.*b[RR]*dxdxp[1]*a +
3333  b[PH]*(-2. + r)*r + b[PH]*pow(a,2.)) -
3334  1.*b[RR]*dxdxp[1]*pow(a,3.) + b[PH]*pow(a,4.) +
3335  2.*b[PH]*pow(r,4.)) +
3336  u[RR]*(bsq + en*gam + rho)*
3337  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) + dxdxp[1]*u[RR]*r) +
3338  u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3339  cos(2.*th)*pow(a,2.)*(-1.*dxdxp[1]*u[RR]*a +
3340  u[PH]*(-2. + r)*r + u[PH]*pow(a,2.)) -
3341  1.*dxdxp[1]*u[RR]*pow(a,3.) + u[PH]*pow(a,4.) +
3342  2.*u[PH]*pow(r,4.)))*pow(sigma,-3.)*
3343  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*sin(th) +
3344  (0.5*bsq + en*(-1. + gam) - 1.*sigma*pow(b[TH],2.)*pow(dxdxp[2],2.) +
3345  (bsq + en*gam + rho)*sigma*pow(dxdxp[2],2.)*pow(u[TH],2.))*
3346  (4.*(M_PI*X[2] - 1.*th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) -
3347  1.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*sin(th)) +
3348  dxdxp[1]*dxdxp[2]*r*pow(a,2.)*pow(sigma,-3.)*
3349  (b[RR]*(-1.*b[TT]*pow(a,2.)*pow(cos(th),2.) +
3350  r*(2.*b[TT] + 2.*b[RR]*dxdxp[1] - 1.*b[TT]*r -
3351  2.*b[PH]*a*pow(sin(th),2.))) +
3352  u[RR]*(bsq + en*gam + rho)*
3353  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3354  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3355  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.))))*sin(2.*th) -
3356  1.*dxdxp[2]*r*pow(a,2.)*pow(sigma,-2.)*
3357  (0.5*bsq + en*(-1. + gam) +
3358  b[TT]*pow(sigma,-1.)*(b[TT]*pow(a,2.)*pow(cos(th),2.) +
3359  r*(-2.*b[TT] - 2.*b[RR]*dxdxp[1] + b[TT]*r +
3360  2.*b[PH]*a*pow(sin(th),2.))) -
3361  1.*u[TT]*(bsq + en*gam + rho)*pow(sigma,-1.)*
3362  (u[TT]*pow(a,2.)*pow(cos(th),2.) +
3363  r*(-2.*dxdxp[1]*u[RR] - 2.*u[TT] + dxdxp[1]*u[TT] +
3364  u[TT]*R0 + 2.*u[PH]*a*pow(sin(th),2.))))*sin(2.*th) +
3365  0.5*dxdxp[2]*(bsq + 2.*en*(-1. + gam) -
3366  1.*b[PH]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3367  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3368  cos(2.*th)*pow(a,2.)*(-1.*b[RR]*dxdxp[1]*a +
3369  b[PH]*(-2. + r)*r + b[PH]*pow(a,2.)) -
3370  1.*b[RR]*dxdxp[1]*pow(a,3.) + b[PH]*pow(a,4.) +
3371  2.*b[PH]*pow(r,4.))*pow(sigma,-1.)*pow(sin(th),2.) +
3372  u[PH]*(bsq + en*gam + rho)*
3373  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) + dxdxp[1]*u[RR]*r) +
3374  u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3375  cos(2.*th)*pow(a,2.)*(-1.*dxdxp[1]*u[RR]*a +
3376  u[PH]*(-2. + r)*r + u[PH]*pow(a,2.)) -
3377  1.*dxdxp[1]*u[RR]*pow(a,3.) + u[PH]*pow(a,4.) +
3378  2.*u[PH]*pow(r,4.))*pow(sigma,-1.)*pow(sin(th),2.))*
3379  (cot(th) + r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th)) +
3380  (0.5*bsq + en*(-1. + gam) - 1.*sigma*pow(b[TH],2.)*pow(dxdxp[2],2.) +
3381  (bsq + en*gam + rho)*sigma*pow(dxdxp[2],2.)*pow(u[TH],2.))*
3382  (-1.*dxdxp[2]*cot(th) + 4.*(-1.*M_PI*X[2] + th)*
3383  pow(dxdxp[2],-1.)*pow(M_PI,2.) +
3384  dxdxp[2]*pow(a,2.)*pow(sigma,-1.)*sin(2.*th));
3385  }
3386  else{
3387 
3388 
3389  dU[U2]+=0.5*(-2.*dxdxp[1]*r*
3390  (b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH])*pow(dxdxp[2],2.) -
3391  1.*a*r*pow(dxdxp[2],2.)*pow(sigma,-2.)*
3392  (4.*a*r*(b[TT]*b[TH] - 1.*(bsq + en*gam + rho)*u[TT]*u[TH]) -
3393  1.*(b[TH]*b[PH] - 1.*(bsq + en*gam + rho)*u[TH]*u[PH])*
3394  (r*(2. + 3.*r)*pow(a,2.) +
3395  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3396  2.*pow(r,4.)) + 2.*dxdxp[1]*a*
3397  (b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH])*
3398  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)))*pow(sin(th),2.) -
3399  4.*pow(dxdxp[2],2.)*pow(r,2.)*pow(sigma,-2.)*
3400  (-2.*dxdxp[1]*r*(b[RR]*b[TH] - 1.*(bsq + en*gam + rho)*u[RR]*u[TH]) -
3401  1.*(b[TT]*b[TH] - 1.*(bsq + en*gam + rho)*u[TT]*u[TH])*
3402  ((2. - 1.*r)*r - 1.*pow(a,2.)*pow(cos(th),2.)) +
3403  2.*a*r*(b[TH]*b[PH] - 1.*(bsq + en*gam + rho)*u[TH]*u[PH])*pow(sin(th),2.)
3404  ) - 2.*r*pow(dxdxp[2],2.)*((-2. + r)*r + pow(a,2.))*pow(sigma,-2.)*
3405  (2.*r*(-1.*b[TT]*b[TH] + (bsq + en*gam + rho)*u[TT]*u[TH]) +
3406  dxdxp[1]*(-1.*b[RR]*b[TH] + (bsq + en*gam + rho)*u[RR]*u[TH])*
3407  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) -
3408  1.*a*(-1.*b[TH]*b[PH] + (bsq + en*gam + rho)*u[TH]*u[PH])*
3409  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.)) +
3410  4.*dxdxp[2]*r*cos(th)*pow(a,3.)*pow(sigma,-3.)*
3411  (r*(a*(bsq + 2.*en*(-1. + gam)) - 2.*dxdxp[1]*sigma*b[RR]*b[PH] +
3412  2.*dxdxp[1]*(bsq + en*gam + rho)*sigma*u[RR]*u[PH])*pow(sigma,-1.)\
3413  - 1.*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH])*
3414  ((2. - 1.*r)*r - 1.*pow(a,2.)*pow(cos(th),2.)) -
3415  2.*a*r*(0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*
3416  pow(csc(th),2.) - 1.*pow(b[PH],2.) +
3417  (bsq + en*gam + rho)*pow(u[PH],2.))*pow(sin(th),2.))*pow(sin(th),3.)
3418  - 2.*dxdxp[2]*a*r*cos(th)*pow(sigma,-4.)*
3419  (-1.*sigma*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH])*
3420  (r*(2. + 3.*r)*pow(a,2.) +
3421  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3422  2.*pow(r,4.)) - 2.*a*
3423  ((bsq + 2.*en*(-1. + gam))*r - 1.*dxdxp[1]*sigma*b[TT]*b[RR] +
3424  0.5*dxdxp[1]*(bsq + en*gam + rho)*u[TT]*u[RR]*
3425  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)))*
3426  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) -
3427  4.*a*r*sigma*(-0.5*(bsq + 2.*en*(-1. + gam))*pow(sigma,-1.)*
3428  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) - 1.*pow(b[TT],2.) +
3429  (bsq + en*gam + rho)*pow(u[TT],2.)))*sin(th) -
3430  1.*dxdxp[1]*dxdxp[2]*a*cos(th)*pow(sigma,-4.)*
3431  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3432  (-4.*a*r*((bsq + 2.*en*(-1. + gam))*r - 1.*dxdxp[1]*sigma*b[TT]*b[RR] +
3433  dxdxp[1]*(bsq + en*gam + rho)*sigma*u[TT]*u[RR])*
3434  pow(dxdxp[1],-1.) +
3435  sigma*(r*(2. + 3.*r)*pow(a,2.) +
3436  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3437  2.*pow(r,4.))*(-1.*b[RR]*b[PH] + (bsq + en*gam + rho)*u[RR]*u[PH] +
3438  0.5*a*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-1.)*pow(sigma,-1.))
3439  - 1.*a*pow(dxdxp[1],-1.)*
3440  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3441  ((bsq + 2.*en*(-1. + gam))*((-2. + r)*r + pow(a,2.)) -
3442  2.*sigma*pow(dxdxp[1],2.)*pow(b[RR],2.) +
3443  (bsq + en*gam + rho)*pow(dxdxp[1],2.)*
3444  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.))*pow(u[RR],2.)))*
3445  sin(th) - 2.*dxdxp[1]*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-2.)*
3446  (2.*r*(-1.*b[TT]*b[RR] + (bsq + en*gam + rho)*u[TT]*u[RR] +
3447  (bsq + 2.*en*(-1. + gam))*r*pow(dxdxp[1],-1.)*pow(sigma,-1.)) +
3448  dxdxp[1]*(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*
3449  (0.5*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-2.)*
3450  ((-2. + r)*r + pow(a,2.))*pow(sigma,-1.) - 1.*pow(b[RR],2.) +
3451  (bsq + en*gam + rho)*pow(u[RR],2.)) -
3452  1.*a*(-1.*b[RR]*b[PH] + (bsq + en*gam + rho)*u[RR]*u[PH] +
3453  0.5*a*(bsq + 2.*en*(-1. + gam))*pow(dxdxp[1],-1.)*pow(sigma,-1.))
3454  *(r*(2. + r) + pow(a,2.)*pow(cos(th),2.))*pow(sin(th),2.))*sin(th)\
3455  + (bsq - 2.*en + 2.*en*gam - 2.*sigma*pow(dxdxp[2],2.)*pow(b[TH],2.) +
3456  (bsq + en*gam + rho)*pow(dxdxp[2],2.)*
3457  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.))*pow(u[TH],2.))*
3458  (4.*(M_PI*X[2] - 1.*th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) -
3459  1.*dxdxp[2]*cos(th)*pow(a,2.)*pow(sigma,-1.)*sin(th)) -
3460  1.*dxdxp[2]*r*pow(a,2.)*pow(sigma,-3.)*
3461  ((bsq + 2.*en*(-1. + gam))*sigma +
3462  2.*((-2. + r)*r + pow(a,2.)*pow(cos(th),2.))*pow(b[TT],2.) -
3463  1.*(bsq + en*gam + rho)*
3464  (2.*(-2. + r)*r + pow(a,2.) + cos(2.*th)*pow(a,2.))*pow(u[TT],2.) +
3465  4.*r*b[TT]*(-1.*dxdxp[1]*b[RR] + a*b[PH]*pow(sin(th),2.)) +
3466  4.*r*(bsq + en*gam + rho)*u[TT]*
3467  (dxdxp[1]*u[RR] - 1.*a*u[PH]*pow(sin(th),2.)))*sin(2.*th) -
3468  2.*dxdxp[1]*dxdxp[2]*r*pow(a,2.)*pow(sigma,-3.)*
3469  (b[TT]*b[RR]*((-2. + r)*r + pow(a,2.)*pow(cos(th),2.)) -
3470  2.*dxdxp[1]*r*pow(b[RR],2.) + 2.*a*r*b[RR]*b[PH]*pow(sin(th),2.) +
3471  0.5*(bsq + en*gam + rho)*u[RR]*
3472  (-1.*u[TT]*(2.*(-2. + r)*r + pow(a,2.) + cos(2.*th)*pow(a,2.)) +
3473  4.*r*(dxdxp[1]*u[RR] - 1.*a*u[PH]*pow(sin(th),2.))))*sin(2.*th) +
3474  dxdxp[2]*pow(sigma,-2.)*
3475  (4.*a*r*sigma*(b[TT]*b[PH] - 1.*(bsq + en*gam + rho)*u[TT]*u[PH]) -
3476  1.*a*(a*(bsq + 2.*en*(-1. + gam)) -
3477  1.*dxdxp[1]*b[RR]*b[PH]*
3478  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)) +
3479  dxdxp[1]*(bsq + en*gam + rho)*u[RR]*u[PH]*
3480  (pow(a,2.) + cos(2.*th)*pow(a,2.) + 2.*pow(r,2.)))*
3481  (r*(2. + r) + pow(a,2.)*pow(cos(th),2.)) +
3482  0.5*(r*(2. + 3.*r)*pow(a,2.) +
3483  cos(2.*th)*pow(a,2.)*((-2. + r)*r + pow(a,2.)) + pow(a,4.) +
3484  2.*pow(r,4.))*((bsq + 2.*en*(-1. + gam))*pow(csc(th),2.) -
3485  2.*sigma*pow(b[PH],2.) + 2.*(bsq + en*gam + rho)*sigma*pow(u[PH],2.))
3486  )*pow(sin(th),2.)*(cot(th) + r*pow(a,2.)*pow(sigma,-2.)*sin(2.*th)))\
3487  ;
3488  }
3489 
3490 
3491 
3492  if((WHICHEOM==WITHNOGDET)&&(NOGDETU3==1)){
3493 
3494  dU[U3]+=0.5*pow(sigma,-1.)*pow(sin(th),2.)*
3495  ((-1.*b[RR]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3496  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3497  cos(2.*th)*pow(a,2.)*
3498  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
3499  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
3500  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.)) +
3501  u[RR]*(bsq + en*gam + rho)*
3502  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
3503  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3504  cos(2.*th)*pow(a,2.)*
3505  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
3506  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
3507  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.)))*
3508  (-1. - 2.*dxdxp[1]*r*pow(sigma,-1.)) +
3509  (-1.*b[TH]*(-2.*a*r*(2.*b[TT] + b[RR]*dxdxp[1]*(2. + r)) +
3510  b[PH]*r*(2. + 3.*r)*pow(a,2.) +
3511  cos(2.*th)*pow(a,2.)*
3512  (-1.*b[RR]*dxdxp[1]*a + b[PH]*(-2. + r)*r +
3513  b[PH]*pow(a,2.)) - 1.*b[RR]*dxdxp[1]*pow(a,3.) +
3514  b[PH]*pow(a,4.) + 2.*b[PH]*pow(r,4.)) +
3515  u[TH]*(bsq + en*gam + rho)*
3516  (-2.*a*r*(2.*(dxdxp[1]*u[RR] + u[TT]) +
3517  dxdxp[1]*u[RR]*r) + u[PH]*r*(2. + 3.*r)*pow(a,2.) +
3518  cos(2.*th)*pow(a,2.)*
3519  (-1.*dxdxp[1]*u[RR]*a + u[PH]*(-2. + r)*r +
3520  u[PH]*pow(a,2.)) - 1.*dxdxp[1]*u[RR]*pow(a,3.) +
3521  u[PH]*pow(a,4.) + 2.*u[PH]*pow(r,4.)))*
3522  (-1.*dxdxp[2]*cot(th) +
3523  4.*(-1.*M_PI*X[2] + th)*pow(dxdxp[2],-1.)*pow(M_PI,2.) +
3524  dxdxp[2]*pow(a,2.)*pow(sigma,-1.)*sin(2.*th)));
3525 
3526  }
3527  // else nothing to add then
3528 
3529 
3530 }
3531 
3532 
3533 
3536 void mks_unitheta_idxvol_func(int i, int j, int k, FTYPE *idxvol)
3537 {
3538 
3539  /*
3540  int k, l;
3541  FTYPE r,th,ph;
3542  FTYPE r1[2],th1[2],r2[2],th2[2];
3543  FTYPE X0[NDIM],X1[2][NDIM],X2[2][NDIM];
3544  // FTYPE cot(FTYPE arg);
3545 
3546 
3547  coord(i, j, CENT, X0);
3548  coord(i, j, FACE1, X1[0]);
3549  coord(i+1, j, FACE1, X1[1]);
3550  coord(i, j, FACE2, X2[0]);
3551  coord(i, j+1, FACE2, X2[1]);
3552 
3553  // get bl coordinates
3554  bl_coord(X0,&r,&th);
3555  bl_coord(X1[0],&r1[0],&th1[0]);
3556  bl_coord(X1[1],&r1[1],&th1[1]);
3557  bl_coord(X2[0],&r2[0],&th2[0]);
3558  bl_coord(X2[1],&r2[1],&th2[1]);
3559 
3560  */
3561 
3562 
3563 
3564 
3565  /* comment out for now until adjust everything
3566 
3567 
3568  if(MBH!=1.0){
3569  dualfprintf(fail_file,"mks_unitheta_idxvold_func not setup for MBH!=1.0\n");
3570  myexit(13);
3571  }
3572 
3573 
3574  // this is not exactly right, since derivative of metric is derivative of absolute values, but shouldn't/doesn't seem to matter much
3575  // follows gcov_func()
3576  if(POSDEFMETRIC){
3577  if(th<0.0){ th=-th;}
3578  if(th>M_PI) { th=M_PI-th; }
3579  }
3580  else{
3581  }
3582  // avoid singularity at polar axis
3583  #if(COORDSINGFIX)
3584  if(fabs(th)<SINGSMALL){
3585  if(th>=0) th=SINGSMALL;
3586  if(th<0) th=-SINGSMALL;
3587  }
3588  if(fabs(M_PI-th)<SINGSMALL){
3589  if(th>=M_PI) th=M_PI+SINGSMALL;
3590  if(th<M_PI) th=M_PI-SINGSMALL;
3591  }
3592  #endif
3593  */
3594 
3595 #define IDXR(a,R0,r,th,rl,rh) ((pow(a,2) + 2.*pow(r,2) + pow(a,2)*cos(2.*th))/((rh - 1.*rl)*(2.*R0 + rh + rl) + (pow(a,2) + 2.*pow(R0,2))*(log(-1.*R0 + rh) - 1.*log(-1.*R0 + rl)) + pow(a,2)*cos(2.*th)*(log(-1.*R0 + rh) - 1.*log(-1.*R0 + rl))))
3596 
3597 #define IDXTH(a,R0,r,th,thl,thh) ((-3.*M_PI*(pow(a,2) + 2.*pow(r,2) + pow(a,2)*cos(2.*th))*sin(th))/(cos(thh)*(pow(a,2) + 6.*pow(r,2) + pow(a,2)*cos(2.*thh)) - 1.*cos(thl)*(pow(a,2) + 6.*pow(r,2) + pow(a,2)*cos(2.*thl))))
3598 
3599 
3600  //#define IDXTH(a,R0,r,th,thl,thh) ((3.*(pow(a,2) + 2.*pow(r,2) + pow(a,2)*cos(2.*th))*sin(th))/((-1.*cos(thh) + cos(thl))*(2.*(pow(a,2) + 3.*pow(r,2)) + pow(a,2)*(cos(2.*thh) + 2.*cos(thh)*cos(thl) + cos(2.*thl)))))
3601 
3602  // fullth integrated -- not used currently
3603 #define FIDXR(a,R0,r,th) ((pow(a,2) + 2.*pow(r,2) + pow(a,2)*cos(2.*th))/(4.*R0*(-1.*R0 + r) + pow(-1.*R0 + r,2) + (pow(a,2) + 2.*pow(R0,2) + pow(a,2)*cos(2.*th))*log(-1.*R0 + r)))
3604 
3605 #define FIDXTH(a,R0,r,th) ((-3.*pow(a,2)*sin(2.*th) - 6.*pow(r,2)*tan(th))/(pow(a,2) + 6.*pow(r,2) + pow(a,2)*cos(2.*th)))
3606 
3607 
3608  /*
3609 
3610  idxvol[TT]=1.0; // really 1/dt, but changes in time
3611  // comment out non-volume regularization
3612  // idxvol[RR]=1.0/dx[1];
3613  //idxvol[TH]=1.0/dx[2];
3614  idxvol[RR]=IDXR(a,R0,r,th,r1[0],r1[1]);
3615  idxvol[TH]=IDXTH(a,R0,r,th,th2[0],th2[1]);
3616  idxvol[PH]=1.0/dx[3];
3617 
3618  dualfprintf(fail_file,"%d %d %21.15g %21.15g\n",i,j,idxvol[RR]*dx[1],idxvol[TH]*dx[2]);
3619 
3620  */
3621 
3622 
3623 
3624 
3625 }
3626 
3627 
3631 void gdetvol_func(struct of_geom *ptrgeom, FTYPE *gdettrue, FTYPE *EOMFUNCNAME, FTYPE *gdetvol)
3632 {
3633  int i,j,k,loc;
3634  FTYPE Xc[NDIM],Vc[NDIM];
3635  FTYPE Xim[NDIM],Vim[NDIM],Xip1[NDIM],Vip1[NDIM];
3636  FTYPE Xjm[NDIM],Vjm[NDIM],Xjp1[NDIM],Vjp1[NDIM];
3637  FTYPE Xkm[NDIM],Vkm[NDIM],Xkp1[NDIM],Vkp1[NDIM];
3638  FTYPE dr,dh,dp;
3639  FTYPE drold,dhold,dpold;
3640  FTYPE newjacvol,oldjacvol;
3641  FTYPE dxdxpc[NDIM][NDIM];
3642  int locarray[4];
3643  int pl,pliter;
3644 
3645 
3646  i=ptrgeom->i;
3647  j=ptrgeom->j;
3648  k=ptrgeom->k;
3649  loc=ptrgeom->p;
3650 
3651 
3652  // loc==CENT
3653 
3654  if(loc==CENT){ // e.g. for conserved quantities at CENT
3655  locarray[0]=loc;
3656  locarray[1]=FACE1; // left-right
3657  locarray[2]=FACE2; // up-down
3658  locarray[3]=FACE3; // in-out
3659  }
3660  else if(loc==FACE1){ // e.g. for fluxes F1
3661  locarray[0]=loc;
3662  locarray[1]=CENT;
3663  locarray[2]=CORN3;
3664  locarray[3]=CORN2;
3665  }
3666  else if(loc==FACE2){ // e.g. for fluxes F2
3667  locarray[0]=loc;
3668  locarray[1]=CORN3;
3669  locarray[2]=CENT;
3670  locarray[3]=CORN1;
3671  }
3672  else if(loc==FACE3){ // e.g. for fluxes F3
3673  locarray[0]=loc;
3674  locarray[1]=CORN2;
3675  locarray[2]=CORN1;
3676  locarray[3]=CENT;
3677  }
3678  else{
3679  // maybe don't need yet :: only for fields -- not used yet.
3680  locarray[0]=loc;
3681  locarray[1]=FACE1; // left-right
3682  locarray[2]=FACE2; // up-down
3683  locarray[3]=FACE3; // in-out
3684  }
3685 
3686 
3687  // rest of code independent of loc
3688 
3689  coord_ijk(i, j, k, locarray[0], Xc);
3690  bl_coord_ijk(i, j, k, locarray[0], Vc);
3691  dxdxprim_ijk(i, j, k, locarray[0], dxdxpc);
3692 
3693 
3694 
3695 
3697  //
3698  // NEW JAC
3699  coord_ijk(i, j, k, locarray[1], Xim);
3700  bl_coord_ijk(i, j, k, locarray[1], Vim);
3701 
3702  coord_ijk(i+1, j, k, locarray[1], Xip1); // ok to have i+1 since coord and bl_coord don't depend upon memory locations
3703  bl_coord_ijk(i+1, j, k, locarray[1], Vip1);
3704 
3705  coord_ijk(i, j, k, locarray[2], Xjm);
3706  bl_coord_ijk(i, j, k, locarray[2], Vjm);
3707 
3708  coord_ijk(i, j+1, k, locarray[2], Xjp1); // ok to have j+1 since coord and bl_coord don't depend upon memory locations
3709  bl_coord_ijk(i, j+1, k, locarray[2], Vjp1);
3710 
3711  coord_ijk(i, j, k, locarray[3], Xkm);
3712  bl_coord_ijk(i, j, k, locarray[3], Vkm);
3713 
3714  coord_ijk(i, j, k+1, locarray[3], Xkp1); // ok to have k+1 since coord and bl_coord don't depend upon memory locations
3715  bl_coord_ijk(i, j, k+1, locarray[3], Vkp1);
3716 
3717  dr = THIRD*(pow(Vip1[1],3.0)-pow(Vim[1],3.0)); // really \delta(r^3/3) = r^2 dr
3718  // dh doesn't reduce to Pi, but 2.0 that is the correct answer for any N2 resolution
3719  dh = -(cos(Vjp1[2])-cos(Vjm[2])); // really \delta(-cos(h)) = sinh dh
3720  // below is inconsistent with rest of code when N2!=1
3721  // dh = (N2==1) ? M_PI : -(cos(Vjp1[2])-cos(Vjm[2])); // really \delta(-cos(h)) = sinh dh
3722  dp = (Vkp1[3]-Vkm[3]); // just d\phi
3723 
3724  // finite volume of cell
3725  newjacvol = dr*dh*dp/(dx[1]*dx[2]*dx[3]);
3726 
3727 
3729  //
3730  // OLD JAC
3731  drold=Vc[1]*Vc[1]*dxdxpc[1][1]*dx[1]; // r^2 dr
3732  if(totalsize[2]==1) dhold=2.0;
3733  else dhold=sin(Vc[2])*dxdxpc[2][2]*dx[2]; // sinh dh
3734  //dhold=sin(Vc[2])*dxdxpc[2][2]*dx[2]; // sinh dh (oldjacvol consistent with gdettrue)
3735  dpold=2.0*M_PI*dx[3];
3736  oldjacvol = drold*dhold*dpold/(dx[1]*dx[2]*dx[3]); // only true if diagonal mapping
3737 
3738  // suppose gdettrue already corrected if wanted to be corrected
3739  // if(totalsize[2]==1) (*gdettrue) /= (M_PI*0.5); // correct eomfunc and gdettrue
3740  // if(WHICHEOM==WITHGDET) PLOOP(pliter,pl) EOMFUNCMAC(pl)=(*gdettrue); // else up to user to make sure right
3741 
3743  //
3744  // use below, works best. Only central conserved quantity operated on (source terms too)
3745  if(loc==CENT){
3746  PLOOP(pliter,pl) EOMFUNCASSIGN(pl)=newjacvol;
3747  *gdetvol=(*gdettrue)=newjacvol;
3748  }
3749  else *gdetvol=(*gdettrue);
3750 
3751 
3753  //
3754  // SOME OTHER ATTEMPTS for 2 lines "best" above:
3755  //
3757 
3758  // horrible
3759  // if(loc==CENT) *gdetvol=newjacvol;
3760  // else *gdetvol=(*gdettrue);
3761 
3762 
3763  // central region undershoots
3764  //if(loc==CENT) *gdetvol=newjacvol;
3765  // else *gdetvol=(*gdettrue);
3766 
3767  // big overshoot, so probably don't want to use for fluxes
3768  // PLOOP(pliter,pl) *gdetvol=(*gdettrue)=EOMFUNCMAC(pl)=newjacvol;
3769 
3770  // with newjac in potential and gdettrue here
3771  // little bit more of a spike at center
3772  // *gdetvol=(*gdettrue);
3773 
3774  // with oldjac in potential and gdettrue here
3775  // same as above : little bit more of a spike at center
3776  // *gdetvol=(*gdettrue);
3777 
3778 
3779  // *gdetvol = newjacvol;
3780  // *gdetvol = dr*dh*dp/5.0;
3781  //*gdetvol = (*gdettrue); // disable gdetvol but keep memory
3782  // PLOOP(pliter,pl) EOMFUNCMAC(pl)=(*gdettrue)=*gdetvol=newjacvol;
3783 
3784 
3785  // *gdetvol=newjacvol;
3786  // if(loc==CENT){
3787  // PLOOP(pliter,pl) EOMFUNCMAC(pl)=(*gdettrue)=*gdetvol=newjacvol;
3788  // PLOOP(pliter,pl) EOMFUNCMAC(pl)=(*gdettrue)=*gdetvol=oldjacvol;
3789  // PLOOP(pliter,pl) EOMFUNCMAC(pl)=*gdetvol=(*gdettrue);
3790  // }
3791  // else{
3792  // PLOOP(pliter,pl) EOMFUNCMAC(pl)=*gdetvol=(*gdettrue)=oldjacvol;
3793  //PLOOP(pliter,pl) EOMFUNCMAC(pl)=*gdetvol=(*gdettrue)=newjacvol;
3794  // }
3795 
3796 
3797  // if(loc==CENT){
3798  // dualfprintf(fail_file,"i=%d loc=%d :: gdettrue=%15.7g oldjac=%15.7g newjac=%15.7g :: dr=%g drold=%g :: dh=%g dhold=%g :: dp=%g dpold=%g :: %g %g %g\n",i,loc,(*gdettrue),oldjacvol,newjacvol,dr,drold,dh,dhold,dp,dpold,dx[1],dx[2],dx[3]);
3799  // }
3800 
3801  //(*gdettrue) = *gdetvol; // replace gdettrue with this version
3802 
3803 
3804 }
3805 
3806 
3807 
3808 
3809 
3810 
3811 
3812 
3813 
3814 
3815 
3816 
3817 
3818 
3819 
3820 
3821 
Generated on Fri May 20 2016 15:52:34 for HARM by   doxygen 1.8.3.1