HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
transforms.c
Go to the documentation of this file.
1 
7 #include "decs.h"
8 
9 // this file includes all coordinate transformations and velocity transformations
10 // No user functions, unless new transformation of coordinates required
11 
12 
13 // NOTE: When dealing with multiple coordinate geom structures, must keep the PRIMECOORD one separate since repoint pointer to storage of PRIMECOORD geom structure, and using that pointer again afterwards for non-PRIMECOORD or other uses would overwrite PRIMECOORD geom structure values used during evolution and elsewhere when get_geometry() called.
14 
15 
17 int bl2met2metp2v(int whichvel, int whichcoord, FTYPE *pr, int ii, int jj, int kk)
18 {
19  int loc;
20 
21  loc=CENT;
22 
23  return(bl2met2metp2v_genloc(whichvel, whichcoord, pr, ii, jj, kk, loc));
24 
25 }
26 
27 
30 int bl2met2metp2v_genloc(int whichvel, int whichcoord, FTYPE *pr, int ii, int jj, int kk, int loc)
31 {
32  int k = 0;
33  FTYPE ucon[NDIM],uradcon[NDIM];
34  struct of_geom geomdontusebl;
35  struct of_geom *ptrgeombl=&geomdontusebl;
36  struct of_geom geomdontuse;
37  struct of_geom *ptrgeom=&geomdontuse;
38  FTYPE Bcon[NDIM];
39 
40 
41 
42  // whichvel==0 means supplied the 4-velocity
43  // whichvel==1 means supplied the 3-velocity
44  // whichvel==2 means supplied the relative 4-velocity
45 
46  // if whichcoord==PRIMECOORDS, then really use uses pr2ucon and ucon2pr, could probably optimize if wanted
47  // effectively this results in changing from one primitive velocity to another within PRIMECOORDS
48 
49 
50  // pr is in whichcoord coordinates
51  // get geometry (non-prime coords)
52  gset_genloc(0,whichcoord,ii,jj,kk,loc,ptrgeombl);
53 
54  // convert whichvel-pr in whichcoord coords to ucon in whichcoord coordinates
55  if (pr2ucon(whichvel,pr, ptrgeombl ,ucon) >= 1) FAILSTATEMENT("transforms.c:bl2met2metp2v_genloc()", "pr2ucon()", 1);
56  if(EOMRADTYPE!=EOMRADNONE){
57  if (pr2ucon(whichvel,&pr[URAD1-U1], ptrgeombl ,uradcon) >= 1) FAILSTATEMENT("transforms.c:bl2met2metp2v_genloc() for radiation", "pr2ucon()", 2);
58  }
59 
60  // convert field
61  Bcon[0]=0.0;
62  SLOOPA(k) Bcon[k]=pr[B1+k-1];
63 
64  // convert from whichcoord to MCOORD, the coordinates of evolution
65  if(whichcoord>=0){
66  coordtrans(whichcoord,MCOORD,ii,jj,kk,loc,ucon);
67  // transform MCOORD ucon from MCOORD non-prime to MCOORD prime coords
68  mettometp_genloc(ii,jj,kk,loc,ucon);
69 
70  if(EOMRADTYPE!=EOMRADNONE){
71  coordtrans(whichcoord,MCOORD,ii,jj,kk,loc,uradcon);
72  mettometp_genloc(ii,jj,kk,loc,uradcon);
73  }
74 
75  // field
76  coordtrans(whichcoord,MCOORD,ii,jj,kk,loc,Bcon);
77  mettometp_genloc(ii,jj,kk,loc,Bcon);
78 
79  }
80  // otherwise already in prime
81 
82  // get prime geometry
83  get_geometry(ii,jj,kk,loc,ptrgeom) ;
84 
85  // convert from MCOORD prime 4-vel to MCOORD prime WHICHVEL-vel(i.e. primitive velocity of evolution)
86  ucon2pr(WHICHVEL,ucon,ptrgeom,pr);
87 
88  if(EOMRADTYPE!=EOMRADNONE){
89  ucon2pr(WHICHVEL,uradcon,ptrgeom,&pr[URAD1-U1]);
90  }
91 
92  // convert field
93  SLOOPA(k) pr[B1+k-1]=Bcon[k];
94 
95  return(0);
96 }
97 
98 
99 
100 
102 int bl2met2metp2v_genloc_fieldonly(int whichvel, int whichcoord, FTYPE *pr, int ii, int jj, int kk, int loc)
103 {
104  int k = 0;
105  FTYPE ucon[NDIM],uradcon[NDIM];
106  struct of_geom geomdontusebl;
107  struct of_geom *ptrgeombl=&geomdontusebl;
108  struct of_geom geomdontuse;
109  struct of_geom *ptrgeom=&geomdontuse;
110  FTYPE Bcon[NDIM];
111 
112 
113 
114  // whichvel==0 means supplied the 4-velocity
115  // whichvel==1 means supplied the 3-velocity
116  // whichvel==2 means supplied the relative 4-velocity
117 
118  // if whichcoord==PRIMECOORDS, then really use uses pr2ucon and ucon2pr, could probably optimize if wanted
119  // effectively this results in changing from one primitive velocity to another within PRIMECOORDS
120 
121 
122  // pr is in whichcoord coordinates
123  // get geometry (non-prime coords)
124  gset_genloc(0,whichcoord,ii,jj,kk,loc,ptrgeombl);
125 
126  // convert field
127  Bcon[0]=0.0;
128  SLOOPA(k) Bcon[k]=pr[B1+k-1];
129 
130  // convert from whichcoord to MCOORD, the coordinates of evolution
131  if(whichcoord>=0){
132  // field
133  coordtrans(whichcoord,MCOORD,ii,jj,kk,loc,Bcon);
134  mettometp_genloc(ii,jj,kk,loc,Bcon);
135 
136  }
137  // otherwise already in prime
138 
139  // get prime geometry
140  get_geometry(ii,jj,kk,loc,ptrgeom) ;
141 
142  // convert field
143  SLOOPA(k) pr[B1+k-1]=Bcon[k];
144 
145  return(0);
146 }
147 
148 
149 
150 
151 
153 int ucov_whichcoord2primecoords(int whichcoord, int ii, int jj, int kk, int loc, FTYPE *ucov)
154 {
155  int k = 0;
156  FTYPE ucon[NDIM],uradcon[NDIM];
157  struct of_geom geomdontuse;
158  struct of_geom *ptrgeom=&geomdontuse;
159  struct of_geom geomprimedontuse;
160  struct of_geom *ptrgeomprime=&geomprimedontuse;
161  FTYPE idxdxp[NDIM][NDIM];
162 
163  // pr is in whichcoord coordinates
164  // get geometry (non-prime coords)
165  gset_genloc(0,whichcoord,ii,jj,kk,loc, ptrgeom);
166 
167  // first raise vector in whichcoord coordinates to MCOORD coordinates
168  raise_vec(ucov,ptrgeom,ucon);
169  coordtrans(whichcoord,MCOORD,ii,jj,kk,loc,ucon);
170 
171  // get idxdxp
172  idxdxprim_ijk(ii, jj, kk, loc, idxdxp);
173 
174  // convert to PRIMECOORDS
175  mettometp_simple(idxdxp, ucon);
176 
177  // get prime geometry
178  get_geometry(ii,jj,kk,loc,ptrgeomprime) ;
179 
180  // lower back
181  lower_vec(ucon,ptrgeomprime,ucov);
182 
183  return(0);
184 }
185 
186 
187 
189 int bl2met2metp2v_gen(int whichvel, int whichcoord, int newwhichvel, int newwhichcoord, FTYPE *pr, int ii, int jj, int kk)
190 {
191  int k = 0;
192  FTYPE ucon[NDIM],uradcon[NDIM];
193  struct of_geom geomdontusebl;
194  struct of_geom *ptrgeombl=&geomdontusebl;
195  struct of_geom geomdontuse;
196  struct of_geom *ptrgeom=&geomdontuse;
197  FTYPE Bcon[NDIM];
198 
199 
200  // whichvel==0 means supplied the 4-velocity
201  // whichvel==1 means supplied the 3-velocity
202  // whichvel==2 means supplied the relative 4-velocity
203 
204  // if whichcoord==PRIMECOORDS, then really use uses pr2ucon and ucon2pr, could probably optimize if wanted
205  // effectively this results in changing from one primitive velocity to another within PRIMECOORDS
206 
207 
208  // pr is in whichcoord coordinates
209  // get geometry (non-prime coords)
210  gset(0,whichcoord,ii,jj,kk,ptrgeombl);
211 
212  // convert whichvel-pr in whichcoord coords to ucon in whichcoord coordinates
213  if (pr2ucon(whichvel,pr, ptrgeombl ,ucon) >= 1) FAILSTATEMENT("transforms.c:bl2met2metp2v_gen()", "pr2ucon()", 1);
214  if(EOMRADTYPE!=EOMRADNONE){
215  if (pr2ucon(whichvel,&pr[URAD1-U1], ptrgeombl ,uradcon) >= 1) FAILSTATEMENT("transforms.c:bl2met2metp2v_gen() for radiation", "pr2ucon()", 2);
216  }
217 
218 
219  // field
220  Bcon[0]=0.0;
221  SLOOPA(k) Bcon[k]=pr[B1+k-1];
222 
223 
224  // convert from whichcoord to MCOORD, the coordinates of evolution
225  if(whichcoord>=0){
226  coordtrans(whichcoord,newwhichcoord,ii,jj,kk,CENT,ucon);
227  // transform MCOORD ucon from MCOORD non-prime to MCOORD prime coords
228  mettometp(ii,jj,kk,ucon);
229 
230  if(EOMRADTYPE!=EOMRADNONE){
231  coordtrans(whichcoord,newwhichcoord,ii,jj,kk,CENT,uradcon);
232  mettometp(ii,jj,kk,uradcon);
233  }
234 
235  coordtrans(whichcoord,newwhichcoord,ii,jj,kk,CENT,Bcon);
236  mettometp(ii,jj,kk,Bcon);
237 
238  }
239  // otherwise already in prime
240 
241  // get prime geometry
242  get_geometry(ii,jj,kk,CENT,ptrgeom) ;
243 
244  // convert from MCOORD prime 4-vel to MCOORD prime WHICHVEL-vel(i.e. primitive velocity of evolution)
245  ucon2pr(newwhichvel,ucon,ptrgeom,pr);
246  if(EOMRADTYPE!=EOMRADNONE){
247  ucon2pr(newwhichvel,uradcon,ptrgeom,&pr[URAD1-U1]);
248  }
249 
250  // convert field
251  SLOOPA(k) pr[B1+k-1]=Bcon[k];
252 
253  return(0);
254 }
255 
256 
258 int metp2met2bl(int whichvel, int whichcoord, FTYPE *pr, int ii, int jj, int kk)
259 {
260  int loc;
261 
262  loc=CENT;
263 
264  return(metp2met2bl_genloc(whichvel, whichcoord, pr, ii, jj, kk, loc));
265 
266 }
267 
269 int metp2met2bl_genloc(int whichvel, int whichcoord, FTYPE *pr, int ii, int jj, int kk, int pos)
270 {
271  int k;
272  FTYPE ucon[NDIM],uradcon[NDIM];
273  struct of_geom geomdontuse;
274  struct of_geom *ptrgeom=&geomdontuse;
275  struct of_geom geomdontusebl;
276  struct of_geom *ptrgeombl=&geomdontusebl;
277  FTYPE Bcon[NDIM];
278 
279 
280 
281  // which=WHICHVEL
282  // which==0 means supplied the 4-velocity
283  // which==1 means supplied the 3-velocity
284  // which==2 means supplied the relative 4-velocity
285 
286  // if whichcood==PRIMECOORDS, then just pr2ucon and ucon2pr
287  // effectively this results in changing from one primitive velocity to another within PRIMECOORDS
288 
289 
290  // get prime MCOORD geometry
291  get_geometry(ii,jj,kk,pos,ptrgeom) ;
292  // transform prime MCOORD primitive to prim MCOORD 4-vel
293 
294  // if (pr2ucon(WHICHVEL,pr, ptrgeom ,ucon) >= 1) FAILSTATEMENT("transforms.c:metp2met2bl()", "pr2ucon()", 1);
295  MYFUN(pr2ucon(WHICHVEL,pr, ptrgeom ,ucon),"transforms.c:pr2ucon()","metp2met2bl",0);
296  if(EOMRADTYPE!=EOMRADNONE){
297  MYFUN(pr2ucon(WHICHVEL,&pr[URAD1-U1], ptrgeom ,uradcon),"transforms.c:pr2ucon() for radiation","metp2met2bl",1);
298  }
299 
300  // field
301  Bcon[0]=0.0;
302  SLOOPA(k) Bcon[k]=pr[B1+k-1];
303 
304 
305 
306  if(whichcoord>=0){
307  // transform from prime MCOORD 4-vel to non-prime MCOORD 4-vel
308  metptomet_genloc(ii,jj,kk,pos,ucon);
309  // transform from non-prime MCOORD to non-prime whichcoord
310  coordtrans(MCOORD,whichcoord,ii,jj,kk,pos,ucon);
311 
312  if(EOMRADTYPE!=EOMRADNONE){
313  metptomet_genloc(ii,jj,kk,pos,uradcon);
314  coordtrans(MCOORD,whichcoord,ii,jj,kk,pos,uradcon);
315  }
316 
317  // convert field
318  metptomet_genloc(ii,jj,kk,pos,Bcon);
319  coordtrans(MCOORD,whichcoord,ii,jj,kk,pos,Bcon);
320 
321  }
322  // else already in prime
323 
324 
325  // transform from non-prime whichcoord 4-vel to non-prime whichcoord whichvel-velocity
326  // can't use same ptrgeom here, since ptrgeom would then the PRIMECOORD version would be overwritten by "whichcoord" version
327  gset_genloc(0,whichcoord,ii,jj,kk,pos,ptrgeombl);
328 
329  ucon2pr(whichvel,ucon,ptrgeombl,pr);
330 
331  if(EOMRADTYPE!=EOMRADNONE){
332  ucon2pr(whichvel,uradcon,ptrgeombl,&pr[URAD1-U1]);
333  }
334 
335  // convert field
336  SLOOPA(k) pr[B1+k-1]=Bcon[k];
337 
338  return(0);
339 }
340 
341 
342 
343 
344 
346 int coordtrans(int whichcoordin, int whichcoordout, int ii, int jj, int kk, int loc, FTYPE *ucon)
347 {
348  // GODMARK: need transformation from BL to KS_JP1 for EP3!=0 or THETAROT!=0
349  if(whichcoordin==whichcoordout){// then no transformation
350  return(0);
351  }
352  else if((whichcoordin==BLCOORDS)&&(whichcoordout==KSCOORDS||whichcoordout==KS_JP1_COORDS)){
353  bltoks(ii,jj,kk ,loc,ucon);
354  }
355  else if((whichcoordin==KSCOORDS || whichcoordin==KS_JP1_COORDS)&&(whichcoordout==BLCOORDS)){
356  kstobl(ii,jj,kk ,loc,ucon);
357  }
358  else{
359  dualfprintf(fail_file,"HARDFAIL: No such transformation: %d -> %d\n",whichcoordin,whichcoordout);
360  myexit(1);
361  }
362 
363  return(0);
364 
365 }
366 
367 // transformation matrix
368 void bltoks_trans(int ii, int jj, int kk, int loc, FTYPE (*bl2ks)[NDIM])
369 {
370  FTYPE V[NDIM], r, th;
371  int j,k;
372 
373 
374  bl_coord_ijk(ii,jj,kk,loc,V) ;
375  r=V[1]; th=V[2];
376 
377  // don't rotate, because assume bltoks() only called to make local transformation of u^i from BL to KS with same alignment, and notice no angle factors in transformation -- so nothing to transform anyways!
378  // FTYPE Vmetric[NDIM];
379  // rotate_V(BLCOORD,V,Vmetric);
380  // r=Vmetric[1]; th=Vmetric[2];
381 
382 
383  // bl2ks for contravariant components
384 #define bl2ks_trans00 (1)
385 #define bl2ks_trans01 (2.*r/(r*r - 2.*r + a*a))
386 #define bl2ks_trans02 (0)
387 #define bl2ks_trans03 (0)
388 #define bl2ks_trans10 (0)
389 #define bl2ks_trans11 (1)
390 #define bl2ks_trans12 (0)
391 #define bl2ks_trans13 (0)
392 #define bl2ks_trans20 (0)
393 #define bl2ks_trans21 (0)
394 #define bl2ks_trans22 (1)
395 #define bl2ks_trans23 (0)
396 #define bl2ks_trans30 (0)
397 #define bl2ks_trans31 (a/(r*r - 2.*r + a*a))
398 #define bl2ks_trans32 (0)
399 #define bl2ks_trans33 (1)
400 
401  /* make transform matrix */
402  // order for trans is [ourmetric][bl]
403  // DLOOP(j,k) bl2ks[j][k] = 0. ;
404  // DLOOPA(j) bl2ks[j][j] = 1. ;
405  bl2ks[0][0] = bl2ks_trans00;
406  bl2ks[0][1] = bl2ks_trans01;
407  bl2ks[0][2] = bl2ks_trans02;
408  bl2ks[0][3] = bl2ks_trans03;
409  bl2ks[1][0] = bl2ks_trans10;
410  bl2ks[1][1] = bl2ks_trans11;
411  bl2ks[1][2] = bl2ks_trans12;
412  bl2ks[1][3] = bl2ks_trans13;
413  bl2ks[2][0] = bl2ks_trans20;
414  bl2ks[2][1] = bl2ks_trans21;
415  bl2ks[2][2] = bl2ks_trans22;
416  bl2ks[2][3] = bl2ks_trans23;
417  bl2ks[3][0] = bl2ks_trans30;
418  bl2ks[3][1] = bl2ks_trans31;
419  bl2ks[3][2] = bl2ks_trans32;
420  bl2ks[3][3] = bl2ks_trans33;
421 
422  /* done! */
423 }
424 
425 
426 /* transforms u^i to our ks from boyer-lindquist */
427 void bltoks(int ii, int jj, int kk, int loc, FTYPE*ucon)
428 {
429  FTYPE tmp[NDIM];
430  FTYPE bl2ks[NDIM][NDIM];
431  FTYPE V[NDIM], r, th;
432  int j,k;
433 
434 
435  bltoks_trans(ii, jj, kk, loc, bl2ks);
436 
437 
438  /* transform ucon; solve for v */
439  // this is u^j[ks] = T^j[ks]_k[bl] u^k[bl]
440  DLOOPA(j) tmp[j] = 0.;
441  DLOOP(j,k) tmp[j] += bl2ks[j][k] * ucon[k];
442  DLOOPA(j) ucon[j] = tmp[j];
443 
444  /* done! */
445 }
446 
447 /* transforms u^i to our ks from boyer-lindquist */
448 void bltoks_ucov(int ii, int jj, int kk, int loc, FTYPE *ucov)
449 {
450  FTYPE tmp[NDIM];
451  FTYPE ks2bl[NDIM][NDIM];
452  FTYPE V[NDIM], r, th;
453  int j,k;
454 
455 
456  kstobl_trans(ii, jj, kk, loc, ks2bl);
457 
458 
459  /* transform ucon; solve for v */
460  // this is u_j[ks] = T_j[ks]^k[bl] u_k[bl] = T^k[bl]_j[ks] u_k[bl] = ks2bl^k[bl]_j[ks] u_k[bl]
461  DLOOPA(j) tmp[j] = 0.;
462  DLOOP(j,k) tmp[j] += ks2bl[k][j] * ucov[k];
463  DLOOPA(j) ucov[j] = tmp[j];
464 
465  /* done! */
466 }
467 
468 
469 // transformation matrix from ks to bl (inverse *and* transpose of bl2ks -- i.e. just Inverse[] in mathematica)
470 void kstobl_trans(int ii, int jj, int kk, int loc, FTYPE (*ks2bl)[NDIM])
471 {
472  FTYPE V[NDIM], r, th;
473  int j,k;
474 
475  bl_coord_ijk(ii,jj,kk,loc,V) ;
476  r=V[1]; th=V[2];
477 
478  // don't rotate, because assume kstobl() only called to make local transformation of u^i from BL to KS with same alignment, and notice no angle factors in transformation -- so nothing to transform anyways!
479  // FTYPE Vmetric[NDIM];
480  // rotate_V(KSCOORD,V,Vmetric);
481  // r=Vmetric[1]; th=Vmetric[2];
482 
483 
484  // just inverse (no transpose) of above
485 #define ks2bl_trans00 (1)
486 #define ks2bl_trans01 (-2.*r/(r*r - 2.*r + a*a))
487 #define ks2bl_trans02 (0)
488 #define ks2bl_trans03 (0)
489 #define ks2bl_trans10 (0)
490 #define ks2bl_trans11 (1)
491 #define ks2bl_trans12 (0)
492 #define ks2bl_trans13 (0)
493 #define ks2bl_trans20 (0)
494 #define ks2bl_trans21 (0)
495 #define ks2bl_trans22 (1)
496 #define ks2bl_trans23 (0)
497 #define ks2bl_trans30 (0)
498 #define ks2bl_trans31 (-a/(r*r - 2.*r + a*a))
499 #define ks2bl_trans32 (0)
500 #define ks2bl_trans33 (1)
501 
502 
503 
504  /* make transform matrix */
505  // order for trans is [bl][ourmetric] (i.e. inverse transpose of bl2ks)
506  // DLOOP(j,k) ks2bl[j][k] = 0. ;
507  // DLOOPA(j) ks2bl[j][j] = 1. ;
508  ks2bl[0][0] = ks2bl_trans00;
509  ks2bl[0][1] = ks2bl_trans01;
510  ks2bl[0][2] = ks2bl_trans02;
511  ks2bl[0][3] = ks2bl_trans03;
512  ks2bl[1][0] = ks2bl_trans10;
513  ks2bl[1][1] = ks2bl_trans11;
514  ks2bl[1][2] = ks2bl_trans12;
515  ks2bl[1][3] = ks2bl_trans13;
516  ks2bl[2][0] = ks2bl_trans20;
517  ks2bl[2][1] = ks2bl_trans21;
518  ks2bl[2][2] = ks2bl_trans22;
519  ks2bl[2][3] = ks2bl_trans23;
520  ks2bl[3][0] = ks2bl_trans30;
521  ks2bl[3][1] = ks2bl_trans31;
522  ks2bl[3][2] = ks2bl_trans32;
523  ks2bl[3][3] = ks2bl_trans33;
524 
525  /* done! */
526 }
527 
528 /* transforms u^i to our ks from boyer-lindquist */
529 void kstobl(int ii, int jj, int kk, int loc, FTYPE*ucon)
530 {
531  FTYPE tmp[NDIM];
532  FTYPE ks2bl[NDIM][NDIM];
533  FTYPE V[NDIM], r, th;
534  int j,k;
535 
536 
537  kstobl_trans(ii, jj, kk, loc, ks2bl);
538 
539  /* transform ucon; solve for v */
540  // this is u^j[bl] = T^j[bl]_k[ks] u^k[ks]
541  DLOOPA(j) tmp[j] = 0.;
542  DLOOP(j,k) tmp[j] += ks2bl[j][k] * ucon[k];
543  DLOOPA(j) ucon[j] = tmp[j];
544 
545  /* done! */
546 }
547 
548 
549 
554 void transV2Vmetric(int whichcoord, int ii, int jj, int kk, int loc, FTYPE ROTANGLE, FTYPE *X, FTYPE *V, FTYPE *Xmetric, FTYPE *Vmetric, FTYPE*gcov, FTYPE *gcovpert)
555 {
556 
557  if(ROTANGLE!=0.0 && ALLOWMETRICROT==1 && ISSPCMCOORD(whichcoord)){
558 
559 #if(0)// if input X,V,Xmetric,Vmetric, no longer need to duplicate doing below.
560  FTYPE V[NDIM];
561  bl_coord_ijk(ii,jj,kk,loc,V) ;
562 
563  // V[X] is not same as Vmetric (used in set_gcov), so get it.
564  FTYPE Vmetric[NDIM];
565  rotate_VtoVmetric(MCOORD,ROTANGLE,V,Vmetric);
566 #endif
567 
568  FTYPE told, r, h, ph;
569  told=Vmetric[0]; r=Vmetric[1]; h=Vmetric[2]; ph=Vmetric[3];
570  // now have Vmetric and can put into trans below
571 
572  /* make transform matrix */
573  // order for trans is [ourmetric][bl]
574  // DLOOP(j,k) trans[j][k] = 0. ;
575  // DLOOPA(j) trans[j][j] = 1. ;
576 
577  FTYPE b0=ROTANGLE;
578 
579 
580  FTYPE trans[NDIM][NDIM];
581  trans[0][0]=1.;
582  trans[0][1]=0.;
583  trans[0][2]=0.;
584  trans[0][3]=0.;
585  trans[1][0]=0.;
586  trans[1][1]=1.;
587  trans[1][2]=0.;
588  trans[1][3]=0.;
589  trans[2][0]=0.;
590  trans[2][1]=0.;
591  trans[2][2]=pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),-0.5)*(-1.*cos(h)*cos(ph)*sin(b0) + cos(b0)*sin(h));
592  trans[2][3]=pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),-0.5)*sin(b0)*sin(h)*sin(ph);
593  trans[3][0]=0.;
594  trans[3][1]=0.;
595  trans[3][2]=-1.*pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),-1.)*sin(b0)*sin(ph);
596  trans[3][3]=pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),-1.)*sin(h)*(-1.*cos(h)*cos(ph)*sin(b0) + cos(b0)*sin(h));
597 
598 
599  // now perform transformation
600  transgcovgcovpertself(gcov,gcovpert,trans);
601  }
602  else{
603  // then no change, and given function format gcov is just not changed.
604  }
605 
606  /* done! */
607 }
608 
609 
610 
611 
612 
616 void transVmetrictoV(int whichcoord, int ii, int jj, int kk, int loc, FTYPE ROTANGLE, FTYPE *X, FTYPE *V, FTYPE *Xmetric, FTYPE *Vmetric, FTYPE*gcov, FTYPE *gcovpert)
617 {
618 
619  if(ROTANGLE!=0.0 && ALLOWMETRICROT==1 && ISSPCMCOORD(whichcoord)){
620 
621 #if(0)
622  // gets V[X]
623  FTYPE V[NDIM];
624  bl_coord_ijk(ii,jj,kk,loc,V) ;
625 
626  // V[X] is not same as Vmetric (used in set_gcov)
627  // trans needs Vmetric to keep expression simple, so get it.
628  FTYPE Vmetric[NDIM];
629  rotate_VtoVmetric(MCOORD,ROTANGLE,V,Vmetric);
630  // now have Vmetric and can put into trans below
631 #endif
632 
633  // get transformation from Vmetric to V (i.e. rotation from BH spin along z-axis to rotated in -x direction by ROTANGLE)
634  FTYPE trans[NDIM][NDIM];
635  transVmetrictoV_trans(ROTANGLE, Vmetric, trans);
636 
637  // now perform transformation
638  transgcovgcovpertself(gcov,gcovpert,trans);
639  }
640  else{
641  // then no change, and given function format gcov is just not changed.
642  }
643 
644 
645  /* done! */
646 }
647 
649 void transVmetrictoV_ucov(FTYPE ROTANGLE, FTYPE *Vmetric, FTYPE*ucov)
650 {
651 
652  // get transformation from Vmetric to V (i.e. rotation from BH spin along z-axis to rotated in -x direction by ROTANGLE
653  FTYPE trans[NDIM][NDIM];
654  transVmetrictoV_trans(ROTANGLE, Vmetric, trans);
655 
656  // now perform transformation
657  FTYPE tmpucov[NDIM];
658  int j,l;
659  DLOOPA(j){
660  tmpucov[j] = 0.;
661  DLOOPA(l){
662  // u_{mup} = u_{mu} T^mu_mup
663  // where T^mu_mup == dx^mu[BL]/dx^mup[KSP uni grid]
664  // where T^\mu[Vmetric]_mup[V] =dx^\mu[z-axis is BH axis]_mup[tilted axis is BH axis]
665  tmpucov[j] += ucov[l] * trans[l][j];
666  }
667  }
668  DLOOPA(j){
669  ucov[j] = tmpucov[j];
670  }
671 
672  /* done! */
673 }
674 
675 
676 
682 // see metricrot.nb
683 void transVmetrictoV_trans(FTYPE ROTANGLE, FTYPE *Vmetric, FTYPE (*trans)[NDIM])
684 {
685 
686  FTYPE told, r, h, ph;
687  told=Vmetric[0]; r=Vmetric[1]; h=Vmetric[2]; ph=Vmetric[3];
688 
689  /* make transform matrix */
690  // order for trans is [ourmetric][bl]
691  // DLOOP(j,k) trans[j][k] = 0. ;
692  // DLOOPA(j) trans[j][j] = 1. ;
693 
694  FTYPE b0=ROTANGLE;
695 
696  extern FTYPE csc(FTYPE arg);
697  extern FTYPE cot(FTYPE arg);
698 
699  // transformation is from SPC[Vmetric] to SPC[V] using Vmetric values of coordinates to identify coordinate location.
700  trans[0][0]=1.;
701  trans[0][1]=0.;
702  trans[0][2]=0.;
703  trans[0][3]=0.;
704  trans[1][0]=0.;
705  trans[1][1]=1.;
706  trans[1][2]=0.;
707  trans[1][3]=0.;
708  trans[2][0]=0.;
709  trans[2][1]=0.;
710  trans[2][2]=pow(pow(cos(h)*cos(ph)*sin(b0) - 1.*cos(b0)*sin(h),2.) + pow(sin(b0),2.)*pow(sin(ph),2.),-1.)*pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),0.5)*(-1.*cos(h)*cos(ph)*sin(b0) + cos(b0)*sin(h));
711  trans[2][3]=-1.*sin(b0)*sin(ph);
712  trans[3][0]=0.;
713  trans[3][1]=0.;
714  trans[3][2]=csc(h)*pow(pow(cos(h)*cos(ph)*sin(b0) - 1.*cos(b0)*sin(h),2.) + pow(sin(b0),2.)*pow(sin(ph),2.),-1.)*pow(pow(cos(h)*sin(b0) - 1.*cos(b0)*cos(ph)*sin(h),2.) + pow(sin(h),2.)*pow(sin(ph),2.),0.5)*sin(b0)*sin(ph);
715  trans[3][3]=cos(b0) - 1.*cos(ph)*cot(h)*sin(b0);
716 
717 
718  /* done! */
719 }
720 
721 
722 
723 
724 
725 // all below stuff independent of metrics
726 
727 
728 
731 int pr2ucon(int whichvel, FTYPE *pr, struct of_geom *geom, FTYPE*ucon)
732 {
733  FTYPE others[NUMOTHERSTATERESULTS];
734 
735  if(whichvel==VEL4){
736  // here pr has true 4-velocities, as supplied by init.c
737  if (ucon_calc_4vel(pr, geom, ucon, others) >= 1) {
738  dualfprintf(fail_file, "pr2ucon(ucon_calc): space-like error: whichvel=%d\n",whichvel);
739  return(1);
740  }
741  }
742  else if(whichvel==VEL3){ // supplied vel's are 3-vels
743  if (ucon_calc_3vel(pr, geom, ucon, others) >= 1) {
744  dualfprintf(fail_file, "pr2ucon(ucon_calc): space-like error: whichvel=%d\n",whichvel);
745  return(1);
746  }
747  }
748  else if(whichvel==VELREL4){ // supplied vel's are relative 4-vels
749  if (ucon_calc_rel4vel(pr, geom, ucon, others) >= 1) {
750  dualfprintf(fail_file, "pr2ucon(ucon_calc): space-like error: whichvel=%d\n",whichvel);
751  return(1);
752  }
753  }
754  else{
755  dualfprintf(fail_file,"No such whichvel=%d\n",whichvel);
756  myexit(1);
757  }
758  return(0);
759 }
760 
761 
762 void mettometp(int ii, int jj, int kk, FTYPE*ucon)
763 {
764  int loc;
765 
766  loc=CENT;
767  mettometp_genloc(ii, jj, kk, loc, ucon);
768 
769 }
770 
772 void mettometp_genloc(int ii, int jj, int kk, int loc, FTYPE*ucon)
773 {
774  int j,k;
775  FTYPE idxdxp[NDIM][NDIM];
776  FTYPE tmp[NDIM];
777 
778  idxdxprim_ijk(ii, jj, kk, loc, idxdxp);
779 
780  // actually gcon_func() takes inverse of first arg and puts result into second arg.
781  // matrix_inverse(PRIMECOORDS, dxdxp,idxdxp);
782 
783  /* transform ucon */
784  // this is u^j = (iT)^j_k u^k, as in mettobl() above
785  DLOOPA(j) tmp[j] = 0.;
786  DLOOP(j,k) tmp[j] += idxdxp[j][k] * ucon[k];
787  DLOOPA(j) ucon[j] = tmp[j];
788 
789  // note that u_{k,BL} = u_{j,KSP} (iT)^j_k
790 
791  // note that u_{k,KSP} = u_{j,BL} T^j_k
792 
793  // note that u^{j,BL} = T^j_k u^{k,KSP} // (T) called ks2bl in grmhd-transforms.nb
794 
795  // note that u^{j,KSP} = (iT)^j_k u^{k,BL} // (iT) called bl2ks in grmhd-transforms.nb
796 
797  // So T=BL2KSP for covariant components and KSP2BL for contravariant components
798  // and (iT)=BL2KSP for contra and KSP2BL for cov
799 
800  // where here T=dxdxp and (iT)=idxdxp (tranpose inverse)
801 
802  /* done! */
803 }
804 
806 void mettometp_simple(FTYPE (*idxdxp)[NDIM], FTYPE*ucon)
807 {
808  int j,k;
809  FTYPE X[NDIM], V[NDIM];
810  FTYPE tmp[NDIM];
811 
812  /* transform ucon */
813  // this is u^j = (iT)^j_k u^k, as in mettobl() above
814  DLOOPA(j) tmp[j] = 0.;
815  DLOOP(j,k) tmp[j] += idxdxp[j][k] * ucon[k];
816  DLOOPA(j) ucon[j] = tmp[j];
817 }
818 
820 void metptomet_ucov_simple(FTYPE (*idxdxp)[NDIM], FTYPE*ucov)
821 {
822  int j,k;
823  FTYPE X[NDIM], V[NDIM];
824  FTYPE tmp[NDIM];
825 
826  /* transform ucov */
827  DLOOPA(j) tmp[j] = 0.;
828  DLOOP(j,k) tmp[j] += idxdxp[k][j] * ucov[k];
829  DLOOPA(j) ucov[j] = tmp[j];
830 }
831 
832 
833 void metptomet(int ii, int jj, int kk, FTYPE*ucon)
834 {
835  int loc;
836 
837  loc=CENT;
838  metptomet_genloc(ii, jj, kk, loc, ucon);
839 
840 }
841 
843 void metptomet_genloc(int ii, int jj, int kk, int loc, FTYPE*ucon)
844 {
845  int j,k;
846  FTYPE dxdxp[NDIM][NDIM];
847  FTYPE tmp[NDIM];
848 
849  dxdxprim_ijk(ii, jj, kk, loc, dxdxp);
850 
851  /* transform ucon */
852  // this is u^j = T^j_k u^k, as in above
853  DLOOPA(j) tmp[j] = 0.;
854  DLOOP(j,k) tmp[j] += dxdxp[j][k] * ucon[k];
855  DLOOPA(j) ucon[j] = tmp[j];
856 
857  /* done! */
858 }
859 
861 void metptomet_simple(FTYPE (*dxdxp)[NDIM], FTYPE*ucon)
862 {
863  int j,k;
864  FTYPE tmp[NDIM];
865 
866  /* transform ucon */
867  // this is u^j = T^j_k u^k, as in above
868  DLOOPA(j) tmp[j] = 0.;
869  DLOOP(j,k) tmp[j] += dxdxp[j][k] * ucon[k];
870  DLOOPA(j) ucon[j] = tmp[j];
871 
872  /* done! */
873 }
874 
876 void mettometp_ucov_simple(FTYPE (*dxdxp)[NDIM], FTYPE*ucov)
877 {
878  int j,k;
879  FTYPE tmp[NDIM];
880 
881  /* transform ucov */
882  DLOOPA(j) tmp[j] = 0.;
883  DLOOP(j,k) tmp[j] += dxdxp[k][j] * ucov[k];
884  DLOOPA(j) ucov[j] = tmp[j];
885 
886  /* done! */
887 }
888 
890 void metptomet_Tud(int ii, int jj, int kk, FTYPE (*Tud)[NDIM])
891 {
892  int j,k;
893  int alpha,beta;
894  FTYPE V[NDIM], X[NDIM];
895  FTYPE dxdxp[NDIM][NDIM];
896  FTYPE idxdxp[NDIM][NDIM];
897  FTYPE tmp[NDIM][NDIM];
898  void metptomet_simple_Tud(FTYPE (*dxdxp)[NDIM], FTYPE (*idxdxp)[NDIM], FTYPE (*Tud)[NDIM]);
899  int loc;
900 
901 
902  loc=CENT;
903  dxdxprim_ijk(ii, jj, kk, loc, dxdxp);
904  idxdxprim_ijk(ii, jj, kk, loc, dxdxp);
905 
906  metptomet_simple_Tud(dxdxp, idxdxp, Tud);
907 
908  /* done! */
909 }
910 
914 void metptomet_simple_Tud(FTYPE (*dxdxp)[NDIM], FTYPE (*idxdxp)[NDIM], FTYPE (*Tud)[NDIM])
915 {
916  int j,k;
917  int alpha,beta;
918  FTYPE tmp[NDIM][NDIM];
919 
920  // Tud^j_k [KS] = T^\beta_\alpha [KS'] ((Lambda)^{-1})^\alpha_k \Lambda^j_\beta
921  DLOOP(j,k) tmp[j][k] = 0.;
922  DLOOP(alpha,beta) DLOOP(j,k) tmp[j][k] += idxdxp[alpha][k] * dxdxp[j][beta] * Tud[beta][alpha];
923  DLOOP(j,k) Tud[j][k] = tmp[j][k];
924 
925  /* done! */
926 }
927 
929 void ucon2pr(int whichvel, FTYPE *ucon, struct of_geom *geom, FTYPE *pr)
930 {
931  FTYPE alphasq,gammaoalpha,beta[NDIM] ;
932  int j;
933 
934  if(whichvel==VEL4){
935  SLOOPA(j) pr[U1+j-1]=ucon[j];
936  }
937  if(whichvel==VEL3){
938  SLOOPA(j) pr[U1+j-1] = ucon[j] / ucon[TT];
939  }
940  else if(whichvel==VELREL4){
941  alphasq = 1./(-geom->gcon[GIND(TT,TT)]) ;
942  gammaoalpha = ucon[TT] ;
943 
944  SLOOPA(j) beta[j] = alphasq*geom->gcon[GIND(TT,j)] ;
945  SLOOPA(j) pr[U1+j-1] = ucon[j] + beta[j]*gammaoalpha ;
946  }
947 }
948 
950 int vcon2pr(int whichvel, FTYPE *vcon, struct of_geom *geom, FTYPE *pr)
951 {
952  FTYPE alphasq,gammaoalpha,beta[NDIM] ;
953  int j;
954  FTYPE ucon[NDIM];
955  FTYPE others[NUMOTHERSTATERESULTS];
956  FTYPE prlocal[NPR];
957 
958  SLOOPA(j) prlocal[U1+j-1]=vcon[j];
959 
960  if(whichvel==VEL4){
961  // vcon->ucon
962 
963  if(ucon_calc_3vel(prlocal, geom, ucon, others)>=1){
964  dualfprintf(fail_file, "vcon2pr(ucon_calc_3vel): space-like error\n");
965  return(1);
966  }
967 
968  // ucon->pr
969  // SLOOPA(j) pr[U1+j-1]=vcon[j]*ucon[TT];
970  // or just directly pr=ucon
971  SLOOPA(j) pr[U1+j-1]=ucon[j];
972  }
973  if(whichvel==VEL3){
974  // vcon=pr
975  SLOOPA(j) pr[U1+j-1] = vcon[j] ;
976  }
977  else if(whichvel==VELREL4){
978  // vcon->ucon
979  if(ucon_calc_3vel(prlocal, geom, ucon, others)>=1){
980  dualfprintf(fail_file, "vcon2pr(ucon_calc_3vel): space-like error\n");
981  return(1);
982  }
983 
984  // now go from ucon->pr
985  alphasq = 1./(-geom->gcon[GIND(TT,TT)]) ;
986  gammaoalpha = ucon[TT] ;
987 
988  SLOOPA(j) beta[j] = alphasq*geom->gcon[GIND(TT,j)] ;
989  SLOOPA(j) pr[U1+j-1] = ucon[j] + beta[j]*gammaoalpha ;
990  }
991  return(0);
992 }
993 
995 //
996 // below not used
997 
998 #define NORMALDENSITY 0
999 #define LOGDENSITY 1
1000 
1001 #define DENSITYTYPE LOGDENSITY
1002 
1003 
1005 void density2pr(FTYPE *density, FTYPE *pr)
1006 {
1007 
1008 #if(DENSITYTYPE==NORMALDENSITY)
1009  density[RHO]=pr[RHO];
1010  density[UU]=pr[UU];
1011 #elif(DENSITYTYPE==LOGDENSITY)
1012  density[RHO]=log(pr[RHO]);
1013  density[UU]=log(pr[UU]);
1014 #endif
1015 
1016 }
1017 
1019 void pr2density(FTYPE *pr, FTYPE *density)
1020 {
1021 
1022 #if(DENSITYTYPE==NORMALDENSITY)
1023  pr[RHO]=density[RHO];
1024  pr[UU]=density[UU];
1025 #elif(DENSITYTYPE==LOGDENSITY)
1026  pr[RHO]=exp(density[RHO]);
1027  pr[UU]=exp(density[UU]);
1028 #endif
1029 
1030 }
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1