Init functions

Creating equatorial Ts and ρ profiles ()

loading Ts and extending by a power-law

imp=Import["G:\\mathematica\\3DMHD\\a098\\Tsmap098.dat"];

lpl=ListLogLogPlot[Transpose[Transpose[imp][[1;;2]]],Joined->True,PlotStyle->Red];
pl=LogLogPlot[(10^12.22)(r^(-0.8)),{r,imp[[1,1]],imp[[200,1]]}];
Show[pl,lpl,PlotRange->All]

imp[[1,1]]

Creating x2 grid

old simulations till 200M radial grid

np=64;θxx=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.h->0.65/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]]);
Cos[θxx]
x2tab
(*ListPlot[Transpose[{x2tab,1/((π/2)(h-7(1-h)(((x2-1))^6)))/.h->0.65/.x2->x2tab}],PlotRange->All]
ListPlot[Transpose[{x2tab,1/(π+(1-h) π Cos[2 π x2])/.h->0.65/.x2->x2tab}]]*)

new simulations till 10^5M radial grid

np=128;θxx=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.h->1.00/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]]);
Cos[θxx]
x2tab

Creating r grid

old simulations till 200M radial grid

dir = NotebookDirectory[];
read11[fnum_] := (namefile ="E:\\GRMHD\\a0\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];
   Close[namefile];);

read11[6950]

stt=StringToStream[strrr];
Do[x1min=Read[stt,Number],{i,5}]
Do[dx1=Read[stt,Number],{i,3}]
Do[off=Read[stt,Number],{i,6}]
Close[stt];
rtab=off+Exp[x1min+dx1 Range[0,255]]

rtab[[60]]

strrr

{x1min,dx1,off}

Dealing with normalizations

n rhonor mp (cc^2)->((B^2)/(4π));  B=b*Bnor
n  rhonor ->number density
n  rhonor mp (cc^2)->energy density ->(b^2)->(B^2)/(4π) => Bnor=(Sqrt[4π rhonor mp])cc

Separating Tp and Te

sp=1;
atab={0.,0.5,0.7,0.9,0.98};nmtab={8,8,8,9,9};ncuttab={145,147,149,153,162};a=atab[[sp]];nm=nmtab[[sp]];ncut=ncuttab[[sp]];
xx=Import["E:\\GRMHD\\a0\\Tsmap0.dat"];
rhonor=1042075.25;heat=0.41929;
rhoout=130;rmax=3.4 (10^5);rcut=xx[[ncut,1]];
rhopo=-Log[rhoout/rhonor/xx[[ncut]][[3]]]/Log[rmax/rcut];
rhocon=rhonor*xx[[ncut,3]];
Ucon=xx[[ncut,2]];
ruru={rg->((2 6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2)),arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.67 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12))};
rategs=7.722 (10^17);(*g/s*)rategs=ratexx/.Solve[0.5==(ratexx/(2 cc mp (rg^2)(r^2)ρfun))/.ruru/.r->3][[1]];(*this is a self-consistency check for the radial velocity*)
rtab=xx[[All,1]];
(*rategs=-NIntegrate[ρρ U1 sqdet/.{lr->Log[rr-off],x2->θ2x2[θ]}/.r->rr,{θ,θmin,θmax}]ρnor (rg^2)cc mp/.ruru
Print["Accretion rate = ",rate=rategs(year/Msun)/.ruru,"  Msun/year"];*)
(*rategs=-(rcut^2)((dat070[[All,32,3]]+dat070[[All,33,3]])/2);*)
rntab=Table[rcut((((rcut/rmax)))^(-(i/100))),{i,100}];ρfun=rhocon((((r/rcut)))^(-rhopo));
vrvr=-(rategs/(2 cc mp (rg^2)(r^2)ρfun ))/.ruru;dcap=(2.-(0.49/(((0.7+ (cser^2)))^2)));
rhofun=rhocon((((r/rcut)))^(-rhopo));Tsfun=Ucon((((r/rcut)))^(-0.98));Tstab=xx[[1;;ncut,1;;2]];(*Tsinter=Exp[Interpolation[Log[Join[Tstab,Transpose[{rntab,Tsfun/.r->#&/@rntab}]]]][Log[r]]];*)
Tsinter=Interpolation[Join[Tstab,Transpose[{rntab,Tsfun/.r->#&/@rntab}]]][r];
sysm={tail=(Dts=D[Tsinter,r];νc=8.9 (10^(-11))(((((Te[r])/(3 (10^10)))))^(-(3/2)))((ρfun)/(10^7)))((2 mp (rg^3)(r^2)ρfun )/rategs)(Tp[r]-Te[r])+2((((Sqrt[Tp[r]])/((Sqrt[Tp[r]])+heat(Sqrt[Te[r]]))))-(1/dcap)(1((heat(Sqrt[Te[r]]))/((Sqrt[Tp[r]])+heat(Sqrt[Te[r]])))))Dts;Te'[r]==((2Dts-tail)/(dcap+1)),Tp'[r]==tail+((2Dts-tail)/(dcap+1))}/.{θe->((kb Te[r])/(me (cc^2))),cser->(Sqrt[((kb Te[r])/(me (cc^2)))]),cspr->(Sqrt[((kb Tp[r])/(mp (cc^2)))])}/.ruru;

Off[FindRoot::"lstol"];
tosol=Flatten[{sysm,{Tp[r]==Tsinter,Te[r]==Tsinter}/.r->rmax}];
solnew=NDSolve[tosol,{Tp,Te},{r,rtab[[1]],rmax}][[1]];
lpl00=LogLogPlot[{Tp[r]/.solnew,Te[r]/.solnew},{r,rtab[[nm]],rmax},FrameLabel->{{"T, K",None},{"r  measured in (Subscript[r, g])","Red - (Subscript[T, p]), blue - (Subscript[T, e]); dashed - old model; solid - new model"}},PlotStyle->{{AbsoluteThickness[2],Red},{AbsoluteThickness[2],Blue}},Frame->True,FrameStyle->AbsoluteThickness[2],Axes->False,ImageSize->650,BaseStyle->{FontFamily->"Arial",FontWeight->"Bold",FontSize->16}]
((Tp[6.])/(Te[6.]))/.solnew
((Tp[3])/(Te[3]))/.solnew

rtab11=xx[[1;;ncut,1]];
ρtab11=rhonor xx[[1;;ncut,3]];
rtab22=Table[rcut((((rmax/rcut)))^(i/500)),{i,500}];
ρtab22=ρfun/.r->rtab22;
rtab=Join[rtab11,rtab22];
ρtab=Join[ρtab11,ρtab22];
ListLogLogPlot[Transpose[{rtab,ρtab}],Joined->True]
Tetab=Te[rtab]/.solnew;
Tptab=Tp[rtab]/.solnew;
solful={Te->Interpolation[Transpose[{rtab/2,Tetab}]],Tp->Interpolation[Transpose[{rtab/2,Tptab}]],ρ->Interpolation[Transpose[{Log[rtab/2],ρtab}]]}
solful>>"C:\\study\\accretion\\! 2012 1spring umd\\NGC3115\\spin0solution.dat"

With cooling

; ((B^2)/(8π))==ρfun kB Tp;
(dTe/dr)=((dTe )/(vr dt ))=(Te/(vr trad))=(Te/(cc vrvr))((4σT γ (B^2))/(3 8π me cc))=(Te/(cc vrvr))((4σT γ ρfun kB Tp)/(3 me cc))=(Te/(cc vrvr))((4σT kB ρfun Tp)/(3 me cc))(((kB Te )/(me (cc^2))));
3ρfun kb (dTe/dt)==-Λ ((ρfun)^2) => (dTe/dr)=(dTe/(vr dt))=-(1/(cc vrvr))((Λ ρfun)/(3kb))

cc vrvr/.r->10/.ruru

imp=Import["C:\\study\\accretion\\! 2012 1spring umd\\NGC3115\\Sutherland_cooling.xls"][[1]];(*for nH=1; ne=1.2; (Subscript[log, 10])((Subscript[Λ, net])[erg (cm^3)(s^(-1))]); *)
Tex=imp[[-1,1]]+0.05;Λex=imp[[-1,-1]];
Tran=Range[Tex,Tex+6,0.1];
Λran=Λex +0.5(Tran-Tex);
Λtot=Join[imp,Transpose[{Tran,Λran}]];
ListPlot[Λtot,Joined->True]
Xcool=Interpolation[Λtot]

sp=1;
atab={0.,0.5,0.7,0.9,0.98};nmtab={8,8,8,9,9};ncuttab={145,147,149,153,162};a=atab[[sp]];nm=nmtab[[sp]];ncut=ncuttab[[sp]];
xx=Import["E:\\GRMHD\\a0\\Tsmap0.dat"];
rhonor=1042075.25;heat=0.41929;
rhoout=130;rmax=3.4 (10^5);rcut=xx[[ncut,1]];
rhopo=-Log[rhoout/rhonor/xx[[ncut]][[3]]]/Log[rmax/rcut];
rhocon=rhonor*xx[[ncut,3]];
Ucon=xx[[ncut,2]];
ruru={rg->((2 6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2)),rs->rg,arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.67 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12)),σT->((8π)/3)(((((ee^2)/(me (cc^2)))))^2)};
ruru=Thread[ruru[[All,1]]->(ruru[[All,2]]//.ruru)];
rategs=7.722 (10^17);(*g/s*)rategs=ratexx/.Solve[0.5==(ratexx/(2 cc mp (rg^2)(r^2)ρfun))/.ruru/.r->3][[1]];(*this is a self-consistency check for the radial velocity*)
rtab=xx[[All,1]];
(*rategs=-NIntegrate[ρρ U1 sqdet/.{lr->Log[rr-off],x2->θ2x2[θ]}/.r->rr,{θ,θmin,θmax}]ρnor (rg^2)cc mp/.ruru
Print["Accretion rate = ",rate=rategs(year/Msun)/.ruru,"  Msun/year"];*)
(*rategs=-(rcut^2)((dat070[[All,32,3]]+dat070[[All,33,3]])/2);*)
rntab=Table[rcut((((rcut/rmax)))^(-(i/100))),{i,100}];ρfun=rhocon((((r/rcut)))^(-rhopo));
vrvr=-(rategs/(2 cc mp (rg^2)(r^2)ρfun ))/.ruru;dcap=(2.-(0.49/(((0.7+ (cser^2)))^2)));
rhofun=rhocon((((r/rcut)))^(-rhopo));Tsfun=Ucon((((r/rcut)))^(-0.98));Tstab=xx[[1;;ncut,1;;2]];(*Tsinter=Exp[Interpolation[Log[Join[Tstab,Transpose[{rntab,Tsfun/.r->#&/@rntab}]]]][Log[r]]];*)
Tsinter=Interpolation[Join[Tstab,Transpose[{rntab,Tsfun/.r->#&/@rntab}]]][r];
sysm={tail=(Dts=D[Tsinter,r];νc=8.9 (10^(-11))(((((Te[r])/(3 (10^10)))))^(-(3/2)))((ρfun)/(10^7)))((2 mp (rg^3)(r^2)ρfun )/rategs)(Tp[r]-Te[r])+2((((Sqrt[Tp[r]])/((Sqrt[Tp[r]])+heat(Sqrt[Te[r]]))))-(1/dcap)(1((heat(Sqrt[Te[r]]))/((Sqrt[Tp[r]])+heat(Sqrt[Te[r]])))))Dts;
Te'[r]==((2Dts-tail)/(dcap+1))-2rg((Te[r])/(cc vrvr))((4σT kb ρfun Tp[r])/(3 me cc))(((kb Te[r] )/(me (cc^2))))-0.195((2rg)/(cc vrvr))((10^Xcool[Log10[Te[r]]]ρfun)/(3 kb))/.ruru,Tp'[r]==tail+((2Dts-tail)/(dcap+1))}/.{θe->((kb Te[r])/(me (cc^2))),cser->(Sqrt[((kb Te[r])/(me (cc^2)))]),cspr->(Sqrt[((kb Tp[r])/(mp (cc^2)))])}/.ruru;

Off[FindRoot::"lstol"];
tosol=Flatten[{sysm,{Tp[r]==Tsinter,Te[r]==Tsinter}/.r->rmax}];
solnew=NDSolve[tosol,{Tp,Te},{r,rtab[[1]],rmax}][[1]];
lpl00=LogLogPlot[{Tp[r]/.solnew,Te[r]/.solnew},{r,rtab[[nm]],rmax},FrameLabel->{{"T, K",None},{"r  measured in (Subscript[r, g])","Red - (Subscript[T, p]), blue - (Subscript[T, e]); dashed - old model; solid - new model"}},PlotStyle->{{AbsoluteThickness[2],Red},{AbsoluteThickness[2],Blue}},Frame->True,FrameStyle->AbsoluteThickness[2],Axes->False,ImageSize->650,BaseStyle->{FontFamily->"Arial",FontWeight->"Bold",FontSize->16}]
{((Tp[6.])/(Te[6.]))/.solnew,Te[6.]/.solnew}
((Tp[3])/(Te[3]))/.solnew

0.15

0.1

NGC3115

with cooling

0.001 cooling

LogLogPlot[-rg(r/(cc vrvr))((4σT kb ρfun Tp[r])/(3 me cc))(((kb Te[r] )/(me (cc^2))))/.ruru/.solnew,{r,2,(10^4)}]

{ρfun,rg,ρfun rg}/.ruru/.r->3(*for Sgr A*)
{10586.346475568054`,5.928888888888888`*^14,5.928888888888888`*^14 10586.346475568054`}

Checking accretion rate

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];
gmn=Simplify[gmn];
invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1];
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm
Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}};
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}](*/.{R0->0.3,h->0.65}*);
AMKS=Inverse[AMKS];
invnew=AMKS.invgmnKS.AMKS
(*Do[g[i-1,j-1]=gmn[[i,j]],{i,4},{j,4}];
Do[invg[i-1,j-1]=invgmn[[i,j]],{i,4},{j,4}];
(*gmn = (g^(μν)); invgmn = (Subscript[g, μν])*)ee[0]={u[0],u[1],u[2],u[3]};
ee[1]=(s/N2){-g[0,0]C1-g[0,3]C2,g[1,1]C0,0,-g[3,3]C2-g[0,3]C1};ee[3]=-(1/N1){l,0,0,1};ee[2]=(1/N3){1,(C1/C0)((g[0,0]-(l^2))/(g[1,1])),(Cmin/C3)-((u[1])/(u[2]))(C1/C0)((g[0,0]-(l^2))/(g[1,1])),-l};
ulo=invgmn.ee[0];
Do[uco[i-1]=ulo[[i]],{i,4}];
l=uco[3]/uco[1];C0=u[0]uco[0]+u[3]uco[3];C1=u[1]uco[0];C2=u[1]uco[3];C3=u[2]uco[1];Cmin=u[3]uco[3]-u[0]uco[0];
s=Sign[C0];j=FullSimplify[(Sqrt[-Det[invgmn]]),{r>2,a>0,π>θ>0}];
(*gmn = (g^(μν)); invgmn = (Subscript[g, μν])*)N1=(((invg[0,0](l^2)-2invg[0,3]l+invg[3,3]))^(1/2));N2=(((g[0,0](C1^2)+g[1,1](C0^2)+g[3,3](C2^2)+2g[0,3]C1 C2))^(1/2));N3=(((g[0,0]-2g[0,3]l+((((C1/C0)))^2)((((g[0,0]-(l^2)))^2)/(g[1,1]))+g[3,3](l^2)+g[2,2]((((Cmin/C3)-((u[1])/(u[2]))(C1/C0)((g[0,0]-(l^2))/(g[1,1]))))^2)))^(1/2));*)

FullSimplify[Sqrt[-Det[invnew]],{0<h<1,a>0,π>θ>0,r>R0>0,0<x2<1}]
Simplify[Sqrt[-Det[invnew]],(h+7 (1-h) (((-1+2 x2))^6))>0]

(π^2)(r-R0)(h + 7 (1 - h) (((-1 + 2 x2))^6))((r^2)+(a^2)((Cos[θ])^2))Sin[θ]dx2 dx3==π(r-R0)((r^2)+(a^2)((Cos[θ])^2))Sin[θ]dθ dx3==-π(r-R0)((r^2)+(a^2)((Cos[θ])^2))dCos[θ] dx3==-(r-R0)((r^2)+(a^2)((Cos[θ])^2))dCos[θ] dϕ==-((2π)/phlen)(r-R0)((r^2)+(a^2)((Cos[θ])^2))dCos[θ] (Underscript[Σ, ϕ])

FullSimplify[Sqrt[-Det[invgmnKS]],{π>θ>0,r>0,a>0}]

dθ/dx2==D[(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2),x2]

Intensity functions

Defining points on the picture plane

beta = atan2 (yg, xg)

Polarized spherical coordinates -> Stokes

(R={{1,0,0,0},{0,Cos[2χ[s]],Sin[2χ[s]],0},{0,-Sin[2χ[s]],Cos[2χ[s]],0},{0,0,0,1}})//MatrixForm
(mat={{aIc,aQc,0,aVc},{aQc,aIc,rVc,0},{0,-rVc,aIc,rQc},{aVc,0,-rQc,aIc}})//MatrixForm
(fs=FullSimplify[Transpose[R].mat.R])//MatrixForm
eqq=Thread[{f[0],f[1],f[2],f[3]}==det({jIc,jQc,0,jVc}.R/frsq-frcu/frsq fs.{y[0],y[1],y[2],y[3]})]/.{Sin[2 χ[s]]->sin2k,Cos[2χ[s]]->cos2k}

doub cosqs=cos(yyy[2]),sinqs=sin(yyy[2]),frsq=fr*fr,sinths=sin(yyy[3]),cosths=cos(yyy[3]);

rull={y[0]->yyy[0],y[1]->yyy[1]Sin[yyy[3]]Cos[yyy[2]],y[2]->yyy[1]Sin[yyy[3]]Sin[yyy[2]],y[3]->yyy[1]Cos[yyy[3]],f[1]->D[yyy[1][s]Sin[yyy[3][s]]Cos[yyy[2][s]],s],f[2]->D[yyy[1][s]Sin[yyy[3][s]]Sin[yyy[2][s]],s],f[3]->D[yyy[1][s]Cos[yyy[3][s]],s],f[0]->ff[0]};
finex=FullSimplify[Simplify[Solve[eqq/.rull/.{(yyy[1])'[s]->ff[1],(yyy[2])'[s]->ff[2],(yyy[3])'[s]->ff[3],yyy[1][s]->yyy[1],yyy[2][s]->yyy[2],yyy[3][s]->yyy[3]},{ff[0],ff[1],ff[2],ff[3]}][[1]]]/.{Cos[yyy[2]]->cosqs,Sin[yyy[2]]->sinqs,Cos[yyy[3]]->cosths,Sin[yyy[3]]->sinths,Csc[yyy[3]]->(1/sinths),Cot[yyy[3]]->(cosths/sinths)}]

ToString[InputForm[finex]]

finex[[3,1]]->rVc D[finex[[3,2]],rVc]/.{frcu->(fr^3),frsq->(fr^2)}

det<0, but integration along negative time => effectively det >0;
dyy[2]/dl -> rVc<0 => yy[2]<0
in[kk][1]-> cos[yy[2]]
in[kk][2]-> sin[yy[2]]

χ=ang[kk]=atan2(in[kk][2],in[kk][1])/2.==ArcTan[in[kk][1],in[kk][2]]/2

Simplify[ArcTan[Cos[yy[2]],Sin[yy[2]]]/2,{0<yy[2]<(π/2)}]

χ ∝ rVc ∝ RM/ν^2<0 (RM<0)
=> higher ν => larger χ

Check the integration routine for Stokes

maxy=2;
susu[nn_]:=(F[a_,b_]:=Cos[a b]Exp[(a^2)+0.1(b^2)+(1/((a^2)+0.25))];snxy=nn;
xtab=maxy Range[-1,1,(2./snxy)];
ftab=F[xtab,#]&/@xtab;
Dimensions[ftab];
ss=Sum[(ftab[[gx+0,gy+0]]+3ftab[[gx+1,gy+0]]+3ftab[[gx+2,gy+0]]+ftab[[gx+3,gy+0]])+3(ftab[[gx+0,gy+1]]+3ftab[[gx+1,gy+1]]+3ftab[[gx+2,gy+1]]+ftab[[gx+3,gy+1]])+3(ftab[[gx+0,gy+2]]+3ftab[[gx+1,gy+2]]+3ftab[[gx+2,gy+2]]+ftab[[gx+3,gy+2]])+
(ftab[[gx+0,gy+3]]+3ftab[[gx+1,gy+3]]+3ftab[[gx+2,gy+3]]+ftab[[gx+3,gy+3]]),{gx,1,snxy-2,3},{gy,1,snxy-2,3}]((((3/8)h))^2)(maxy^2)/.h->(2/snxy);
tt=(Sum[(ftab[[gx+0,gy+0]]+4ftab[[gx+1,gy+0]]+ftab[[gx+2,gy+0]])+4(ftab[[gx+0,gy+1]]+4ftab[[gx+1,gy+1]]+ftab[[gx+2,gy+1]])+(ftab[[gx+0,gy+2]]+4ftab[[gx+1,gy+2]]+ftab[[gx+2,gy+2]]),{gx,1,snxy-2,2},{gy,1,snxy-2,2}]((((h/3)))^2)+Sum[ftab[[gx+0,snxy+0]]+ftab[[gx+1,snxy+0]]+ftab[[gx+0,snxy+1]]+ftab[[gx+1,snxy+1]],{gx,1,snxy}]((((h/2)))^2)+Sum[ftab[[snxy+0,gy+0]]+ftab[[snxy+1,gy+0]]+ftab[[snxy+0,gy+1]]+ftab[[snxy+1,gy+1]],{gy,1,snxy-1}]((((h/2)))^2))(maxy^2)/.h->(2/snxy);
vv=Sum[ftab[[gx+0,gy+0]]+ftab[[gx+1,gy+0]]+ftab[[gx+0,gy+1]]+ftab[[gx+1,gy+1]],{gx,1,snxy},{gy,1,snxy}]((((h/2)))^2)(maxy^2)/.h->(2/snxy);)//Quiet;
susu[45];Print[InputForm[ss]," ",InputForm[tt]," ",InputForm[vv]]
NIntegrate[F[a,b],{a,-maxy,maxy},{b,-maxy,maxy}]

Clear[snxy];StringReplace[ToString[InputForm[{(ftab[[gx+0,gy+0]]+4ftab[[gx+1,gy+0]]+ftab[[gx+2,gy+0]])+4(ftab[[gx+0,gy+1]]+4ftab[[gx+1,gy+1]]+ftab[[gx+2,gy+1]])+(ftab[[gx+0,gy+2]]+4ftab[[gx+1,gy+2]]+ftab[[gx+2,gy+2]]),ftab[[gx+0,snxy+0]]+ftab[[gx+1,snxy+0]]+ftab[[gx+0,snxy+1]]+ftab[[gx+1,snxy+1]],ftab[[snxy+0,gy+0]]+ftab[[snxy+1,gy+0]]+ftab[[snxy+0,gy+1]]+ftab[[snxy+1,gy+1]]}]/.{A_[[n_,m_]]->ausin[n][m][il]}/.{gx->ix,gy->iy}],{"ausin"->"(*ausin)"}]

Check definitions of distance to Sgr A*, percentages, position angles

percentages OK

ruru={rg->(( 6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2))(*=1M, not 2M*),arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.26 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12))}
unit=(10^23)((((rg))^2)/(DD^2))/.DD->8.4 (10^3)pc/.ruru

(rg/(cc 60))/.ruru(*min in rg*)

Mod[(180/π)ArcTan[Q,U]/2+180,180]/.{U->-7.,Q->-244}(*compare with Marrone 2007*)
Mod[(180/π)ArcTan[Q,U]/2+180,180]/.{U->28,Q->-224.}(*Note that ArcTan[Q,U] is the same as atan2(U,Q) in C++ --- the arguments exchange*)
Mod[(180/π)ArcTan[Q,U]/2+180,180]/.{U->-176,Q->100.}
Mod[(180/π)ArcTan[Q,U]/2+180,180]/.{U->-14,Q->-12.}

atan2 calculates the arctangent of y/x
atan2( y, x)
ArcTan[x,y]
gives the arc tangent of (y/x), taking into account which quadrant the point (x,y) is in.

Check xisq definition & residual calculations

Check that file output is correct for poliresa

Geoint functions

Initializing the equations and defining quantities

Defining initial vector to certain accuracy

b->(Sqrt[(xpl^2)+(ypl^2)]);rA=(Sqrt[(r0^2)+(b^2)]);θA=ArcCos[((r0 Cos[θ0]-b Cos[β] Sin[θ0])/rA)];ϕA=ArcTan[r0 Sin[θ0]+b Cos[β] Cos[θ0],b Sin[β]];
initss=Solve[FullSimplify[Normal[Series[init={t[0]==0,ϕ[0]==ϕA,r[0]==rA,μ[0]==((r0 Cos[θ0]-b Cos[β] Sin[θ0])/rA),r'[0]==-r0,μ'[0]==-Sin[θA]((r0 (Sin[θA]-Cos[θA] Cos[ϕA] Tan[θ0]))/(rA (Cos[θA]+Cos[ϕA] Sin[θA] Tan[θ0])))},{b,0,2}]],{r0>0,π>β>0,π>θ0>0}],{t[0],ϕ[0],r[0],μ[0],r'[0],μ'[0]}][[1]]

StringReplace[ToString[InputForm[initss]],{"Cos[β]"->"cosbeta","Sin[β]"->"sinbeta","Sin[θ0]"->"sinth","r0^2"->"r0sq","b^2"->"bsq","Cos[θ0]"->"costh","*Csc[θ0]"->"/sinth","Cot[θ0]"->"costh/sinth"}]

FullSimplify[ArcTan[(b/a)]==ArcTan[a,b],{(π/2)>a>0,(π/2)>b>0}]

beta = atan2 (yg, xg) - corresponds to the picture
xpl=b Cos[β]; ypl=b Sin[β];rA Cos[θA] == r Cos[θ0] + ypl Sin[θ0] => θA=ArcCos[((r0 Cos[θ0]+ypl Sin[θ0])/rA)];
 Tan[ϕA]==xpl/(r Sin[θ0]-ypl Cos[θ0]) => ϕA=ArcTan[r0 Sin[θ0]-ypl Cos[θ0], xpl]; 
 
k={0, Sin[θA] Cos[ϕA] Sin[θ0] + Cos[θA] Cos[θ0], Cos[θA] Cos[ϕA] Sin[θ0] - Sin[θA] Cos[θ0], Sin[ϕA] Sin[θ0]};
(k^μ)={kt,1,0,0}; a0={0,0,-1,0}; b={0,0,0,-1}; 
z=r Cos[θ];
x=r Sin[θ]Sin[ϕ];
y=r Sin[θ]Cos[ϕ];
r=(Sqrt[(x^2)+(y^2)+(z^2)])=> dr== ((dx+dy+dz)/r); 
θ=ArcTan[(Sqrt[(x^2)+(y^2)])/z] => dθ ==
k={kt,0,Sin[θ],Cos[θ]}={kt,kx,ky,kz} at a point {x,y,z}=rA{ Sin[θA]Sin[ϕA], Sin[θA]Cos[ϕA],Cos[ϕA]}; line: ((x-xA)/kx)==((y-yA)/ky)==((z-zA)/kz)==t;

Clear[b,rA,β,xpl,ypl];rsub=(Sqrt[(x^2)+(y^2)+(z^2)]);θsub=ArcTan[z,(Sqrt[(x^2)+(y^2)])];ϕsub=ArcTan[x,y];
zsub=r0 Cos[θ];xsub=r0 Sin[θ]Cos[ϕ];ysub=r0 Sin[θ]Sin[ϕ];kxyz=-{Sin[θ0],0,Cos[θ0]};axyz={-Cos[θ0],0,Sin[θ0]};
Print["Checking substitutions: ",FullSimplify[{z==zsub,x==xsub,y==ysub}/.{r0->rsub,θ->θsub,ϕ->ϕsub}]]
krr=Thread[Simplify[{Dt[rsub],Dt[θsub],Dt[ϕsub]}/.{Dt[x]->kx,Dt[y]->ky,Dt[z]->kz}/.{z->zsub,x->xsub,y->ysub}/.{r0->rA,θ->θA,ϕ->ϕA},{rA>0,π>θA>0}]]/.Thread[{kx,ky,kz}->kxyz];dkrr=r0 Flatten[{1,krr}]
arr=Thread[Simplify[{Dt[rsub],Dt[θsub],Dt[ϕsub]}/.{Dt[x]->kx,Dt[y]->ky,Dt[z]->kz}/.{z->zsub,x->xsub,y->ysub}/.{r0->rA,θ->θA,ϕ->ϕA},{rA>0,π>θA>0}]]/.Thread[{kx,ky,kz}->axyz];darr=Flatten[{0,arr}]
(*FullSimplify[{r==rsub,θ==θsub,ϕ==ϕsub}/.{z->r Cos[θ],x->r Sin[θ]Cos[ϕ],y->r Sin[θ]Sin[ϕ]},{0<θ<π,r>0,-(π/2)<ϕ<(π/2)}]*)

rA=(Sqrt[(r0^2)+(b^2)]);θA=ArcCos[((r0 Cos[θ0]-b Cos[β] Sin[θ0])/rA)];ϕA=ArcTan[r0 Sin[θ0]+b Cos[β] Cos[θ0],b Sin[β]];
xkrr=FullSimplify[Normal[Series[{dkrr[[1]],dkrr[[2]],-Sin[θA]dkrr[[3]],dkrr[[4]]},{b,0,2}]],{r0>0,0<θ0<π}];
xarr=FullSimplify[Normal[Series[{darr[[1]],darr[[2]],-Sin[θA]darr[[3]],darr[[4]]},{b,0,2}]],{r0>0,0<θ0<π}];
StringReplace[ToString[InputForm[{tprin==xkrr[[1]],rprin==xkrr[[2]],muprin==xkrr[[3]],phprin==xkrr[[4]]}]],{"Cos[β]"->"cosbeta","Sin[β]"->"sinbeta","Sin[θ0]"->"sinth","r0^2"->"r0sq","b^2"->"bsq","Cos[θ0]"->"costh","*Csc[θ0]"->"/sinth","Cot[θ0]"->"costh/sinth"}]
StringReplace[StringReplace[ToString[InputForm[{tain==xarr[[1]],rain==xarr[[2]],muain==xarr[[3]],pain==xarr[[4]]}]],{"Cos[β]"->"cosbeta","Sin[β]"->"sinbeta","Sin[θ0]"->"sinth","r0^2"->"r0sq","r0^3"->"(r0sq*r0)","b^2"->"bsq","Cos[θ0]"->"costh","*Csc[θ0]"->"/sinth","Cot[θ0]"->"costh/sinth"}],"(2*(b - r0)*(b + r0) + bsq*Cos[2*β])*sinth)"->"(bsq-2*r0sq+2*bsq*cosbetasq)*sinth)"]
(*{t'[0]->dkrr[[1]],r'[0]->dkrr[[2]],μ'[0]->-Sin[θA]dkrr[[3]],ϕ'[0]->dkrr[[4]]}*)

Simplify[(2*(bsq - r0sq) + bsq*Cos[2*β])/.Cos[2β]->2((Cos[β])^2)-1]

Simplify[(b*costh/sinth*(r0-2*b*cosbeta*costh/sinth)/sinth*sinbeta)/(r0sq*r0)==(
(b*costh*sinbeta*(-2*b*cosbeta*costh+r0*sinth))/(r0sq*r0*(sinth^2)*sinth))]
Simplify[(2*(b - r0)*(b + r0) + bsq*Cos[2*β])*sinth==(bsq-2*r0sq+2*bsq*cosbetasq)*sinth/.{r0sq->(r0^2),bsq->(b^2),cosbetasq->((Cos[β])^2)}]
FullSimplify[(bsq*(-1+cosbetasq*(-1+2*costhsq))-2*(-1+costhsq)*r0sq)/(2*r0sq*r0*sinthsq)/((2*b*r0*cosbeta*costh - (bsq-2*r0sq+2*bsq*cosbetasq)*sinth)/(2*r0sq*r0))/.{cosbetasq->((Cos[β])^2),costhsq->((Cos[θ])^2),costh->Cos[θ],cosbeta->Cos[β],sinthsq->((Sin[θ])^2)}](*comparison to the incorrect previous result for muprin*)

-Sin[θ] θ'[0] == μ'[0]

diag=DiagonalMatrix[{1,rA,rA Sin[θA]}];ddiag=DiagonalMatrix[{I,1,rA,rA Sin[θA]}];
Simplify[dkrr.ddiag.ddiag.dkrr](*(k^μ)(Subscript[k, μ])==0*)
Simplify[krr.diag.diag.krr](*(Subscript[k, i])(Subscript[k, i])==1 - spatial*)
Simplify[arr.diag.diag.arr](*(Subscript[a, i])(Subscript[a, i])==1 - spatial*)
Simplify[arr.diag.diag.krr](*(Subscript[a, i])(Subscript[k, i])==0 - orthogonality*)
Simplify[krr/.θA->θ0/.ϕA->0](*radial vector*)

subs={θ0->1.1,r0->10};
b=(Sqrt[(xpl^2)+(ypl^2)]);β=ArcTan[xpl,ypl];
rA=(Sqrt[(r0^2)+(b^2)]);
θA=FullSimplify[ArcCos[((r0 Cos[θ0]-b Cos[β] Sin[θ0])/rA)]];
ϕA=FullSimplify[ArcTan[r0 Sin[θ0]+b Cos[β] Cos[θ0],b Sin[β]]];nn=5;xmax=3;
xtab=Range[-xmax,xmax,2(xmax/nn)];ytab=Range[-xmax,xmax,2(xmax/nn)];
lpoints=ListPointPlot3D[Join[{{0,0,0}},coords=Flatten[Transpose[({xsub,ysub,zsub}/.{r0->rA,θ->θA,ϕ->ϕA}/.subs/.xpl->xtab)/.ypl->ytab,{3,2,1}],1]],BoxRatios->{1,1,1},ViewPoint->{1000,1000,1000},PlotStyle->Black,ImageSize->500,PlotStyle->PointSize[Large]];

(*vect=Transpose[(krr/.subs/.xpl->Range[-xmax,xmax,2(xmax/nn)])/.ypl->Range[-xmax,xmax,2(xmax/nn)],{3,2,1}];*)len=0.7;
lkvect=Graphics3D[{Blue,Line[Flatten[Table[{{xsub,ysub,zsub}/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]
,{xsub,ysub,zsub}/.Thread[{r0,θ,ϕ}->({rA,θA,ϕA}+len krr/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]})]}/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},
{i,Length[xtab]},{k,Length[ytab]}],1]]},BoxRatios->{1,1,1},Axes->True,PlotRange->All];
lavect=Graphics3D[{Blue,Line[Flatten[Table[{{xsub,ysub,zsub}/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]
,{xsub,ysub,zsub}/.Thread[{r0,θ,ϕ}->({rA,θA,ϕA}+len arr/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]})]}/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},
{i,Length[xtab]},{k,Length[ytab]}],1]]},BoxRatios->{1,1,1},Axes->True,PlotRange->All];
gkvect=Graphics3D[{Red,Line[Flatten[Table[{{xsub,ysub,zsub},{xsub,ysub,zsub}+len kxyz}/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},{i,Length[xtab]},{k,Length[ytab]}],1]]},BoxRatios->{1,1,1},Axes->True,PlotRange->All];
gavect=Graphics3D[{Red,Line[Flatten[Table[{{xsub,ysub,zsub},{xsub,ysub,zsub}+len axyz}/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},{i,Length[xtab]},{k,Length[ytab]}],1]]},BoxRatios->{1,1,1},Axes->True,PlotRange->All];
Show[lpoints,lkvect,gkvect,lavect,gavect]

Total[Abs[Flatten[Table[kxyz/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},{i,Length[xtab]},{k,Length[ytab]}],1].
Transpose[Flatten[Table[axyz/.Thread[{r0,θ,ϕ}->{rA,θA,ϕA}]/.subs/.{xpl->xtab[[i]],ypl->ytab[[k]]},{i,Length[xtab]},{k,Length[ytab]}],1]]],2]

Correct Dsolve, definition of lmin, lmax, interpolation

lsol=Simplify[TrigToExp[Solve[1+3*Sinh[llog]/r0==t,llog][[1]]],t>0]
lsol/.t->0
FullSimplify[(llog/.lsol/.t->0)==-Log[(r0/3)+(Sqrt[1+((r0^2)/9)])],r0>0]

Geodes functions

Definitions of auxilliary and Christoffel symbols

polarization basis propagation

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);ρmin=(r^2)-(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);gmn={{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}};
invgmn=FullSimplify[Inverse[gmn]];
Print["Upper index (g^(μν))=",gmn//MatrixForm]
Print["Lower index (Subscript[g, μυ])=",invgmn//MatrixForm](*lower indices*)
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
Gtab=Table[G[i,k,l]=FullSimplify[TrigExpand[FullSimplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]]/.{(r^2)+(a^2) ((Cos[θ])^2)->ρxsq,(a^2)+2 (r^2)+(a^2) Cos[2 θ]->2ρxsq,(a^2)+(-2+r) r->Δx,(a^2)-2 (r^2)+(a^2) Cos[2 θ]->-2ρxmin,(a^2)+(r^2)->rasq}]/.{Cos[θ]->μ,Sin[θ]->(Sqrt[1-(μ^2)]),Cot[θ]->(μ/(Sqrt[1-(μ^2)]))}]/.a (r^2)-(a^3) (μ^2)->a ρxmin/.{1/(Sqrt[1-(μ^2)])->(sinmu^(-1)),(Sqrt[1-(μ^2)])->sinmu,(((1-(μ^2)))^(3/2))->sinmu sinmusq,(-1+(μ^2))->-sinmusq},{i,1,4},{k,1,4},{l,1,4}]
extab={Δ,ρsq,ρmin,Σsq}

litt=Complement[Tuples[{1,2,3,4},{3}],pos0=Position[Gtab,0]]
zerrul=Thread[gg@@@pos0->0];
uniq=Select[litt,#[[2]]≤#[[3]]&]
compl=Complement[litt,uniq]
StringReplace[ToString[(sublg=Thread[(gg@@@compl)->gg@@@({#[[1]],#[[3]],#[[2]]}&/@compl)])/.gg[i_,j_,k_]->gg[i-1][j-1][k-1]],{"["->"","]"->"",","->","," "->"","->"->"="}]
StringReplace[StringReplace[StringReplace[ToString[InputForm[Thread[(gg@@@uniq)->(G@@@uniq)]/.gg[i_,j_,k_]->gg[i-1][j-1][k-1]]],{"["->"","]"->"",","->","," "->"","->"->"="}],{"ρxsq"->"rhosq","Δx"->"dx","ρxmin"->"rhomin"}],{"μ^4"->"musq*musq","μ^2"->"musq","μ"->"y[5]","r^2"->"rsq","r^4"->"rsq*rsq","r^5"->"rsq*rsq*r","a^2"->"asq","a^3"->"asq*a","a^4"->"asq*asq","rhosq^2"->"rhosqsq","rhosq^3"->"(rhosqsq*rhosq)","rasq^2"->"rasq*rasq","r*"->"y[4]*","r/"->"y[4]/","r)"->"y[4])"}]

eqs=StringReplace[StringReplace[ToString[InputForm[Expand[Solve[Table[0==(∂_({λ,2}))x[k][λ]+Sum[(gg[k,σ,ν]/.{μ->μ[λ],θ->θ[λ],ϕ->ϕ[λ],t->t[λ],r->r[λ]})(∂_λ)x[σ][λ](∂_λ)x[ν][λ],{σ,1,4},{ν,1,4}],{k,1,4}]/.{θ''[λ]->-((μ[λ] ((μ'[λ])^2))/(sinmu^3))-((μ''[λ])/sinmu),θ'[λ]->-((μ'[λ])/sinmu)}/.sublg/.zerrul,{t''[λ],r''[λ],μ''[λ],ϕ''[λ]}]][[1]]/.{(sinmu^(-2))->(((1-musq))^(-1))}/.{r'[λ]->y[2],μ'[λ]->y[3],t'[λ]->y[6],ϕ'[λ]->y[7],t''[λ]->f[6],ϕ''[λ]->f[7],r''[λ]->f[2],μ''[λ]->f[3],μ[λ]->y[5]}/.{gg[o_,p_,s_]->gg[o-1][p-1][s-1]}]],{"y[2]^2"->"y[2]*y[2]","y[6]^2"->"y[6]*y[6]","y[7]^2"->"y[7]*y[7]","y[3]^2"->"y[3]*y[3]"," "->"","->"->"=","]["->"","gg["->"gg"}],{"23]"->"23","13]"->"13","02]"->"02","01]"->"01","11]"->"11","03]"->"03","12]"->"12","22]"->"22","33]"->"33"}]
(*//y:t, phi,  rp,mup,  r,mu, tpr,ϕpr
  //f:t',phi',r'',mu'',r',mu',t'',ϕ''
*)

uniq

D[-((μ'[λ])/(Sin[θ[λ]])),λ]/.{Csc[θ[λ]]->(1/sinmu),θ'[λ]->-((μ'[λ])/sinmu),Cot[θ[λ]]->((μ[λ])/sinmu)}

Clear[gg];eq11=Thread[Table[eph11[i]'[λ],{i,1,4}]->-Table[Sum[(gg[i,j,k]eph11[j][λ]/.{θ->θ[λ],r->r[λ]})x[k]'[λ]/.θ'[λ]->-((μ'[λ])/(Sqrt[1-(μ^2)])),{j,1,4},{k,1,4}],{i,1,4}]]/.sublg/.zerrul/.{(eph11[i_])'[λ]->f[i+7],eph11[i_][λ]->y[i+7],t'[λ]->y[6],ϕ'[λ]->y[7], μ'[λ]->y[3],r'[λ]->y[2],1/(Sqrt[1-(μ^2)])->(sinmu^(-1))}/.gg[i_,j_,k_]->gg[i-1][j-1][k-1];
StringReplace[StringReplace[ToString[InputForm[eq11]],{" "->"","->"->"=","]["->"","gg["->"gg"}],{"]+"->"+","])"->")","]-"->"-",","->";"}]
(*eph11[i_][λ]->y[i+7],...,eph11[4][λ]->y[7]*)
(*//y:t, phi,  rp,mup,  r,mu, tpr,ϕpr
  //f:t',phi',r'',mu'',r',mu',t'',ϕ''
*)

ComplexExpand[(A-I B)(r-I a Cos[θ])]
K2->A r-a B Cos[θ]
K1->(B r+a A Cos[θ])
A->kBL[0]*e1BL[1]-kBL[1]*e1BL[0]+a*(kBL[1]*e1BL[3]-kBL[3]*e1BL[1])*sinsq;
B->((rr*rr+asq)*(kBL[3]e1BL[2]-kBL[2]e1BL[3])-a(kBL[0]e1BL[2]-kBL[2]e1BL[0]))*sqrt(sinsq)

t->0.9999
r->9.
K2->-704470.42291774042
K1->636218.82084316690

r->19
K2->-705233.30073512404
K1->636206.49786268233

r->9.
K2->-704470.42291774042
K1->636218.82084316690

r->62
K1->26187649.977335777
K2->-27847388.231093597
r->207.
K1->26187863.498674039
K2->-28122333.039226286

Solvetrans functions

Correct initialization

Correct initial II and convergent solutions

Transnew functions

Correct auxilliary definitions and equation

(R={{1,0,0,0},{0,Cos[2χ[s]],-Sin[2χ[s]],0},{0,Sin[2χ[s]],Cos[2χ[s]],0},{0,0,0,1}})//MatrixForm
(mat={{aIc,aQc,0,aVc},{aQc,aIc,rVc,0},{0,-rVc,aIc,rQc},{aVc,0,-rQc,aIc}})//MatrixForm
(fs=FullSimplify[R.mat.Transpose[R]])//MatrixForm
eqq=Thread[{f[0],f[1],f[2],f[3]}==det(R.{jIc,jQc,0,jVc}/frsq-frcu/frsq fs.{y[0],y[1],y[2],y[3]})]/.{Sin[2 χ[s]]->sin2k,Cos[2χ[s]]->cos2k}

doub cosqs=cos(yyy[2]),sinqs=sin(yyy[2]),frsq=fr*fr,sinths=sin(yyy[3]),cosths=cos(yyy[3]);

rull={y[0]->yyy[0],y[1]->yyy[1]Sin[yyy[3]]Cos[yyy[2]],y[2]->yyy[1]Sin[yyy[3]]Sin[yyy[2]],y[3]->yyy[1]Cos[yyy[3]],f[1]->D[yyy[1][s]Sin[yyy[3][s]]Cos[yyy[2][s]],s],f[2]->D[yyy[1][s]Sin[yyy[3][s]]Sin[yyy[2][s]],s],f[3]->D[yyy[1][s]Cos[yyy[3][s]],s],f[0]->ff[0]};
finex=FullSimplify[Simplify[Solve[eqq/.rull/.{(yyy[1])'[s]->ff[1],(yyy[2])'[s]->ff[2],(yyy[3])'[s]->ff[3],yyy[1][s]->yyy[1],yyy[2][s]->yyy[2],yyy[3][s]->yyy[3]},{ff[0],ff[1],ff[2],ff[3]}][[1]]]/.{Cos[yyy[2]]->cosqs,Sin[yyy[2]]->sinqs,Cos[yyy[3]]->cosths,Sin[yyy[3]]->sinths,Csc[yyy[3]]->(1/sinths),Cot[yyy[3]]->(cosths/sinths)}]

D[ff[2]/.finex,rVc]/.{frcu->(fr^3),frsq->(fr^2)}

{3.200,0.1665,0.060722,0.53051} at coeff=20;EVPA=24.0deg;
{2.6999,-0.29622,-0.127,0.3655}at coeff=10 ;EVPA=19.7deg;
magnetic field is in the opposite direction (kbpar<0, rVc<0,)
{-3.58,1.3375 (10^(-5)),0.23816}
{-3.67,1.3684 (10^(-5)),0.22091}
{-3.76,1.3920 (10^(-5)),0.210575} =>II[2] decrease along the ray
{1.433665 (10^(-5)),0.166509} at coeff=20;
{1.3937 (10^(-5)), 0.20823} at coeff=21;
{1.34405 (10^(-5)),0.26189} at coeff=22;
stronger B field extension: Bpo=2 => 24deg
stronger B field extension: Bpo=1.9 => 26.2deg
stronger B field extension: Bpo=1.8 => 20.9deg
stronger B field extension: Bpo=1.7 => 19.9deg
stronger B field extension: Bpo=1.6 => 17.0deg
stronger B field extension: Bpo=1.5 => 10.5deg - correct behavior! => magnetic field reversal is present

Larger ff[2] from larger kbpar, larger coeff at rVc
coeff=1:   {{4.73058, 1.0154}, {4.260062, 0.89517}, {1.801617, 1.11190}, {5.46000, 1.59975}} => 
coeff=1.1:{{4.3664,   1.0138}, {5.090768, 1.01840}, {1.905887, 1.10751}, {5.506377, 1.5992}}
ausin[1]==Cos[ff[2]]
ausin[2]== Sin[ff[2]]
ang==atan2[in[2],in[1]]/2==atan[Sin[ff[2]]/Cos[ff[2]]]/2=ff[2]/2.

coef=1.05: {{4.91083,1.01699},{4.3144,0.960277},   {1.85309,1.10971},    {5.49337, 1.600488}}

ToString[InputForm[finex]]

Symmetries

For  B->-B: ;
θ->π-θ ; (Subscript[ϵ, I])->(Subscript[ϵ, I]); (Subscript[ϵ, Q])->(Subscript[ϵ, Q]); (Subscript[ϵ, V])->-(Subscript[ϵ, V]); (Subscript[α, V])->-(Subscript[α, V]); (Subscript[ρ, V])->-(Subscript[ρ, V]);(Subscript[ρ, Q])->(Subscript[ρ, Q])
(e^2)->-(e^2) =>χ->χ, but transformation  is  different

matrot={{Cos[χ],Sin[χ]},{-Sin[χ],Cos[χ]}};inv={{1,0},{0,-1}};
Thread[{a,b}=={e1,e2}.matrot]
Thread[{a,b}=={e1,-e2}.inv.matrot]
Simplify[Inverse[matrot.matrot]]//MatrixForm
Simplify[Inverse[(inv.matrot).(inv.matrot)]]

RRx={{1,0,0,0},{0,Cos[2χ],Sin[2χ],0},{0,Sin[2χ],-Cos[2χ],0},{0,0,0,-1}};
Simplify[RRx==Inverse[RRx]]
RR={{1,0,0,0},{0,Cos[2χ],-Sin[2χ],0},{0,Sin[2χ],Cos[2χ],0},{0,0,0,1}};
InvRR=Simplify[Inverse[RR]];Simplify[InvRR==(RR/.χ->-χ)]

Print["regular 0<θ<(π/2)"]
RR.{ϵI,ϵQ,0,ϵV}//MatrixForm
Simplify[RR.{{αI,αQ,0,αV},{αQ,αI,ρV,0},{0,-ρV,αI,ρQ},{αV,0,-ρQ,αI}}.InvRR]//MatrixForm

Print["Inverted magnetic field B-> -B"]
RRx.{ϵI,ϵQ,0,ϵV}//MatrixForm
Simplify[RRx.{{αI,αQ,0,αV},{αQ,αI,ρV,0},{0,-ρV,αI,ρQ},{αV,0,-ρQ,αI}}.RRx]//MatrixForm

Thus, ϵV -> -ϵV; αV-> -αV; ρV -> -ρV - everything else is identical!!! => eqs. of radiative transfer are not symmetric w/ respect to change of B direction
(χ - same angle)

Emis functions

Computing emissivities/rotativities

ttx=SessionTime[];Clear[inum,jnum,θetab,νΩtab,θe,νΩ,partab];
inum=100;jnum=60;
θetab=Table[0.4 (1.05^i),{i,0,inum}];νΩtab=Table[12000 (1.1^j),{j,-jnum,jnum}];
SetSharedVariable[inum,jnum,θetab,νΩtab];
lookuprV=Log[Table[((BesselK[0,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{i,Length[θetab]},{j,Length[νΩtab]}]];
lookuprQ=Log[Table[((BesselK[1,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))+6θe/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{i,Length[θetab]},{j,Length[νΩtab]}]];
lookupjI=Log[ParallelTable[Quiet[NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{γ,1,3000},{t,((νΩ)/(γ^2))/.{νΩ->νΩtab[[j]]},+∞}]],{i,Length[θetab]},{j,Length[νΩtab]}]];
lookupjQ=Log[ParallelTable[Quiet[NIntegrate[((νΩ)/(γ^2))BesselK[(2/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{γ,1,3000}]],{i,Length[θetab]},{j,Length[νΩtab]}]];
lookupjV=Log[ParallelTable[Quiet[NIntegrate[((νΩ)/(γ^3))BesselK[(1/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{γ,1,3000}]+NIntegrate[((νΩ)/(γ^3))BesselK[(1/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{γ,1,3000},{t,((νΩ)/(γ^2))/.{νΩ->νΩtab[[j]]},+∞}]],{i,Length[θetab]},{j,Length[νΩtab]}]];
Print[SessionTime[]-ttx]
(*149s for 10x12 table \[Implies] 7450s for 100x60 table*)

Export[NotebookDirectory[]<>"lookupjIy.dat",lookupjI,"Table"]
Export[NotebookDirectory[]<>"lookupjQy.dat",lookupjQ,"Table"]
Export[NotebookDirectory[]<>"lookupjVy.dat",lookupjV,"Table"]
Export[NotebookDirectory[]<>"lookuprVa.dat",lookuprQ,"Table"]
Export[NotebookDirectory[]<>"lookuprQa.dat",lookuprV,"Table"]

Correct lookup of emissivities/rotativities

Clear[dx,dy];ta={{{0,0},f[0,0]},{{1,1},f[1,1]},{{1,0},f[1,0]},{{0,1},f[0,1]}};
fuss=InterpolatingPolynomial[ta,{lTeman,lnWman}];
InputForm[FullSimplify[fuss/.{f[0,0]->jI[j][k],f[1,0]->jI[j+1][k],f[0,1]->jI[j][k+1],f[1,1]->jI[j+1][k+1]}]]
(*ToString[InputForm[FullSimplify[fuss]/.f[o_,m_]:>xx[k+o][j+m][i]]]
ToString[InputForm[FullSimplify[Collect[fuss,f[___]]]/.f[o_,m_]:>xx[k+o][j+m][i]]]*)

temp=(1+lnWman*(-1+lTeman)-lTeman);
Print[InputForm[FullSimplify[temp*jI[j][k]+lTeman*jI[1+j][k]+lnWman*(jI[j][1+k]+lTeman*(-jI[1+j][k]-jI[j][1+k]+jI[1+j][1+k]))]]]
temp=lnWman*lTeman;
InputForm[FullSimplify[(-lnWman-lTeman+1+temp)*jI[j][k]+(lnWman-temp)*jI[j][1+k]+(lTeman-temp)*jI[1+j][k]+temp*jI[1+j][1+k]]]

Evalpointnew (obsolete) functions

Lookup of coordinate, k, e1 vectors in BL metric

Definition of MKS metric (see GRMHD analysis shots.nb)

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];
gmn=Simplify[gmn];
invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1];
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm
Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}};
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}];
(*AMKS=Inverse[AMKS];*)
invnewxx=Inverse[AMKS].invgmnKS.Inverse[AMKS]
invnew=invnewxx/.{R0->0.3,h->0.65};
(*Do[g[i-1,j-1]=gmn[[i,j]],{i,4},{j,4}];
Do[invg[i-1,j-1]=invgmn[[i,j]],{i,4},{j,4}];
(*gmn = (g^(μν)); invgmn = (Subscript[g, μν])*)ee[0]={u[0],u[1],u[2],u[3]};
ee[1]=(s/N2){-g[0,0]C1-g[0,3]C2,g[1,1]C0,0,-g[3,3]C2-g[0,3]C1};ee[3]=-(1/N1){l,0,0,1};ee[2]=(1/N3){1,(C1/C0)((g[0,0]-(l^2))/(g[1,1])),(Cmin/C3)-((u[1])/(u[2]))(C1/C0)((g[0,0]-(l^2))/(g[1,1])),-l};
ulo=invgmn.ee[0];
Do[uco[i-1]=ulo[[i]],{i,4}];
l=uco[3]/uco[1];C0=u[0]uco[0]+u[3]uco[3];C1=u[1]uco[0];C2=u[1]uco[3];C3=u[2]uco[1];Cmin=u[3]uco[3]-u[0]uco[0];
s=Sign[C0];j=FullSimplify[(Sqrt[-Det[invgmn]]),{r>2,a>0,π>θ>0}];
(*gmn = (g^(μν)); invgmn = (Subscript[g, μν])*)N1=(((invg[0,0](l^2)-2invg[0,3]l+invg[3,3]))^(1/2));N2=(((g[0,0](C1^2)+g[1,1](C0^2)+g[3,3](C2^2)+2g[0,3]C1 C2))^(1/2));N3=(((g[0,0]-2g[0,3]l+((((C1/C0)))^2)((((g[0,0]-(l^2)))^2)/(g[1,1]))+g[3,3](l^2)+g[2,2]((((Cmin/C3)-((u[1])/(u[2]))(C1/C0)((g[0,0]-(l^2))/(g[1,1]))))^2)))^(1/2));*)

invnewxx//MatrixForm
(matr={{ig[0][0],ig[0][1],0,ig[0][3]},{ig[0][1],ig[1][1],0,ig[1][3]},{0,0,ig[2][2],0},{ig[0][3],ig[1][3],0,ig[3][3]}})//MatrixForm

Clear[off];InputForm[Delete[Union[Thread[Flatten[matr]->Flatten[invnewxx]]]/.{R0->off,(r^2)->rsq,(a^2)->asq,((Cos[θ])^2)->cossq,((Sin[θ])^2)->sinsq,h->h1}/.{asq cossq+rsq->rhosq}/.((A_)/(asq cossq+rsq))->A/rhosq/.r->rr,1]]

Definition and lookup of quantities on the grid

Correct power-law extension of B etc.

dir="G:\\mathematica\\3DMHD\\a";atab={"0","05","07","09","098"};
read11[n_]:=(namefile=dir<>atab[[n]]<>"\\fieldline00-1.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)
res2m=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*32*64*256],11],256],64],{3,2,1,4}];
Close[namefile];
namefile =dir<>atab[[n]]<>"\\fieldline0000.bin";
   siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   (*res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];*)
   Close[namefile];);
read11[5]
stt=StringToStream[strrr];
Do[x1min=Read[stt,Number],{i,5}]
Do[dx1=Read[stt,Number],{i,3}]
Do[off=Read[stt,Number],{i,6}]
Close[stt];
rtab=off+Exp[x1min+dx1 Range[0,255]]

(*{x1min,x2min,x3min}={BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"]};
{dx1,dx2,dx3}={BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"]};
Do[BinaryRead[fi11,"Real32"],{i,8}];BinaryRead[fi11,"Real32"];
ress=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",30*32*64*256],30],256],64],{3,2,1,4}];
Close[namefile];*)

Dimensions[avgres=((Total[res2m[[All,All,All,All]],{3}])/32)]
(*uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]*)
lpl3=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,6]]}]]];pl3=ListLogLogPlot[Transpose[{rtab,0.7(rtab^(-2.5))}],PlotStyle->Red];
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,8]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,0.3(rtab^(-1.5))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,7]]}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.03(rtab^(-1.5))}],PlotStyle->Red];
GraphicsRow[{Show[lpl3,pl3,PlotLabel->"(u^r)-upper index"],Show[lpl1,pl1,PlotLabel->"(u^ϕ)-upper index"],Show[lpl2,pl2,PlotLabel->"(u^θ)-upper index"]},ImageSize->1200]
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,9]]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,0.8(rtab^(-2.5))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,10]]]}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.01(rtab^(-1.5))}],PlotStyle->Red];
lpl3=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,11]]]}]]];pl3=ListLogLogPlot[Transpose[{rtab,0.8(rtab^(-2.5))}],PlotStyle->Red];
GraphicsRow[{Show[lpl1,pl1,PlotLabel->"(B^r)-upper index"],Show[lpl2,pl2,PlotLabel->"(B^θ)-upper index"],Show[lpl3,pl3,PlotLabel->"(B^ϕ)-upper index"]},ImageSize->1200]
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,1]]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,6(rtab^(-0.8))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,((avgres[[All,34,2]])/(avgres[[All,34,1]]))}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.25(rtab^(-0.8))}],PlotStyle->Red];
GraphicsRow[{Show[lpl1,pl1,PlotLabel->"ρ"],Show[lpl2,pl2,PlotLabel->"U"]},ImageSize->800]

Lookup of Te, u[0]

(matr={{ig[0][0],ig[0][1],0,ig[0][3]},{ig[0][1],ig[1][1],0,ig[1][3]},{0,0,ig[2][2],0},{ig[0][3],ig[1][3],0,ig[3][3]}})//MatrixForm

Clear[ulo];StringReplace[ToString[InputForm[Thread[{ulo[0],ulo[1],ulo[2],ulo[3]}==matr.{u[0],u[1],u[2],u[3]}]]],{"=="->"=",","->";"}]

Clear[u];InputForm[FullSimplify[Solve[-1=={u[0],u[1],u[2],u[3]}.matr.{u[0],u[1],u[2],u[3]},{u[0]}]][[1]]]

BL to MKS transformation

Old simulations

strrr

dir="E:\\GRMHD\\a09\\fieldlines\\";
read11[]:=(namefile=dir<>"fieldline6950.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
res2m=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*32*64*256],11],256],64],{3,2,1,4}];
Close[namefile];);
read11[]
stt=StringToStream[strrr];
Do[nx1=Read[stt,Number],{i,2}];Do[nx2=Read[stt,Number],{i,1}];Do[nx3=Read[stt,Number],{i,1}];{nx1,nx2,nx3}
Do[x1min=Read[stt,Number],{i,1}];
Do[dx1=Read[stt,Number],{i,3}];
Do[off=Read[stt,Number],{i,6}];
Do[h1=Read[stt,Number],{i,3}];
Close[stt];{x1min,dx1,off,h1}
rtab=off+Exp[x1tab=x1min+dx1 Range[0,255]];
np=nx2;h=h1;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]]);
(*read22[]:=(namefilexx=dir<>"usdump0000";
fi22=OpenRead[namefilexx,BinaryFormat->True];
BinaryRead[fi22,"Integer32"];
tt=BinaryRead[fi22,"Real32"];
Do[Print[BinaryRead[fi22,"Integer32"]],{i,3}];
{x1min,x2min,x3min}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
{dx1,dx2,dx3}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
Do[Print[BinaryRead[fi22,"Real32"]],{i,8}];BinaryRead[fi22,"Real32"];
ress=Transpose[Partition[Partition[Partition[BinaryReadList[fi22,"Real32",30*32*64*256],30],256],64],{3,2,1,4}];);
read22[]*)

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];
gmn=Simplify[gmn];
invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1];
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm
Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}};
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}];
invnew=Inverse[AMKS].invgmnKS.Inverse[AMKS];
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm;
(iMBL=AMKS.Inverse[invAnumu]/.{(a^2)+(-2+r) r->Del,1/(π (h+7 (1-h) (((-1+2 x2))^6)))->1/x2temp})//MatrixForm
(*Anumu//MatrixForm
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm
Inverse[invAnumu]//MatrixForm
Inverse[invAnumu].invAnumu//MatrixForm*)

R0=off;h=h1;
tadd=Table[Thread[{i1,i2,i3}={cv,45,21}];(*ext=ress[[i1,i2,i3]];*)fie=res2m[[i1,i2,i3]];Print["Coord (x1,x2): " ,xxx={x1tab[[i1]],x2tab[[i2]]},"; (r,θ): ", rtp={rtab[[i1]],θtab[[i2]]}];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.9};
Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];
ulo=invnew.uup/.sub;
Print["Lab (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}]];
(*MKS->FFFF*)
Thread[{u[0],u[1],u[2],u[3]}=uup];Thread[{uco[0],uco[1],uco[2],uco[3]}=ulo];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invnew.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invnew.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invnew.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invnew.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Norm (u^μ)(Subscript[u, μ])=",ulo.uup];
Print["(Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Co-moving (Subscript[B, i])=",coB=trxx.bup/.sub,";  0-th component should be 0"];
Print["Back to lab frame: (b^μ)=",trcoxx.coB];{rtp[[1]],Norm[coB],ρ},{cv,66,66,10}];

vec={1,0.2,-0.01,F};BFsol=FindInstance[0==vec.invgmn.vec/.sub,{B,F},Reals][[1]];
Print["Vector in BL: ",kBL=vec/.BFsol]
Print["Vector in MKS: ",lll=iMBL.kBL/.{Del->(a^2)+(-2+r) r,1/x2temp->1/(π (h+7 (1-h) (((-1+2 x2))^6)))}/.{R0->0.3,h->0.65}/.sub]

Print["Orthogonality in MKS to  ",lll.invnew.lll/.{R0->0.3,h->0.65}/.sub(*orthogonal to (10^(-5))*)]
Print["Orthogonality in BL to   ",kBL.invgmn.kBL/.sub(*orthogonal to (10^(-5))*)]

Clear[Del,kBL,R0];expr=iMBL.{kBL[0],kBL[1],kBL[2],kBL[3]}
expr==Table[Sum[iMBL[[m,j]]kBL[j-1],{j,4}],{m,4}]

New simulations

imp=Partition[Import["E:\\GRMHD\\new\\a0\\checkmri.dat"],128];
Dimensions[imp]
rtab=imp[[All,7,1]]
θtab=imp[[All,All,2]]

nn=180;rtab[[nn]]
ytab=θtab[[nn]]-(θtab[[nn,-1]]-θtab[[nn,1]])((Range[1,128])/127.5)
ListPlot[ytab,Joined->True]

dir="E:\\GRMHD\\new\\a0\\fieldlines\\";
read11[]:=(namefile=dir<>"fieldline9900.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)siz = FileByteCount[namefile]; xs = 11*256*128*64*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
res2m=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*64*128*256],11],256],128],{3,2,1,4}];
Close[namefile];);
read11[]
stt=StringToStream[strrr];
Do[nx1=Read[stt,Number],{i,2}];Do[nx2=Read[stt,Number],{i,1}];Do[nx3=Read[stt,Number],{i,1}];{nx1,nx2,nx3}
Do[x1min=Read[stt,Number],{i,1}];
Do[dx1=Read[stt,Number],{i,3}];
Do[off=Read[stt,Number],{i,6}];
Do[h1=Read[stt,Number],{i,3}];
Close[stt];{x1min,dx1,off,h1};
x1min=0.4220076586361168;
rtabx=off+Exp[x1tab=x1min+dx1 Range[0,nx1-1]];
(*np=nx2;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]]);*)
(*read22[]:=(namefilexx=dir<>"usdump0000";
fi22=OpenRead[namefilexx,BinaryFormat->True];
BinaryRead[fi22,"Integer32"];
tt=BinaryRead[fi22,"Real32"];
Do[Print[BinaryRead[fi22,"Integer32"]],{i,3}];
{x1min,x2min,x3min}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
{dx1,dx2,dx3}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
Do[Print[BinaryRead[fi22,"Real32"]],{i,8}];BinaryRead[fi22,"Real32"];
ress=Transpose[Partition[Partition[Partition[BinaryReadList[fi22,"Real32",30*32*64*256],30],256],64],{3,2,1,4}];);
read22[]*)

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];
gmn=Simplify[gmn];
invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1];
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm
Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}};
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}];
invnew=Inverse[AMKS].invgmnKS.Inverse[AMKS];
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm
(iMBL=AMKS.Inverse[invAnumu]/.{(a^2)+(-2+r) r->Del,1/(π (h+7 (1-h) (((-1+2 x2))^6)))->1/x2temp})//MatrixForm
(*Anumu//MatrixForm
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm
Inverse[invAnumu]//MatrixForm
Inverse[invAnumu].invAnumu//MatrixForm*)

Thread[{i1,i2,i3}={25,23,21}];(*ext=ress[[i1,i2,i3]];*)fie=res2m[[i1,i2,i3]];
uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]
sub={r->rtab[[i1]],θ->θtab[[i1,i2]],a->0.}
ulo=invgmnKS.uup/.sub
Print["Norm (u^μ)(Subscript[u, μ])=",ulo.uup]

R0=off;
tadd=Table[Thread[{i1,i2,i3}={cv,45,21}];(*ext=ress[[i1,i2,i3]];*)fie=res2m[[i1,i2,i3]];Print["Coord (x1,x2): " ,xxx={x1tab[[i1]],x2tab[[i2]]},"; (r,θ): ", rtp={rtab[[i1]],θtab[[i2]]}];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.9};
Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];
ulo=invnew.uup/.sub;
Print["Lab (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}]];
(*MKS->FFFF*)
Thread[{u[0],u[1],u[2],u[3]}=uup];Thread[{uco[0],uco[1],uco[2],uco[3]}=ulo];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invnew.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invnew.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invnew.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invnew.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Norm (u^μ)(Subscript[u, μ])=",ulo.uup];
Print["(Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Co-moving (Subscript[B, i])=",coB=trxx.bup/.sub,";  0-th component should be 0"];
Print["Back to lab frame: (b^μ)=",trcoxx.coB];{rtp[[1]],Norm[coB],ρ},{cv,45,45,10}];

vec={1,0.2,-0.01,F};BFsol=FindInstance[0==vec.invgmn.vec/.sub,{B,F},Reals][[1]];
Print["Vector in BL: ",kBL=vec/.BFsol]
Print["Vector in MKS: ",lll=iMBL.kBL/.{Del->(a^2)+(-2+r) r,1/x2temp->1/(π (h+7 (1-h) (((-1+2 x2))^6)))}/.{R0->0.3,h->0.65}/.sub]

Print["Orthogonality in MKS to  ",lll.invnew.lll/.{R0->0.3,h->0.65}/.sub(*orthogonal to (10^(-5))*)]
Print["Orthogonality in BL to   ",kBL.invgmn.kBL/.sub(*orthogonal to (10^(-5))*)]

Clear[Del,kBL,R0];expr=iMBL.{kBL[0],kBL[1],kBL[2],kBL[3]}
expr==Table[Sum[iMBL[[m,j]]kBL[j-1],{j,4}],{m,4}]

Shift of e1 by k due to Lorenz gauge

Transformation matrix to locally-flat co-moving frame from MKS

dir="G:\\GRMHD\\a07\\fieldlines\\";
read11[]:=(namefile=dir<>"fieldline0000.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
res2m=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*32*64*256],11],256],64],{3,2,1,4}];
Close[namefile];);
read22[]:=(namefilexx=dir<>"usdump0000";
fi22=OpenRead[namefilexx,BinaryFormat->True];
BinaryRead[fi22,"Integer32"];
tt=BinaryRead[fi22,"Real32"];
Do[Print[BinaryRead[fi22,"Integer32"]],{i,3}];
{x1min,x2min,x3min}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
{dx1,dx2,dx3}={BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"],BinaryRead[fi22,"Real32"]};
Do[Print[BinaryRead[fi22,"Real32"]],{i,8}];BinaryRead[fi22,"Real32"];
ress=Transpose[Partition[Partition[Partition[BinaryReadList[fi22,"Real32",30*32*64*256],30],256],64],{3,2,1,4}];);
read11[]
read22[]

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];
gmn=Simplify[gmn];
invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1];
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm
Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}};
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}];
invnew=Inverse[AMKS].invgmnKS.Inverse[AMKS];
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm;
(iMBL=AMKS.Inverse[invAnumu]/.{(a^2)+(-2+r) r->Del,1/(π (h+7 (1-h) (((-1+2 x2))^6)))->1/x2temp})//MatrixForm
(*Anumu//MatrixForm
(invAnumu=FullSimplify[gmn.Anumu.invgmnKS])//MatrixForm
Inverse[invAnumu]//MatrixForm
Inverse[invAnumu].invAnumu//MatrixForm*)

R0=0.3;h=0.65;
tadd=Table[Thread[{i1,i2,i3}={cv,34,21}];ext=ress[[i1,i2,i3]];fie=res2m[[i1,i2,i3]];Print["Coord (x1,x2,x3): " ,xxx=ext[[6;;8]],"; (r,θ,ϕ): ", rtp=ext[[9;;11]]];
sub={r->rtp[[1]],θ->rtp[[2]],ϕ->rtp[[3]],x2->xxx[[2]],a->0.7};
Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];
ulo=invnew.uup/.sub;
Print["Lab (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}]];
(*MKS->FFFF*)
Thread[{u[0],u[1],u[2],u[3]}=uup];Thread[{uco[0],uco[1],uco[2],uco[3]}=ulo];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invnew.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invnew.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invnew.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
(*xxlo[1]=invnew.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;*)
xlopr[1]=invnew.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Norm (u^μ)(Subscript[u, μ])=",InputForm[ulo.uup]];
Print["(Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Co-moving (Subscript[B, i])=",coB=trxx.bup/.sub,";  0-th component should be 0"];
Print["Back to lab frame: (b^μ)=",trcoxx.coB];{rtp[[1]],Norm[coB],ρ},{cv,175,175,10}];

sub
uup

Thread[{u[0],u[1],u[2],u[3]}=uup];Thread[{uco[0],uco[1],uco[2],uco[3]}=ulo];

xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invnew.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invnew.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invnew.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invnew.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invnew.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;

trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];


Print["Metrics in FLAT frame ",Chop[trxx.Inverse[invnew].Transpose[trxx]/.sub]//MatrixForm,"; should be diag(-1,1,1,1)"]
Print["Orthonormal condition ",Chop[trcoxx.trxx/.sub]//MatrixForm,"; should be diag(1,1,1,1)"]
Print["Transformation from FLAT to LAB frame ",trcoxx/.sub//MatrixForm]
Print["Transformation from LAB to FLAT frame ",trxx/.sub//MatrixForm]
Print["Inverse of transform from LAB to FLAT = transform from FLAT to LAB=",Inverse[trxx]//MatrixForm]

Print["Switching from (X^μ) to physical vector (Subscript[x, 0])=(Sqrt[-(Subscript[g, μ0ν])])(X^0); (Subscript[x, i])=(Sqrt[(Subscript[g, μν])])(X^i):  {(Subscript[x, 0,])(Subscript[x, i])}=",Xinv=Simplify[Partition[Flatten[{(Sqrt[-invnew/.a->0])[[1]],(Sqrt[invnew/.a->0])[[2;;4]]}],4]/.sub,{r>0,Sin[θ]>0}]]
Print["Transforming the physical vector to physical vector:",Re[Diagonal[trcoxx.Xinv/.sub]],"; should be (1,1,1,1)"]

Print["Velocity from FLAT to LAB frame ",trcoxx.{1,0,0,0},";  actual (u^μ)=",{u[0],u[1],u[2],u[3]}]
Print["Velocity from LAB to FLAT frame ",trxx.{u[0],u[1],u[2],u[3]},"; should be {1,0,0,0}"]

Vector product for e2 reconstruction

Definitions of frequency, magnetic fields, emissivities

ruru={rg->((2 6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2)),arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.26 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12))};
νH=((ee B)/(2π me cc ));νc=(3/2)Sin[θ]νH (γ^2);IjcSs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}];
QjcSs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(2/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,3000}];
VjcSs[θe_?NumberQ,νΩ_?NumberQ]:=(NIntegrate[((νΩ)/(γ^3))BesselK[(1/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}]+NIntegrate[((νΩ)/(γ^3))BesselK[(1/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,3000}]);
Iexpr=νH Sin[θ](((Sqrt[3]) (ee^2))/(2 cc)) IjcSs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];Qexpr=νH Sin[θ](((Sqrt[3]) (ee^2))/(2 cc)) QjcSs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];Vexpr=νH Cos[θ]((2(Sqrt[3]) (ee^2))/(3 cc)) VjcSs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];


quant={θ->1.,B->20,Te->3.8 (10^10),ν->1. (10^11),ne->1};AppendTo[quant,β->((ee B)/(me cc 2π ν))/.quant/.ruru];Sν=((2kb Te (ν^2))/(cc^2))/.ruru/.quant;
{IjS,QjS,VjS}={Iexpr,Qexpr,Vexpr}/.θe->((kb Te)/(me (cc^2)))/.ruru/.quant
{αIS,αQS,αVS}={IjS,QjS,VjS}/Sν(*absorptivities from emissivities *)
XX=θe(Sqrt[(Sqrt[2])Sin[θ]((10^3)((ee B)/(me cc 2π ν)))]);
ρQSs=(ω/(2cc))((4π ne (ee^2))/(me (ω^4)))(((((ee B)/(me cc))))^2)(((BesselK[1,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))+6θe)(fx=2.011 Exp[-((XX^1.035)/4.7)]-Cos[(XX/2)]Exp[-((XX^1.2)/2.73)]-0.011Exp[-(XX/47.2)])((Sin[θ])^2);
ρVSs=(*-*)(ω/cc)I (-I)((4π ne (ee^2))/(me (ω^3)))((ee B)/(me cc))((BesselK[0,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))(1-0.11 Log[1+0.035 XX])Cos[θ];
Thread[{ρQS,ρVS}={ρQSs,ρVSs}/.{ω->2π ν,θe->((kb Te)/(me (cc^2)))}/.ruru/.quant](*rotativities*)
fx/.θe->((kb Te)/(me (cc^2)))/.ruru/.quant

(ω/(2cc))((4π ne (ee^2))/(me (ω^4)))(((((ee B)/(me cc))))^2)/.{ω->2π nufr,ne->rho}
(ω/cc)I (-I)((4π ne (ee^2))/(me (ω^3)))((ee B)/(me cc))/.{ω->2π nufr,ne->rho}

kkk={0,Sin[θ],Cos[θ]};η=.;
xmat=(({{Cos[θ], 0, Sin[θ]}, {0, 1, 0}, {-Sin[θ], 0, Cos[θ]}}));mat=(({{1, 0, 0}, {0, Cos[θ], Sin[θ]}, {0, -Sin[θ], Cos[θ]}}));m90=(({{Cos[(π/2)], Sin[(π/2)], 0}, {-Sin[(π/2)], Cos[(π/2)], 0}, {0, 0, 1}}));kx={Sin[θ],0,Cos[θ]};(invmat=Simplify[Inverse[mat]])//MatrixForm;(invm90=Simplify[Inverse[m90]])//MatrixForm;Simplify[invmat.kkk];invm90.kx;
T11=(({{Cos[α β], η Sin[α β], 0}, {-η Sin[α β], Cos[α β], 0}, {0, 0, 1}}));
T22=(({{-((Sin[θ])^2)(((1-Cos[α β]))^2), η ((Sin[θ])^2)Sin[α β](1-Cos[α β]), η Sin[θ]Cos[θ]α β(1-Cos[α β])}, {-η ((Sin[θ])^2)Sin[α β](1-Cos[α β]), ((Sin[θ])^2)((Sin[α β])^2), Sin[θ]Cos[θ]α β Sin[α β]}, {-η Sin[θ]Cos[θ]α β(1-Cos[α β]), Sin[θ]Cos[θ]α β Sin[α β], ((Cos[θ])^2)(α^2)(β^2)}}));
(rot11=FullSimplify[invmat.T11.mat][[1;;2,1;;2]]/.α->ξ )//MatrixForm
(rot22=FullSimplify[invmat.T22.mat][[1;;2,1;;2]]/.α->ξ)//MatrixForm
η=-1;(*total=((I (ee^2)ne ω (ρ^2))/((cc^2)me BesselK[2,ρ]))(rot11 ((BesselK[2,Rex])/(Rex^2))-rot22 ((BesselK[3,Rex])/((β^2) (Rex^3))));*)(*final expression to integrate*)
Rex=((((ρ^2)-2I ρ ξ-(ξ^2)((Sin[θ])^2)+2((Sin[θ])^2)(1-Cos[β ξ])/(β^2)))^(1/2));
total=((I (ee^2)ne (ρ^2))/((cc^2)me BesselK[2,ρ]))(rot11 ((BesselK[2,Rex])/(Rex^2))-rot22 ((BesselK[3,Rex])/((β^2) (Rex^3))));
nonmagn=Normal[Series[total,{β,0,0}]];(*non-magnetic part of the response*)
magn=total-nonmagn;

acc=7;roro={ω->1,q->ee,ne->1,m->1,ρ->((me (cc^2))/(kb Te))}/.ruru/.quant;resx=Quiet[NIntegrate[(Cos[α]+I Sin[α])magn/.roro/.ruru/.quant/.ξ->(Cos[α]+I Sin[α])x/.α->(π/4.),{x,0,+∞},MaxRecursion->40,PrecisionGoal->acc]];
res=(cc/ν)resx/.ruru/.quant;
αIT=Im[res[[1,1]]+res[[2,2]]];αQT=Im[res[[1,1]]-res[[2,2]]];αVT=-2Re[res[[1,2]]];ρQT=Re[res[[1,1]]-res[[2,2]]];ρVT=2Im[res[[1,2]]];(*old definitions*)
αIT=Im[res[[1,1]]+res[[2,2]]];αQT=Im[res[[1,1]]-res[[2,2]]];αVT=2Re[res[[1,2]]];ρQT=Re[res[[2,2]]-res[[1,1]]];ρVT=2Im[res[[1,2]]];(*new definitions - all quantities are positive*)
{αIT,αQT,αVT,ρQT,ρVT}
{((αIS)/(αIT)),((αQS)/(αQT)),((αVS)/(αVT)),((ρQS)/(ρQT)),((ρVS)/(ρVT))}

correct inclusion of low-temperature regime

inum=100;θetab=Table[0.4 (1.05^i),{i,-inum,inum}];
lookuprV=Log[Table[((BesselK[0,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))/.θe->θetab[[i]],{i,Length[θetab]}]];
lookuprQ=Log[Table[((BesselK[1,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))+6θe/.θe->θetab[[i]],{i,Length[θetab]}]];

Export["E:\\mathematica\\3DMHD\\lookuprQa.dat",lookuprQ]
Export["E:\\mathematica\\3DMHD\\lookuprVa.dat",lookuprV]

inum=100;θetab=Table[0.4 (1.05^i),{i,-inum,inum}]6 (10^9)

Length[θetab]-(2inum+1)

{((BesselK[0,((θe)^(-1))])/(BesselK[2,((θe)^(-1))])),((BesselK[1,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))+6θe}/.θe->Exp[0.12575]

r1=20030;B1={-7.64733*(10^(-5)),0.000137425,0.0050052867};
r2=2891.2;B2={-0.0004173620,0.0009503498,0.034529621};
(r1/r2)(B1/B2)(*for 2.0 extension*)

Definitions of cos2k, sin2k

see Lei, Roman 2010 paper

a - North, b - East

Evalpointzero functions

Checking gdump/fieldline, making dxdxp

for old computation

(9+16+16+1+4+16+8)

nsize=256 64 (9+16+16+1+4+16)8
issize=256 64 (9+16+16+1+4+16+7)8
actsize=9175208;
actsize-issize

dir = NotebookDirectory[];
sp=1;
ndd=650;ncuttab={145,147,149,153,162};
dirin="E:\\GRMHD\\a";anum={0.,0.5,0.7,0.9,0.98};atab={"0","05","07","09","098"};
readgd[sp_] := (namefile ="E:\\GRMHD\\a"<>atab[[sp]]<>"\\usgdump2d";
   siz = FileByteCount[namefile]; xs = 256*64*70*8;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", 168]];
resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],70],256],{2,1,3}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[namefile];);

dir = NotebookDirectory[];
read11[sp_,fnum_] := (namefile ="E:\\GRMHD\\a"<>atab[[sp]]<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];
   Close[namefile];);
readgd[sp]
read11[sp,9950]
Dimensions[resgd]

resgd[[All,1,8]]

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm;
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm;
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
(*G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];*)
gmn=Simplify[gmn];
(invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1])//MatrixForm;
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;
(Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}});
(invAnumu=Transpose[Inverse[Anumu]])//MatrixForm;



i1=4;i2=2;i3=21;
rtp=resgd[[i1,i2,8;;10]];
xxx=resgd[[i1,i2,5;;7]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->anum[[sp]]};

gcon=Partition[resgd[[i1,i2,11;;26]],4];gcov=Partition[resgd[[i1,i2,27;;42]],4];
gdet=resgd[[i1,i2,43]];dxdxp=Partition[resgd[[i1,i2,48;;63]],4];
fie=res2m[[i1,i2,i3]];Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];ulo=gcov.uup;
Print["MKS (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}],"; (Subscript[(b^μ), KS])=",bupKS=dxdxp.bup];
Print["In KS coordinates: (Subscript[(u^μ), KS])=", uupKS=dxdxp.uup,";  (Subscript[u, μKS])=",uloKS=invgmnKS.uupKS/.sub];
Print["In BL coordinates: (Subscript[(u^μ), BL])=", uupBL=invAnumu.uupKS/.sub,";  (Subscript[u, μBL])=",uloBL=invgmn.uupBL/.sub];
Print["Test 1: (g^(μν))(Subscript[g, νσ])==(Subscript[(δ^μ), σ]); should be identity matrix:   ",Chop[gcov.gcon//MatrixForm]]
Print["Test 2: gdet==(Sqrt[-Det[gcov]]);  should be unity:      ",((Sqrt[-Det[gcov]])/gdet)]
Print["Test 3: normalization of velocity in MKS: (u^μ)(Subscript[u, μ])=-1:  ",ulo.uup]
Print["Test 4: normalization of velocity in KS: (Subscript[(u^μ), KS])(Subscript[u, (Subscript[μ, KS])])=-1:  ",uupKS.uloKS]
Print["Test 5: normalization of velocity in BL: (Subscript[(u^μ), BL])(Subscript[u, (Subscript[μ, BL])])=-1:  ",uupBL.uloBL]

Thread[{u[0],u[1],u[2],u[3]}=uupKS];Thread[{uco[0],uco[1],uco[2],uco[3]}=uloKS];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invgmnKS.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invgmnKS.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invgmnKS.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invgmnKS.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Test 6: (Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Test 7: co-moving (Subscript[B, i])=",coB=trxx.bupKS/.sub,";  0-th component should be 0"];
Print["Back to KS frame: (Subscript[(b^μ), KS])=",trcoxx.coB];

i1=34;i2=44;i3=21;
rtp=resgd[[i1,i2,8;;10]];
xxx=resgd[[i1,i2,5;;7]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.0};
DiagonalMatrix[{1,r-R0,(π(h+7(1-h)(((2x2-1))^6))),π}]/.{r->resgd[[i1,i2,8]],x2->resgd[[i1,i2,6]]}/.{R0->0.3,h->0.65}

np=64;h=.;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]])/.{h->0.65};
ncut=ncuttab[[sp]];
(*rtab11=resgd[[1;;ncut,1,8]];rcut=rtab11[[-1]];*)
rtab11=Transpose[Import["E:\\GRMHD\\a05\\Tsmap05.dat"]][[1,1;;ncut]];rcut=rtab11[[-1]];

rrmax=340000.;
rtab22=Table[rcut*(((rrmax/rcut))^((k-ncut+1)/(ndd-ncut))),{k,ncut,ndd-1}];
rtab=Join[rtab11,rtab22];rn=Length[rtab];tn=Total[x2tab];
ListLogPlot[rtab];
i1=.;i2=.;
dxdxptab=Table[DiagonalMatrix[{1.,r-R0,(π(h+7(1-h)(((2x2-1))^6))),1.π}]/.{r->rtab[[i1]],x2->x2tab[[i2]]}/.{R0->0.3,h->0.65},{i1,Length[rtab]},{i2,Length[x2tab]}];
coord=Outer[List,Log[rtab],-Cos[θtab]];
fif=OpenWrite["E:\\GRMHD\\a"<>atab[[sp]]<>"\\dxdxp.dat",BinaryFormat->True];
BinaryWrite[fif,Flatten[coord],"Real32"]
BinaryWrite[fif,Flatten[dxdxptab],"Real32"]
Close[fif];

Dimensions[grids]

grids[[-1,2]]

Dimensions[dxdxptab]
Total[Flatten[dxdxptab]]

for new computations

(9+16+16+1+4+16+7)

nsize=256 128 (9+16+16+1+4+16)8
issize=256 128 (9+16+16+1+4+16+7)8
actsize=18088096;
actsize-issize

xs = 256*128*69*8

256 128 69 8

dir = NotebookDirectory[];
sp=5;
ndd=650;ncuttab={100,110,116,128,138};
dirin="E:\\GRMHD\\a";anum={0.,0.5,0.7,0.9,0.98};atab={"0","05","07","09","098"};
readgd[sp_] := (namefile ="E:\\GRMHD\\new\\a"<>atab[[sp]]<>"\\usgdump2d";
   siz = FileByteCount[namefile]; xs = 256*128*69*8;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],69],256],{2,1,3}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[namefile];);

read11[sp_,fnum_] := (namefile ="E:\\GRMHD\\new\\a"<>atab[[sp]]<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*256*128*64*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];
   Close[namefile];);
readgd[sp]
read11[sp,2571]
Dimensions[resgd]
(*extension to large radius: 
quantities extended as power-laws; 
elements of dxdxp matrices are extended as power-laws; 
grids are extended as power-laws; theta grid is fixed at 40M boundary;*)

avtab=Import["E:\\GRMHD\\new\\a07\\Tsmap07.dat"];
lpl1=ListLogLogPlot[Transpose[Delete[Transpose[avtab],2]][[1;;160]],Frame->True];
lpl2=ListLogLogPlot[avtab[[1;;160,1;;2]],Frame->True];
pl1=LogLogPlot[0.047((((r/10)))^(-0.3)),{r,avtab[[1,1]],avtab[[160,1]]},PlotStyle->Red];
Show[pl1,lpl1]

avtab[[1,1]]
avtab[[100,1]]
avtab[[160,1]]
newpar={2.02245567, 4.23384276 (10^11),0.0402070586}
oldpar={1.99835404, 1.18096645 (10^12), 5.41663751}
(*thus, need to increase rhonor by x100, increase heat to 0.55 *)

Dimensions[res2m]
res2m[[1,1,1]]
res2m[[1,2,1]]

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm;
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm;
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
(*G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];*)
gmn=Simplify[gmn];
(invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1])//MatrixForm;
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;
(Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}});
(invAnumu=Transpose[Inverse[Anumu]])//MatrixForm;


i1=54;i2=87;i3=21;
rtp=resgd[[i1,i2,8;;10]];
xxx=resgd[[i1,i2,5;;7]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.98};

gcon=Partition[resgd[[i1,i2,11;;26]],4];gcov=Partition[resgd[[i1,i2,27;;42]],4];
gdet=resgd[[i1,i2,43]];dxdxp=Partition[resgd[[i1,i2,48;;63]],4];
fie=res2m[[i1,i2,i3]];Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];ulo=gcov.uup;
Print["MKS (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}],"; (Subscript[(b^μ), KS])=",bupKS=dxdxp.bup];
Print["In KS coordinates: (Subscript[(u^μ), KS])=", uupKS=dxdxp.uup,";  (Subscript[u, μKS])=",uloKS=invgmnKS.uupKS/.sub];
Print["In BL coordinates: (Subscript[(u^μ), BL])=", uupBL=invAnumu.uupKS/.sub,";  (Subscript[u, μBL])=",uloBL=invgmn.uupBL/.sub];
Print["Back to KS coordinates: (Subscript[(u^μ), KS])=", Transpose[Anumu].uupBL/.sub];
Print["Test 1: (g^(μν))(Subscript[g, νσ])==(Subscript[(δ^μ), σ]); should be identity matrix:   ",Chop[gcov.gcon//MatrixForm]]
Print["Test 2: gdet==(Sqrt[-Det[gcov]]);  should be unity:      ",((Sqrt[-Det[gcov]])/gdet)]
Print["Test 3: normalization of velocity in MKS: (u^μ)(Subscript[u, μ])=-1:  ",ulo.uup]
Print["Test 4: normalization of velocity in KS: (Subscript[(u^μ), KS])(Subscript[u, (Subscript[μ, KS])])=-1:  ",uupKS.uloKS]
Print["Test 5: normalization of velocity in BL: (Subscript[(u^μ), BL])(Subscript[u, (Subscript[μ, BL])])=-1:  ",uupBL.uloBL]

Thread[{u[0],u[1],u[2],u[3]}=uupKS];Thread[{uco[0],uco[1],uco[2],uco[3]}=uloKS];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invgmnKS.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invgmnKS.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invgmnKS.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invgmnKS.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Test 6: (Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Test 7: co-moving (Subscript[B, i])=",coB=trxx.bupKS/.sub,";  0-th component should be 0"];
Print["Back to KS frame: (Subscript[(b^μ), KS])=",trcoxx.coB];

(*np=64;h=.;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]])/.{h->0.65};*)
ncut=ncuttab[[sp]];
rtab11=resgd[[1;;ncut,1,8]];rcut=rtab11[[-1]];
(*rtab11=Transpose[Import["E:\\GRMHD\\a05\\Tsmap05.dat"]][[1,1;;ncut]];rcut=rtab11[[-1]];*)
rrmax=340000.;
rtab22=Table[rcut*(((rrmax/rcut))^((k-ncut+1)/(ndd-ncut))),{k,ncut,ndd-1}];
rtab=Join[rtab11,rtab22];
(*rn=Length[rtab];tn=Total[x2tab];ListLogPlot[rtab];i1=.;i2=.;*)
(*dxdxptab=Table[DiagonalMatrix[{1.,r-R0,(π(h+7(1-h)(((2x2-1))^6))),1.π}]/.{r->rtab[[i1]],x2->x2tab[[i2]]}/.{R0->0.3,h->0.65},{i1,Length[rtab]},{i2,Length[x2tab]}];*)

(*coord=Outer[List,Log[rtab],-Cos[θtab]];*)

rcmax=resgd[[256,1,8]];
coordx11=Table[Transpose[{Log[resgd[[num,All,8]]],-Cos[resgd[[num,All,9]]]}],{num,ncut}];
dxdxptab11=Table[Partition[resgd[[i1,i2,48;;63]],4],{i1,ncut},{i2,128}];
dxdxptab22=Table[{{1.,0.,0.,0.},{0.,rtab[[i1+1]],0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}},{i1,ncut,ndd-1},{i2,128}];
thtab22=-Cos[resgd[[ncut,All,9]]];
coordx22=Outer[List,Log[rtab22],thtab22];
Dimensions[coordfin=Join[coordx11,coordx22]]
Dimensions[dxdxpfin=Join[dxdxptab11,dxdxptab22]]
fif=OpenWrite["E:\\GRMHD\\new\\a"<>atab[[sp]]<>"\\dxdxp.dat",BinaryFormat->True];
BinaryWrite[fif,Flatten[coordfin],"Real32"]
BinaryWrite[fif,Flatten[dxdxpfin],"Real32"]
Close[fif];
(*fullcoord=Table[Transpose[{Log[resgd[[num,All,8]]],-Cos[resgd[[num,All,9]]]}],{num,256}];*)
(*fullcoord=Join[fullcoord,Abs[Reverse[fullcoord,2]],2];*)

tn=23;
dxdxpfin[[ncut-1,tn]]
dxdxpfin[[ncut,tn]]
dxdxpfin[[ncut+1,tn]]
coordfin[[ncut-1,tn]]
coordfin[[ncut,tn]]
coordfin[[ncut+1,tn]]

tn=64;
{lpl22=ListLogLogPlot[Transpose[{rtab11,dxdxptab11[[All,tn,2,2]]}]],ListLogPlot[dxdxptab11[[All,9,3,3]]],ListLogPlot[Abs[dxdxptab11[[All,tn,3,2]]]],ListLogPlot[dxdxptab11[[All,tn,4,4]]]}
pl22=LogLogPlot[r,{r,rtab11[[1]],rtab11[[-1]]},PlotStyle->Red];
Show[lpl22,pl22]

dxdxp//MatrixForm

ListPlot[dxdxptab11[[100,All,3,3]]]
ListPlot[dxdxptab11[[100,All,3,2]]]

(*(0,0) -> 1; (1,1) -> r; (2,2) -> π/2; (3,3) -> 2π; (3,2) -> 0*)
{{1.,0.,0.,0.},{0.,r,0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}}

gcov
invgmnKS/.sub
Transpose[dxdxp].invgmnKS.dxdxp/.sub

(*given u, iKS, dxdxp;
1. Find iMKS=Transpose[dxdxp].iKS.dxdxp;
2. Find ulo=iMKS.u;*)

for ADAF computations

(9+16+16+1+4+16+7)

256*128*69*8+156
264 126 69 8+160
87817358-11*264*126*60*4

dir = NotebookDirectory[];
sp=4;
ndd=650;ncuttab={110,110,126,134,138};
dirin="E:\\GRMHD\\a";anum={0.,0.5,0.7,0.9,0.98};atab={"0","05","07","09","098"};
readgd[sp_] := (namefile ="E:\\GRMHD\\ADAFs\\a"<>atab[[sp]]<>"\\usgdump2d";
   siz = FileByteCount[namefile]; xs = 264*126*69*8;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],69],264],{2,1,3}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[namefile];);

read11[sp_,fnum_] := (namefile ="E:\\GRMHD\\ADAFs\\a"<>atab[[sp]]<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*264*126*60*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*60*126*264], 11], 264], 126], {3, 2, 1, 4}];
   Close[namefile];);
readgd[sp]
read11[sp,3149]
Dimensions[resgd]
(*extension to large radius: 
quantities extended as power-laws; 
elements of dxdxp matrices are extended as power-laws; 
grids are extended as power-laws; theta grid is fixed at 40M boundary;*)

Dimensions[resgd]
resgd[[1,2]]

avtab=Import["E:\\GRMHD\\ADAFs\\a0\\Tsmap0.dat"];
lpl1=ListLogLogPlot[Transpose[Delete[Transpose[avtab],2]][[1;;160]],Frame->True];
lpl2=ListLogLogPlot[avtab[[1;;160,1;;2]],Frame->True];
pl1=LogLogPlot[0.047((((r/10)))^(-1.1)),{r,avtab[[1,1]],avtab[[160,1]]},PlotStyle->Red];
Show[pl1,lpl1]

avtab[[1,1]]
avtab[[100,1]]
avtab[[160,1]]
newpar={2.02245567, 4.23384276 (10^11),0.0402070586}
oldpar={1.99835404, 1.18096645 (10^12), 5.41663751}
(*thus, need to increase rhonor by x100, increase heat to 0.55 *)

Dimensions[res2m]
res2m[[1,1,1]]
res2m[[1,2,1]]

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm;
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm;
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
(*G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];*)
gmn=Simplify[gmn];
(invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1])//MatrixForm;
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;
(Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}});
(invAnumu=Transpose[Inverse[Anumu]])//MatrixForm;


i1=24;i2=47;i3=21;
rtp=resgd[[i1,i2,8;;10]];
xxx=resgd[[i1,i2,5;;7]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.9};

gcon=Partition[resgd[[i1,i2,11;;26]],4];gcov=Partition[resgd[[i1,i2,27;;42]],4];
gdet=resgd[[i1,i2,43]];dxdxp=Partition[resgd[[i1,i2,48;;63]],4];
fie=res2m[[i1,i2,i3]];Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];ulo=gcov.uup;
Print["MKS (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}],"; (Subscript[(b^μ), KS])=",bupKS=dxdxp.bup];
Print["In KS coordinates: (Subscript[(u^μ), KS])=", uupKS=dxdxp.uup,";  (Subscript[u, μKS])=",uloKS=invgmnKS.uupKS/.sub];
Print["In BL coordinates: (Subscript[(u^μ), BL])=", uupBL=invAnumu.uupKS/.sub,";  (Subscript[u, μBL])=",uloBL=invgmn.uupBL/.sub];
Print["Back to KS coordinates: (Subscript[(u^μ), KS])=", Transpose[Anumu].uupBL/.sub];
Print["Test 1: (g^(μν))(Subscript[g, νσ])==(Subscript[(δ^μ), σ]); should be identity matrix:   ",Chop[gcov.gcon//MatrixForm]]
Print["Test 2: gdet==(Sqrt[-Det[gcov]]);  should be unity:      ",((Sqrt[-Det[gcov]])/gdet)]
Print["Test 3: normalization of velocity in MKS: (u^μ)(Subscript[u, μ])=-1:  ",ulo.uup]
Print["Test 4: normalization of velocity in KS: (Subscript[(u^μ), KS])(Subscript[u, (Subscript[μ, KS])])=-1:  ",uupKS.uloKS]
Print["Test 5: normalization of velocity in BL: (Subscript[(u^μ), BL])(Subscript[u, (Subscript[μ, BL])])=-1:  ",uupBL.uloBL]

Thread[{u[0],u[1],u[2],u[3]}=uupKS];Thread[{uco[0],uco[1],uco[2],uco[3]}=uloKS];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invgmnKS.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invgmnKS.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invgmnKS.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invgmnKS.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Test 6: (Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Test 7: co-moving (Subscript[B, i])=",coB=trxx.bupKS/.sub,";  0-th component should be 0"];
Print["Back to KS frame: (Subscript[(b^μ), KS])=",trcoxx.coB];

(*np=64;h=.;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]])/.{h->0.65};*)
ncut=ncuttab[[sp]];
rtab11=resgd[[1;;ncut,1,8]];rcut=rtab11[[-1]];
(*rtab11=Transpose[Import["E:\\GRMHD\\a05\\Tsmap05.dat"]][[1,1;;ncut]];rcut=rtab11[[-1]];*)
rrmax=340000.;
rtab22=Table[rcut*(((rrmax/rcut))^((k-ncut+1)/(ndd-ncut))),{k,ncut,ndd-1}];
rtab=Join[rtab11,rtab22];
(*rn=Length[rtab];tn=Total[x2tab];ListLogPlot[rtab];i1=.;i2=.;*)
(*dxdxptab=Table[DiagonalMatrix[{1.,r-R0,(π(h+7(1-h)(((2x2-1))^6))),1.π}]/.{r->rtab[[i1]],x2->x2tab[[i2]]}/.{R0->0.3,h->0.65},{i1,Length[rtab]},{i2,Length[x2tab]}];*)

(*coord=Outer[List,Log[rtab],-Cos[θtab]];*)

rcmax=resgd[[264,1,8]];
coordx11=Table[Transpose[{Log[resgd[[num,All,8]]],-Cos[resgd[[num,All,9]]]}],{num,ncut}];
dxdxptab11=Table[Partition[resgd[[i1,i2,48;;63]],4],{i1,ncut},{i2,126}];
dxdxptab22=Table[{{1.,0.,0.,0.},{0.,rtab[[i1+1]],0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}},{i1,ncut,ndd-1},{i2,126}];
thtab22=-Cos[resgd[[ncut,All,9]]];
coordx22=Outer[List,Log[rtab22],thtab22];
Dimensions[coordfin=Join[coordx11,coordx22]]
Dimensions[dxdxpfin=Join[dxdxptab11,dxdxptab22]]
fif=OpenWrite["E:\\GRMHD\\ADAFs\\a"<>atab[[sp]]<>"\\dxdxp.dat",BinaryFormat->True];
BinaryWrite[fif,Flatten[coordfin],"Real32"]
BinaryWrite[fif,Flatten[dxdxpfin],"Real32"]
Close[fif];
(*fullcoord=Table[Transpose[{Log[resgd[[num,All,8]]],-Cos[resgd[[num,All,9]]]}],{num,256}];*)
(*fullcoord=Join[fullcoord,Abs[Reverse[fullcoord,2]],2];*)

tn=23;
dxdxpfin[[ncut-1,tn]]
dxdxpfin[[ncut,tn]]
dxdxpfin[[ncut+1,tn]]
coordfin[[ncut-1,tn]]
coordfin[[ncut,tn]]
coordfin[[ncut+1,tn]]

tn=64;
{lpl22=ListLogLogPlot[Transpose[{rtab11,dxdxptab11[[All,tn,2,2]]}]],ListLogPlot[dxdxptab11[[All,9,3,3]]],ListLogPlot[Abs[dxdxptab11[[All,tn,3,2]]]],ListLogPlot[dxdxptab11[[All,tn,4,4]]]}
pl22=LogLogPlot[r,{r,rtab11[[1]],rtab11[[-1]]},PlotStyle->Red];
Show[lpl22,pl22]

dxdxp//MatrixForm

ListPlot[dxdxptab11[[100,All,3,3]]]
ListPlot[dxdxptab11[[100,All,3,2]]]

(*(0,0) -> 1; (1,1) -> r; (2,2) -> π/2; (3,3) -> 2π; (3,2) -> 0*)
{{1.,0.,0.,0.},{0.,r,0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}}

gcov
invgmnKS/.sub
Transpose[dxdxp].invgmnKS.dxdxp/.sub

(*given u, iKS, dxdxp;
1. Find iMKS=Transpose[dxdxp].iKS.dxdxp;
2. Find ulo=iMKS.u;*)

for Jon’s jet simulations with QPOs (i.e. thickdisk7)

dir = NotebookDirectory[];
ndd=650;ncut=120;NDIM=4;
dirin="E:\\GRMHD\\Jon\\";anum=0.9375;
readgd[spstr_] := (namefile =dirin<>atab<>"\\gdump.bin";
   siz = FileByteCount[namefile]; xs = 272*128*256*8*126;
   fi11 = OpenRead[namefile, BinaryFormat -> True];strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
SetStreamPosition[fi11,siz-xs];
(*resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],(NDIM^2)],272],{2,1,3}];*)
resgd = Transpose[Partition[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],126],272],128],{3,2,1,4}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[fi11];);
readgd2d[spstr_] := (namefile =dirin<>atab<>"\\gdump.bin";
   siz = FileByteCount[namefile]; phlen=256;xs = 272*128*phlen*8*126;
   fi11 = OpenRead[namefile, BinaryFormat -> True];strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
SetStreamPosition[fi11,siz-xs];
(*resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],(NDIM^2)],272],{2,1,3}];*)
resgd2d = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8/phlen],126],272],{2,1,3}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[fi11];);

read11[spstr_,fnum_] := (namefile =dirin<>spstr<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*272 128 256*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
str=StringToStream[strrr];param=ReadList[str,Number];Close[str];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", xs/4], 11], 272], 128], {3, 2, 1, 4}];
   Close[namefile];);
atab="thickdisk7";
(*readgd[atab]*)
readgd2d[atab]
read11[atab,7034]
Dimensions[res2m]
Dimensions[resgd]
Dimensions[resgd2d]
(*extension to large radius: 
quantities extended as power-laws; 
elements of dxdxp matrices are extended as power-laws; 
grids are extended as power-laws; theta grid is fixed at 40M boundary;*)

(*rtab=resgd[[All,1,1,7]];
thtab=resgd[[1,All,1,8]];
phtab=resgd[[1,1,All,9]];*)
rtab=resgd2d[[All,1,7]];
thtab=resgd2d[[1,All,8]];
{rdim=param[[2]],thdim=param[[3]],phdim=param[[4]],anum=param[[13]]}
ncut=143;ncut=rdim;
rtab[[ncut]]
(*dxdxptab11=Partition[resgd[[All,All,1,-(NDIM^2);;-1]],{272,128,4}];*)
(*dxdxptab11=Table[Partition[resgd[[i1,i2,1,-(NDIM^2);;-1]],4],{i1,ncut},{i2,128}];*)
dxdxptab11=Table[Partition[resgd2d[[i1,i2,-(NDIM^2);;-1]],4],{i1,ncut},{i2,128}];
rtab11=rtab[[1;;ncut]];
Dimensions[dxdxptab11]

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm;
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm;
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
(*G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];*)
gmn=Simplify[gmn];
(invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1])//MatrixForm;
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;
(Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}});
(invAnumu=Transpose[Inverse[Anumu]])//MatrixForm;

i1=24;i2=47;i3=21;
(*rtp=resgd[[i1,i2,i3,7;;9]];
xxx=resgd[[i1,i2,i3,4;;6]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.9375};

gcon=Partition[resgd[[i1,i2,i3,74;;89]],4];
gcov=Partition[resgd[[i1,i2,i3,90;;105]],4];
gdet=resgd[[i1,i2,i3,106]];
dxdxp=Partition[resgd[[i1,i2,i3,111;;126]],4];*)
rtp=resgd2d[[i1,i2,7;;9]];
xxx=resgd2d[[i1,i2,4;;6]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.9375};

gcon=Partition[resgd2d[[i1,i2,74;;89]],4];
gcov=Partition[resgd2d[[i1,i2,90;;105]],4];
gdet=resgd2d[[i1,i2,106]];
dxdxp=Partition[resgd2d[[i1,i2,111;;126]],4];
fie=res2m[[i1,i2,i3]];Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];ulo=gcov.uup;
Print["MKS (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}],"; (Subscript[(b^μ), KS])=",bupKS=dxdxp.bup];
Print["In KS coordinates: (Subscript[(u^μ), KS])=", uupKS=dxdxp.uup,";  (Subscript[u, μKS])=",uloKS=invgmnKS.uupKS/.sub];
Print["In BL coordinates: (Subscript[(u^μ), BL])=", uupBL=invAnumu.uupKS/.sub,";  (Subscript[u, μBL])=",uloBL=invgmn.uupBL/.sub];
Print["Back to KS coordinates: (Subscript[(u^μ), KS])=", Transpose[Anumu].uupBL/.sub];
Print["Test 1: (g^(μν))(Subscript[g, νσ])==(Subscript[(δ^μ), σ]); should be identity matrix:   ",Chop[gcov.gcon//MatrixForm]]
Print["Test 2: gdet==(Sqrt[-Det[gcov]]);  should be unity:      ",((Sqrt[-Det[gcov]])/gdet)]
Print["Test 3: normalization of velocity in MKS: (u^μ)(Subscript[u, μ])=-1:  ",ulo.uup]
Print["Test 4: normalization of velocity in KS: (Subscript[(u^μ), KS])(Subscript[u, (Subscript[μ, KS])])=-1:  ",uupKS.uloKS]
Print["Test 5: normalization of velocity in BL: (Subscript[(u^μ), BL])(Subscript[u, (Subscript[μ, BL])])=-1:  ",uupBL.uloBL]

Thread[{u[0],u[1],u[2],u[3]}=uupKS];Thread[{uco[0],uco[1],uco[2],uco[3]}=uloKS];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invgmnKS.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invgmnKS.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invgmnKS.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invgmnKS.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Test 6: (Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Test 7: co-moving (Subscript[B, i])=",coB=trxx.bupKS/.sub,";  0-th component should be 0"];
Print["Back to KS frame: (Subscript[(b^μ), KS])=",trcoxx.coB];

uupKS
uloKS
invgmnKS.{1,0,0,0}
dxdxp.{1,0,0,0}
invgmnKS

rtab11=resgd[[1;;ncut,1,1,7]];rcut=rtab11[[-1]]
(*rtab11=Transpose[Import["E:\\GRMHD\\a05\\Tsmap05.dat"]][[1,1;;ncut]];rcut=rtab11[[-1]];*)
rrmax=340000.;
rtab22=Table[rcut*(((rrmax/rcut))^((k-ncut+1)/(ndd-ncut))),{k,ncut,ndd-1}];
rtab=Join[rtab11,rtab22];
(*dxdxptab=Table[DiagonalMatrix[{1.,r-R0,(π(h+7(1-h)(((2x2-1))^6))),1.π}]/.{r->rtab[[i1]],x2->x2tab[[i2]]}/.{R0->0.3,h->0.65},{i1,Length[rtab]},{i2,Length[x2tab]}];*)
rcmax=resgd[[-1,1,1,7]];

(*np=64;h=.;θtab=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.x2->(x2tab=N[Range[(1/(2np)),1-(1/(2np)),(1/np)]])/.{h->0.65};*)
(*coordx11=Table[Transpose[{Log[resgd[[num,All,1,7]]],-Cos[resgd[[num,All,1,8]]]}],{num,ncut}];
dxdxptab11=Table[Partition[resgd[[i1,i2,1,111;;126]],4],{i1,ncut},{i2,thdim}];
dxdxptab22=Table[{{1.,0.,0.,0.},{0.,rtab[[i1+1]],0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}},{i1,ncut,ndd-1},{i2,thdim}];
thtab22=-Cos[resgd[[ncut,All,1,8]]];*)
coordx11=Table[Transpose[{Log[resgd2d[[num,All,7]]],-Cos[resgd2d[[num,All,8]]]}],{num,ncut}];
dxdxptab11=Table[Partition[resgd2d[[i1,i2,111;;126]],4],{i1,ncut},{i2,thdim}];
dxdxptab22=Table[{{1.,0.,0.,0.},{0.,rtab[[i1+1]],0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}},{i1,ncut,ndd-1},{i2,thdim}];
thtab22=-Cos[resgd2d[[ncut,All,8]]];

coordx22=Outer[List,Log[rtab22],thtab22];
Dimensions[coordfin=Join[coordx11,coordx22]]
Dimensions[dxdxpfin=Join[dxdxptab11,dxdxptab22]]
fif=OpenWrite["E:\\GRMHD\\Jon\\"<>atab<>"\\dxdxp.dat",BinaryFormat->True];
BinaryWrite[fif,Flatten[coordfin],"Real32"]
BinaryWrite[fif,Flatten[dxdxpfin],"Real32"]
Close[fif];
(*fullcoord=Table[Transpose[{Log[resgd[[num,All,8]]],-Cos[resgd[[num,All,9]]]}],{num,256}];*)
(*fullcoord=Join[fullcoord,Abs[Reverse[fullcoord,2]],2];*)

tn=64;
{lpl22=ListLogLogPlot[Transpose[{rtab11,dxdxptab11[[All,tn,2,2]]}]],ListLogPlot[dxdxptab11[[All,9,3,3]]],ListLogPlot[Abs[dxdxptab11[[All,tn,3,2]]]],ListLogPlot[dxdxptab11[[All,tn,4,4]]]}
pl22=LogLogPlot[r,{r,rtab11[[1]],rtab11[[-1]]},PlotStyle->Red];
Show[lpl22,pl22]

dxdxptab11[[ncut,50]]//MatrixForm

ncut=143;
GraphicsRow[{ListPlot[dxdxptab11[[ncut,All,3,3]]],ListPlot[dxdxptab11[[ncut-10,All,3,3]],ImageSize->400]}](*th-th*)
GraphicsRow[{ListPlot[dxdxptab11[[ncut,All,3,2]]],ListPlot[dxdxptab11[[ncut-10,All,3,2]]]}](*r-th*)

(*(0,0) -> 1; (1,1) -> r; (2,2) -> π/2; (3,3) -> 2π; (3,2) -> 0*)
{{1.,0.,0.,0.},{0.,r,0.,0.},{0.,0.,(π/2.),0.},{0.,0.,0.,2.π}}

11*272 128 256*4
382976 1024
8984199561
NDIM=4;
272 128 256*8(3+3+3+(NDIM^3)+(NDIM^2)+(NDIM^2)+1+NDIM+(NDIM^2))
3+3+3+(NDIM^3)+(NDIM^2)+(NDIM^2)+1+NDIM+(NDIM^2)

testing for zeros in Jon’s simulations

dir = NotebookDirectory[];
ndd=650;ncut=120;NDIM=4;
dirin="E:\\GRMHD\\Jon\\";anum=0.9375;
read11[spstr_,fnum_] := (namefile =dirin<>spstr<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*272 128 256*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
str=StringToStream[strrr];param=ReadList[str,Number];Close[str];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", xs/4], 11], 272], 128], {3, 2, 1, 4}];
   Close[namefile];);
atab="thickdisk7";
Dimensions[res2m]

arr={3786,3790,3791,3888,3891,3892,3893,3894,3896,5000,5001,5002};
ftest[fnum_]:=(read11[atab,fnum];
Print["fnum=",fnum];Print[Union[Position[res2m[[50]],0.][[All,2]]]];
Print[Union[Position[res2m[[1,All,All,1]],0.][[All,2]]]];)

ftest[#]&/@arr

ftest[3893]

Length[p1=Position[res2m[[-1,All,All,1]],0.]]
Length[p2=Position[res2m[[-10,All,All,1]],0.]]
Length[p3=Position[res2m[[1,All,All,1]],0.]]
Length[p4=Position[res2m[[130,All,All,1]],0.]]

p1

Complement[p2,p1]
Complement[p4,p1]

Poloidal grid

multilinear inverse lookup and interpolation

grid=resgd[[All,All,8;;10]];
rt=5.6;θt=2.3;ϕt=1.2;rlen=Dimensions[grid][[1]];thlen=Dimensions[grid][[2]];phlen=Dimensions[res2m][[3]];rmin=grid[[1,1,1]];rmax=grid[[rlen,1,1]];
(*Find out if the point is on the grid; *)
(*radgrid=grid[[All,1,1]];*)
(*for tests*)intercomp[rt_,θt_,ϕt_]:=(If[(rt<rmin)||rt>rmax,Print["Outside of radial extent"]];
ra=1;rb=rlen;
While[rb-ra>1,rx=Quotient[rb+ra,2];If[rt>(rfin=grid[[rx,1,1]]),ra=rx,rb=rx]];
{rlo=grid[[ra,1,1]],rhi=grid[[rb,1,1]]};
rman=((rt-rlo)/(rhi-rlo));{thmin=rman grid[[rb,1,2]]+(1-rman)grid[[ra,1,2]],thmax=π-thmin};
If[(θt<thmin)||θt>thmax,Print["Outside of polar extent"]];
tha=1;thb=thlen;
While[thb-tha>1,thx=Quotient[thb+tha,2];If[θt>(thfin=rman grid[[rb,thx,2]]+(1-rman)grid[[ra,thx,2]]),tha=thx,thb=thx]];
{thb,tha};
{thlo=rman grid[[rb,tha,2]]+(1-rman)grid[[ra,tha,2]],thhi=rman grid[[rb,thb,2]]+(1-rman)grid[[ra,thb,2]]};
thman=((θt-thlo)/(thhi-thlo));phx=((ϕt)/(2π))thlen;pha=Floor[phx];phb=pha+1;phman=phx-pha;
ρ=(1-thman)(rman res2m[[rb,tha,5,1]]+(1-rman)res2m[[ra,tha,5,1]])+thman(rman res2m[[rb,thb,5,1]]+(1-rman)res2m[[ra,thb,5,1]]);(*not ϕ-dependent*)
ρ=(1-phman)((1-thman)(rman res2m[[rb,tha,pha,1]]+(1-rman)res2m[[ra,tha,pha,1]])+thman(rman res2m[[rb,thb,pha,1]]+(1-rman)res2m[[ra,thb,pha,1]]))+phman
((1-thman)(rman res2m[[rb,tha,phb,1]]+(1-rman)res2m[[ra,tha,phb,1]])+thman(rman res2m[[rb,thb,phb,1]]+(1-rman)res2m[[ra,thb,phb,1]]));
ρx=((3-rman-thman-phman )res2m[[ra,tha,pha,1]]+rman res2m[[rb,tha,pha,1]]+thman res2m[[ra,thb,pha,1]]+phman res2m[[ra,tha,phb,1]])/3(*ϕ-dependent*))
(*ρ=(1-thman)(rman res2m[[rb,tha,5,1]]+(1-rman)res2m[[ra,tha,5,1]])+thman(rman res2m[[rb,thb,5,1]]+(1-rman)res2m[[ra,thb,5,1]]);(*not ϕ-dependent*)*)
tab1x=Table[intercomp[rt,θt,ϕt];{ϕt,ρx},{ϕt,0.9,0.9+2π(4/phlen),(0.1/phlen)}];
tab2x=Table[intercomp[rt,θt,ϕt];{θt,ρx},{θt,0.9,0.9+2π(4/thlen),(0.1/thlen)}];
tab3x=Table[intercomp[rt,θt,ϕt];{rt,ρx},{rt,1.9,2.9,0.02}];
tab1=Table[intercomp[rt,θt,ϕt];{ϕt,ρ},{ϕt,0.9,0.9+2π(4/phlen),(0.1/phlen)}];
tab2=Table[intercomp[rt,θt,ϕt];{θt,ρ},{θt,0.9,0.9+2π(4/thlen),(0.1/thlen)}];
tab3=Table[intercomp[rt,θt,ϕt];{rt,ρ},{rt,1.9,2.9,0.02}];
ListPlot[{tab1,tab1x},Joined->True,PlotRange->All,PlotStyle->{Blue,Red}]
ListPlot[{tab2,tab2x},Joined->True,PlotRange->All,PlotStyle->{Blue,Red}]
ListPlot[{tab3,tab3x},Joined->True,PlotRange->All,PlotStyle->{Blue,Red}]

(*computation of all matrices*)Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);
(gmn=-{{((Σsq)/(ρsq Δ)),0,0,((2a r)/(ρsq Δ))},{0,-(Δ/(ρsq)),0,0},{0,0,-((ρsq)^(-1)),0},{((2a r)/(ρsq Δ)),0,0,-(((Δ-(a^2)((Sin[θ])^2)))/(ρsq Δ ((Sin[θ])^2)))}})//MatrixForm;
(invgmn=FullSimplify[Inverse[gmn],{r>0}])//MatrixForm;
x[1]=t;x[2]=r;x[3]=θ;x[4]=ϕ;
(*G[1,1,1]=Sum[gmn[[1,m]]*(D[invgmn[[m,1]],x[1]]+D[invgmn[[m,1]],x[1]]-D[invgmn[[1,1]],x[m]]),{m,1,4}];
Table[G[i,k,l]=Simplify[(1/2)Sum[gmn[[i,m]]*(D[invgmn[[m,k]],x[l]]+D[invgmn[[m,l]],x[k]]-D[invgmn[[k,l]],x[m]]),{m,1,4}]],{i,1,4},{k,1,4},{l,1,4}];*)
gmn=Simplify[gmn];
(invgmn=Simplify[invgmn/.Cos[2θ]->2((Cos[θ])^2)-1])//MatrixForm;
(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;
(Anumu={{1,0,0,0},{((2r)/Δ),1,0,(a/Δ)},{0,0,1,0},{0,0,0,1}});
(invAnumu=Transpose[Inverse[Anumu]])//MatrixForm;

(*smart inverse lookup*)sinter[rt_,θt_,ϕt_]:=(ra=1;rb=rlen;While[rb-ra>1,rx=Quotient[rb+ra,2];If[rt>(rfin=grid[[rx,1,1]]),ra=rx,rb=rx]];
{rlo=grid[[ra,1,1]],rhi=grid[[rb,1,1]]};
rman=((rt-rlo)/(rhi-rlo));{thmin=rman grid[[rb,1,2]]+(1-rman)grid[[ra,1,2]],thmax=π-thmin};
tha=1;thb=thlen;While[thb-tha>1,thx=Quotient[thb+tha,2];If[θt>(thfin=rman grid[[rb,thx,2]]+(1-rman)grid[[ra,thx,2]]),tha=thx,thb=thx]];
{thb,tha};{thlo=rman grid[[rb,tha,2]]+(1-rman)grid[[ra,tha,2]],thhi=rman grid[[rb,thb,2]]+(1-rman)grid[[ra,thb,2]]};
thman=((θt-thlo)/(thhi-thlo));phx=((ϕt)/(2π))thlen;pha=Floor[phx];phb=pha+1;phman=phx-pha;)

rt=2.3;θt=2.0;ϕt=1.2;sinter[rt,θt,ϕt];
avggd=(1-thman)(rman resgd[[rb,tha,All]]+(1-rman)resgd[[ra,tha,All]])+thman(rman resgd[[rb,thb,All]]+(1-rman)resgd[[ra,thb,All]])

rtp=avggd[[8;;10]];
xxx=avggd[[5;;7]];
Print["Coord (x1,x2): " ,xxx[[1;;2]],"; (r,θ): ", rtp[[1;;2]]];
sub={r->rtp[[1]],x1->xxx[[1]],θ->rtp[[2]],x2->xxx[[2]],a->0.0};

gcon=Partition[avggd[[11;;26]],4];gcov=Partition[avggd[[27;;42]],4];
gdet=avggd[[43]];dxdxp=Partition[avggd[[48;;63]],4];
fie=((3-rman-thman-phman )res2m[[ra,tha,pha]]+rman res2m[[rb,tha,pha]]+thman res2m[[ra,thb,pha]]+phman res2m[[ra,tha,phb]])/3;
(*imprecise \[Implies] 1*(10^(-2)) for rt=2.6;θt=2.1;ϕt=1.2; when matrices except transformation are pre-computed*)
fie=(1-phman)((1-thman)(rman res2m[[rb,tha,pha]]+(1-rman)res2m[[ra,tha,pha]])+thman(rman res2m[[rb,thb,pha]]+(1-rman)res2m[[ra,thb,pha]]))+phman((1-thman)(rman res2m[[rb,tha,phb]]+(1-rman)res2m[[ra,tha,phb]])+thman(rman res2m[[rb,thb,phb]]+(1-rman)res2m[[ra,thb,phb]]));(*precise multilinear \[Implies] 3*(10^(-3)) for rt=2.6;θt=2.1;ϕt=1.2;  when matrices except transformation are pre-computed*)

Print["ρ=", ρ=fie[[1]],"; U=",U=fie[[2]]];
Print["(u^μ)=",uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]];ulo=gcov.uup;
Print["MKS (B^i)=",bi=fie[[9;;11]]];
Print["(b^μ)=",bup=Flatten[{bi.ulo[[2;;4]],((bi+bi.ulo[[2;;4]]uup[[2;;4]])/(uup[[1]]))}],"; (Subscript[(b^μ), KS])=",bupKS=dxdxp.bup];
Print["In KS coordinates: (Subscript[(u^μ), KS])=", uupKS=dxdxp.uup,";  (Subscript[u, μKS])=",uloKS=invgmnKS.uupKS/.sub];
Print["In BL coordinates: (Subscript[(u^μ), BL])=", uupBL=invAnumu.uupKS/.sub,";  (Subscript[u, μBL])=",uloBL=invgmn.uupBL/.sub];
Print["Back to KS coordinates: (Subscript[(u^μ), KS])=", Transpose[Anumu].uupBL/.sub];

Print["Test 1: (g^(μν))(Subscript[g, νσ])==(Subscript[(δ^μ), σ]); should be identity matrix:   ",Chop[gcov.gcon//MatrixForm]]
Print["Test 2: gdet==(Sqrt[-Det[gcov]]);  should be unity:      ",((Sqrt[-Det[gcov]])/gdet)]
Print["Test 3: normalization of velocity in MKS: (u^μ)(Subscript[u, μ])=-1:  ",ulo.uup]
Print["Test 4: normalization of velocity in KS: (Subscript[(u^μ), KS])(Subscript[u, (Subscript[μ, KS])])=-1:  ",uupKS.uloKS]
Print["Test 5: normalization of velocity in BL: (Subscript[(u^μ), BL])(Subscript[u, (Subscript[μ, BL])])=-1:  ",uupBL.uloBL]

Thread[{u[0],u[1],u[2],u[3]}=uupKS];Thread[{uco[0],uco[1],uco[2],uco[3]}=uloKS];
xx[0]={u[0],u[1],u[2],u[3]};xxlo[0]=-invgmnKS.xx[0]/.sub;
xx[2]={uco[2]u[0],uco[2]u[1],uco[2]u[2]+1,uco[2]u[3]};xxlo[2]=invgmnKS.xx[2]/.sub;
xx[3]={-((uco[3])/(uco[0])),0,0,1};xxlo[3]=invgmnKS.xx[3]/.sub;
sc[0]=(Sqrt[xx[0].xxlo[0]])/.sub;sc[3]=(Sqrt[xx[3].xxlo[3]])/.sub;sc[2]=(Sqrt[xx[2].xxlo[2]])/.sub;
xx[0]/=sc[0];xxlo[0]/=sc[0]/.sub;xx[2]/=sc[2];xxlo[2]/=sc[2]/.sub;xx[3]/=sc[3];xxlo[3]/=sc[3]/.sub;

xpr[1]={u[0]uco[1] ,-(u[0]uco[0]+u[3]uco[3]),0,u[3]uco[1]}/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]=xpr[1]-xlopr[1].xx[3]*xx[3];
xxlo[1]=invgmnKS.xx[1]/.sub;sc[1]=(Sqrt[xx[1].xxlo[1]])/.sub;
xlopr[1]=invgmnKS.xpr[1]/.sub;xx[1]/=sc[1];xxlo[1]/=sc[1]/.sub;
trxx={xxlo[0],xxlo[1],xxlo[2],xxlo[3]};trcoxx=Transpose[{xx[0],xx[1],xx[2],xx[3]}];
(*tests*)
Print["Test 6: (Subscript[u, i])=",trxx.{u[0],u[1],u[2],u[3]},";  should be {1,0,0,0}"];
Print["Test 7: co-moving (Subscript[B, i])=",coB=trxx.bupKS/.sub,";  0-th component should be 0"];
Print["Back to KS frame: (Subscript[(b^μ), KS])=",trcoxx.coB];

dxdxp
Inverse[dxdxp]

Transpose[Anumu].{kBL[0],kBL[1],kBL[2],kBL[3]}
Δ

3. (10^10)/(4 (10^12))

new simulations till 10^5M radial grid

xxx[[1;;2]]

dir = NotebookDirectory[];
read11[fnum_] := (namefile ="E:\\GRMHD\\new\\a0\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*256*128*64*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];
   Close[namefile];);
read11[9900]

stt=StringToStream[strrr];
Do[x1min=Read[stt,Number],{i,5}]
Do[dx1=Read[stt,Number],{i,3}]
Do[off=Read[stt,Number],{i,6}]
Close[stt];
rtab=off+Exp[x1min+dx1 Range[0,255]]

strrr

{x1min,dx1,off}

Limiting the polar region

dir = NotebookDirectory[];
ndd=650;ncut=120;NDIM=4;
dirin="E:\\GRMHD\\Jon\\";anum=0.9375;
readgd[spstr_] := (namefile =dirin<>atab<>"\\gdump.bin";
   siz = FileByteCount[namefile]; xs = 272*128*256*8*126;
   fi11 = OpenRead[namefile, BinaryFormat -> True];strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
SetStreamPosition[fi11,siz-xs];
(*resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],(NDIM^2)],272],{2,1,3}];*)
resgd = Transpose[Partition[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],126],272],128],{3,2,1,4}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[fi11];);
readgd2d[spstr_] := (namefile =dirin<>atab<>"\\gdump.bin";
   siz = FileByteCount[namefile]; phlen=256;xs = 272*128*phlen*8*126;
   fi11 = OpenRead[namefile, BinaryFormat -> True];strrr = StringJoin[BinaryReadList[fi11, "Character8", 156]];
SetStreamPosition[fi11,siz-xs];
(*resgd = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8],(NDIM^2)],272],{2,1,3}];*)
resgd2d = Transpose[Partition[Partition[BinaryReadList[fi11, "Real64",xs/8/phlen],126],272],{2,1,3}];
   (*resgd = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*64*128*256], 11], 256], 128], {3, 2, 1, 4}];*)
   Close[fi11];);

read11[spstr_,fnum_] := (namefile =dirin<>spstr<>"\\fieldlines\\fieldline" <> StringJoin[ToString[#] & /@ PadLeft[IntegerDigits[fnum], 4]] <> ".bin";
   siz = FileByteCount[namefile]; xs = 11*272 128 256*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
str=StringToStream[strrr];param=ReadList[str,Number];Close[str];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", xs/4], 11], 272], 128], {3, 2, 1, 4}];
   Close[namefile];);
atab="thickdisk7";
(*readgd[atab]*)
readgd2d[atab]
read11[atab,7034]
Dimensions[res2m]
Dimensions[resgd]
Dimensions[resgd2d]
(*extension to large radius: 
quantities extended as power-laws; 
elements of dxdxp matrices are extended as power-laws; 
grids are extended as power-laws; theta grid is fixed at 40M boundary;*)

rtab=resgd2d[[All,1,7]];
thtab=resgd2d[[1,All,8]];

nm=10;nx=130;
rtab[[1;;nx]]
thtab
Dimensions[ytab=Outer[Times,rtab[[nm;;nx]],Cos[thtab]]]
Dimensions[xtab=Outer[Times,rtab[[nm;;nx]],Sin[thtab]]]
ListPlot[Transpose[{Flatten[xtab],Flatten[ytab]}],PlotRange->All,AspectRatio->2,ImageSize->300]

ListDensityPlot[Transpose[{Flatten[xtab],Flatten[ytab],Flatten[res2m[[nm;;nx,All,1,1]]]}],PlotRange->All,AspectRatio->2,ImageSize->300]

(180/π)thtab[[10]]

nm=1;nx=130;
ListLogLogPlot[{Transpose[{rtab[[nm;;nx]],((res2m[[nm;;nx,10,1,2]])/(res2m[[nm;;nx,10,1,1]]))}],Transpose[{rtab[[nm;;nx]],Log[1+((res2m[[nm;;nx,10,1,2]])/(res2m[[nm;;nx,10,1,1]]))]}]},ImageSize->450,Joined->True,PlotStyle->{Blue,Red}]

ListContourPlot[Transpose[{Flatten[xtab],Flatten[ytab],Log[Flatten[res2m[[nm;;nx,All,1,1]]]]}],PlotRange->All,AspectRatio->2,ImageSize->300]

ListPlot3D[Transpose[{Flatten[xtab],Flatten[ytab],Log[Flatten[res2m[[nm;;nx,All,1,2]]]]}],PlotRange->All,AspectRatio->2,ImageSize->300]

imp=Import["E:\\GRMHD\\Jon\\thickdisk7\\Tsmapthickdisk7.dat"][[All,1;;2]];
rtabx=imp[[All,1]];
Uρtab=((imp[[All,2]])/(mp*cc*cc/3/kb/.ruru));

Max[res2m[[1,All,1,2]]]
Position[res2m[[1,All,1,2]],Max[res2m[[1,All,1,2]]]][[1]]
ListLogPlot[{Transpose[{thtab,res2m[[1,All,1,2]]}],Transpose[{thtab,res2m[[7,All,1,2]]}],Transpose[{thtab,res2m[[12,All,1,2]]}]},Joined->True];
ListLogPlot[{Transpose[{thtab,((res2m[[1,All,1,2]])/(res2m[[1,All,1,1]]))}],Transpose[{thtab,((res2m[[7,All,1,2]])/(res2m[[7,All,1,1]]))}],Transpose[{thtab,((res2m[[12,All,1,2]])/(res2m[[12,All,1,1]]))}],Transpose[{thtab,((res2m[[20,All,1,2]])/(res2m[[20,All,1,1]]))}]},Joined->True,PlotStyle->{Black,Blue,Green,Red}]
Print[" ",Uρtab[[1]]," ",Uρtab[[7]]," ",Uρtab[[12]]," ",Uρtab[[20]]]

rtab[[1]]
rtab[[7]]
rtab[[12]]
rtab[[20]]

{cpl1=ListContourPlot[Transpose[{Flatten[xtab],Flatten[ytab],Log[Flatten[((res2m[[nm;;nx,All,1,2]])/(res2m[[nm;;nx,All,1,1]]))]]}],PlotRange->All,AspectRatio->2,ImageSize->300,Contours->30],
cpl2=ListContourPlot[Transpose[{Flatten[xtab],Flatten[ytab],Log[Log[1+Flatten[((res2m[[nm;;nx,All,1,2]])/(res2m[[nm;;nx,All,1,1]]))]]]}],PlotRange->All,AspectRatio->2,ImageSize->300,Contours->30]}

dcos1=0.02;Print["Angle = ",(180/π)ArcCos[1-dcos1],"°"]
dcos2=0.05;Print["Angle = ",(180/π)ArcCos[1-dcos2],"°"]
dcos3=0.1;Print["Angle = ",(180/π)ArcCos[1-dcos3],"°"]
Print["Jet boundary from numerical simulations at 50M radius: ",(180/π)ArcTan[(16./50)],"°"]
rx=rtab[[nx]];
pl1=Plot[{Cot[ArcCos[1-dcos1]]t,-Cot[ArcCos[1-dcos1]]t},{t,0,rx},PlotRange->{{0,rx},{-rx,rx}},AspectRatio->2,PlotStyle->{{Blue,AbsoluteThickness[3]},{Blue,AbsoluteThickness[3]}}];
pl2=Plot[{Cot[ArcCos[1-dcos2]]t,-Cot[ArcCos[1-dcos2]]t},{t,0,rx},PlotRange->{{0,rx},{-rx,rx}},AspectRatio->2,PlotStyle->{{Red,AbsoluteThickness[3]},{Red,AbsoluteThickness[3]}}];
pl3=Plot[{Cot[ArcCos[1-dcos3]]t,-Cot[ArcCos[1-dcos3]]t},{t,0,rx},PlotRange->{{0,rx},{-rx,rx}},AspectRatio->2,PlotStyle->{{Green,AbsoluteThickness[3]},{Green,AbsoluteThickness[3]}}];
Show[cpl1,pl1,pl2,pl3]

Dimensions[res2m]

(180/π)ArcCos[1-0.02]
(180/π)ArcCos[1-0.05]
(180/π)ArcCos[1-0.1]

(0.249/0.194)
(1.484/1.101)

((1000^2)/(4π (10^6.5)mp (cc^2)))/.ruru

m_space functions

Minimization of χ^2

(*matrix inversion*)ord=2;ma=Table[matr[i][j],{i,0,ord-1},{j,0,ord-1}];
bta=Table[bb[i],{i,0,ord-1}];
-Inverse[ma].bta//MatrixForm

ord=3;ma=Table[matr[i][j],{i,0,ord-1},{j,0,ord-1}];
bta=Table[bb[i],{i,0,ord-1}];
den=Denominator[Inverse[ma][[1,1]]];
fuss=FullSimplify[-det Inverse[ma].bta/.{1/den->1/det}]
Print[InputForm[den]]
Do[Print[InputForm[fuss[[i]]]],{i,ord}]

m_conv functions

test steps for minimization of square sums

averaged over 3

(*dheat step*)
stheat={{0.192,6.67601},{0.096,2.67034},{0.048,0.69097},{0.024,0.17334},{0.012,0.04627},{0.006,0.01321},{0.003,0.00402},{0.0015,0.00114}};
afin=a/.FindFit[stheat,a (x^2),{a},x];
lpl=ListLogLogPlot[stheat,Joined->True,PlotStyle->Red];
pl=LogLogPlot[1.6afin (x^2),{x,stheat[[-1,1]],stheat[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs dheat"]

(*drho step*)
strho={{0.16,1.49231},{0.08,0.56627},{0.04,0.17515},{0.02,0.04414},{0.01,0.01266},{0.005,0.00302},{0.0025,0.00154}};
afin=a/.FindFit[strho,a (x^2),{a},x];
lpl=ListLogLogPlot[strho,Joined->True,PlotStyle->Red];
pl=LogLogPlot[1.85afin (x^2),{x,strho[[-1,1]],strho[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs drho"]

(*dtheta step*)
sttheta={{0.08,2.42186},{0.04,0.53909},{0.02,0.53535},{0.01,0.18599},{0.005,0.07712},{0.0025,0.20468},{0.00125,0.09617},{0.000625,0.01878}};
afin=a/.FindFit[sttheta,a (x^2),{a},x];
lpl=ListLogLogPlot[sttheta,Joined->True,PlotStyle->Red];
pl=LogLogPlot[4afin (x^2),{x,sttheta[[-1,1]],sttheta[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs dtheta"]

averaged over 11

(*dheat step*)
stheat={{0.192,7.48061},{0.096,1.80863},{0.048,0.42979},{0.024,0.10569},{0.012,0.02651},{0.006,0.00697},{0.003,0.00186},{0.0015,0.00044}};
afin=a/.FindFit[stheat,a (x^2),{a},x];
lpl=ListLogLogPlot[stheat,Joined->True,PlotStyle->Red];
pl=LogLogPlot[0.95afin (x^2),{x,stheat[[-1,1]],stheat[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs dheat"]

(*drho step*)
strho={{0.32,6.57651},{0.16,1.48534},{0.08,0.33530},{0.04,0.08119},{0.02,0.01985},{0.01,0.00540},{0.005,0.00151},{0.0025,0.00041}};
afin=a/.FindFit[strho,a (x^2),{a},x];
lpl=ListLogLogPlot[strho,Joined->True,PlotStyle->Red];
pl=LogLogPlot[0.8afin (x^2),{x,strho[[-1,1]],strho[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs drho"]

(*dtheta step*)
sttheta={{0.08,1.06372},{0.04,0.52826},{0.02,0.13402},{0.01,0.06745},{0.005,0.03317},{0.0025,0.03555},{0.00125,0.01417},{0.000625,0.00258}};
afin=a/.FindFit[sttheta,a (x^2),{a},x];
lpl=ListLogLogPlot[sttheta,Joined->True,PlotStyle->Red];
pl=LogLogPlot[4afin (x^2),{x,sttheta[[-1,1]],sttheta[[1,1]]},PlotStyle->Blue];
Show[pl,lpl,Frame->True,ImageSize->500,PlotLabel->"(χ^2) vs dtheta"]

Evolving the simulation along the geodesics Feb2012

6850 - 6950=6900+/-50  - simulation numbers, whereas times are 13700-13900=13800+/-100 M with time step 2M; can also try time step 4M.
Geodesics start at t=0 at certain location and are near the BH some time around t=1;
timescales should be multiplied by 20,000 to get time in M, and multiplied by 10,000 to get time in ”fnum”.
Thus, need to consider times from tmin=0.995 till tmax=1.005
fdiff=50;//subject to change from 0 to 50; at fdiff=0 - dealing with a single snapshot.
fnmin=fnum-fdiff;// corresponds to tmax; for t>tmax - use fnmin snapshot
fnmax=fnum+fdiff;//corresponds to tmin; for t<tmin - use fnmax snapshot.
Multilinear in time - find nearest 2 neighbours and interpolate.

1+((3 Sinh[-9.498])/20000)
1+((3 Sinh[9.5])/20000)

81 4 32 64 256 11

2147483647

Print[N[Integrate[(1/(Sqrt[1-(rg/r)]))-1,{r,rg,20000rg},GenerateConditions->False]]," - time difference for ray propagation with/without the BH"]

(92/1.5)

Throughput

((5(*spins*)7(*models*)21(*periods*)120(*snapshots*)4(*computations of spectrum*)23(*MB - size of snapshot*))/(3 60 60.))(*MB/sec - actual*)
5 7((120 23)/60.)(*MB/sec - needed*)
Print[CPUhrs=(1./3600)4(*cores*)41(*s - per single computation on 4 cores; *)21(*evaluations*)40(*values of temperature*)(7+6)(*values of density + convergence*)50(*values of inclination angle*)(5+6)(*values of spin + convergence*)," CPU-hours"]
Print["Can run ",jobs=20(*4-core jobs*)," 4-core jobs in parallel"]
Print["Then computation is done in ",(CPUhrs/(4. jobs 24(*hrs a day*)))," days"]
Print["One thread needs ",(1./(24 3600))41(*s - per single computation on 4 cores; *)21(*evaluations*) 40(*values of temperature*)(7+6)(*values of density + convergence*)," days"]
(*3hr10min for 13x21 evaluations \[Implies] 41s per evaluation*)

60((3 60+10)/(13 21.))

Characterizing best fits

if(sp==0){rhonor=1030713.25;heat=0.42107;th=0.7328;fdiff=60;}//4.2802
if(sp==1){rhonor=988988.82;heat=0.37012;th=1.8407;fdiff=60;}//2.9363
if(sp==2){rhonor=817113.52;heat=0.37239;th=2.0159;fdiff=60;}//3.0837
if(sp==3){rhonor=656228.94;heat=0.39849;th=2.2074;fdiff=60;}//4.1912
if(sp==4){rhonor=397281.19;heat=0.41343;th=2.1439;fdiff=60;}//3.0059

(180/π)(π-{π-0.7328,1.8407,2.0159,2.2074,2.1439})

de

π/2.*(2*cas-1)/30/.cas->21

dir="G:\\mathematica\\3DMHD\\a";atab={"0","05","07","09","098"};
read11[n_]:=(namefile=dir<>atab[[n]]<>"\\fieldline00-1.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)
res2m=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*32*64*256],11],256],64],{3,2,1,4}];
Close[namefile];
namefile =dir<>atab[[n]]<>"\\fieldline0000.bin";
   siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   (*res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];*)
   Close[namefile];);
read11[3]

Total[res2m[[All,All,All,All]],3]
Union[Flatten[res2m]][[-100;;-1]]

Check spin 0.5 simulation

dir="G:\\GRMHD\\a05\\fieldlines\\1\\fieldline";
dir="G:\\GRMHD\\a05\\fieldlines\\new2\\fieldline";
dir="G:\\GRMHD\\a05\\fieldlines\\new\\fieldline";
dir="G:\\GRMHD\\a05\\fieldlines\\new3\\fieldline";
read11[n_]:=(namefile=dir<>ToString[n]<>".bin";
fi11=OpenRead[namefile,BinaryFormat->True];siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];
   Close[namefile];);

Dimensions[res2m]

ListDensityPlot[Log[Transpose[res2m[[1;;160,32,All,1]]]],ColorFunction->"Rainbow"]

Dynamic[i]
Table[read11[i];ListDensityPlot[((Transpose[res2m[[1;;160,All,1,1]]])/20),ColorFunctionScaling->False],{i,4950,9950,100}]

ListDensityPlot[((Transpose[res2m[[1;;160,All,1,1]]])/20),Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]}]

ListAnimate

Dynamic[i](*ρ and U near the jump*)
Do[read11[i];pl=ListDensityPlot[0.5+0.2Log[Transpose[res2m[[1;;160,All,1,1]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];
plen=ListDensityPlot[1+0.2Log[Transpose[res2m[[1;;160,All,1,2]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];Export["G:\\GRMHD\\movies\\rate05\\img"<>ToString[i]<>".jpg",pl];Export["G:\\GRMHD\\movies\\rate05en\\img"<>ToString[i]<>".jpg",plen],{i,6550(*4950*),6950(*9950*),1}]

Dynamic[i](*ρ and U over the whole simulation*)
Do[read11[i];pl=ListDensityPlot[0.5+0.2Log[Transpose[res2m[[1;;160,All,1,1]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];
plen=ListDensityPlot[1+0.2Log[Transpose[res2m[[1;;160,All,1,2]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];Export["G:\\GRMHD\\movies\\rate05\\img"<>ToString[i]<>".jpg",pl];Export["G:\\GRMHD\\movies\\rate05en\\img"<>ToString[i]<>".jpg",plen],{i,4950,9950,10}]

(*gamma factor near the jump and in the whole simulation*)Do[read11[i];plgam=ListDensityPlot[Log[Transpose[res2m[[1;;160,All,1,5]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];
Export["G:\\GRMHD\\movies\\rate05gam\\img"<>ToString[i]<>".jpg",plgam],{i,6550(*4950*),6950(*9950*),1}]
Do[read11[i];plgam=ListDensityPlot[Log[Transpose[res2m[[1;;160,All,1,5]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False];
Export["G:\\GRMHD\\movies\\rate05gam\\img"<>ToString[i]<>".jpg",plgam],{i,4950,9950,10}]

ListDensityPlot[0.5+0.2Log[Transpose[res2m[[1;;160,All,1,1]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False]
ListDensityPlot[1+0.2Log[Transpose[res2m[[1;;160,All,1,2]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False]

Dynamic[i]
Do[read11[i];pl=ListDensityPlot[((Transpose[res2m[[1;;160,32,All,1]]])/20),Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]}];
Export["G:\\GRMHD\\movies\\rate05phi\\img"<>ToString[i]<>".jpg",pl],{i,4950,9950,10}]

read11[550];plgam=ListDensityPlot[Log[Transpose[res2m[[1;;160,All,1,5]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False]

read11[6650];ListDensityPlot[((Transpose[res2m[[1;;160,All,1,8]]])/3),Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow"]
read11[6651];ListDensityPlot[((Transpose[res2m[[1;;160,All,1,8]]])/3),Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow"]

np=32;θxx=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.h->0.65/.x2->(x2tab=N[Range[(1/(4np)),1-(1/(4np)),(1/(2np))]]);
theta=Cos[θxx];rate=0.;nx=10;diff=Differences[theta];

accr=Table[read11[i];rate=res2m[[nx,1;;-2,All,1]]*res2m[[nx,1;;-2,All,5]]*res2m[[nx,1;;-2,All,6]];
rate=Total[diff.rate];{i,rate},{i,8020(*4950*),8138(*9950*),2}]

ListLogPlot[accr,Joined->True]

5

4

3

2

1

ListLogPlot[accr,Joined->True]

ListLogPlot[accr,Joined->True]

Table[read11[i];plgam=ListDensityPlot[Log[Transpose[res2m[[1;;160,All,1,5]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False],{i,6547(*4950*),6563(*9950*),2}]

Table[read11[i];plgam=ListDensityPlot[Log[Transpose[res2m[[1;;160,All,1,5]]]],Epilog->{Inset[Style[i,White,Bold,16],ImageScaled[{0.5,0.9}]]},ColorFunction->"Rainbow",ColorFunctionScaling->False],{i,7347,7353,2}]

lis=%

(*anim=ListAnimate[lis];*)
Export["G:\\GRMHD\\movie05.avi",lis]

read11[6597];

Dimensions[res2m]

nm=256;
read11[6597];res6597=res2m[[30;;nm,32,1,5]];
read11[6599];res6599=res2m[[30;;nm,32,1,5]];
read11[6601];res6601=res2m[[30;;nm,32,1,5]];

{ListLogPlot[{res6597,res6599,res6601}-1,Joined->True],
ListLogPlot[1/((res6601-1)/(res6599-1)),Joined->True]}

rho,u,-u^t,-T^r_t/(rho u^r),u^t,v^r,v^theta,v^phi,B^r,B^theta,B^phi

ruru={rg->((6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2)),arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.67 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12))};
torb=40((2π)/(Sqrt[rg (cc^2)]))(((20rg))^(3/2))/.ruru
tunit=(rg/cc)/.ruru
torb/tunit

Dynamical model for Lei

Te=8*10^10/(Sqrt[rr]);
θ=2.146755 - inclination angle

rho->n
ut -> u[0]
ur -> u[1]
utheta -> u[2]
uphi -> u[3]
Bi1-> Bi[1]
Bi2-> Bi[2]
Bi3-> Bi[3]
doub blo=(ulo[1]*Bi[1]+ulo[2]*Bi[2]+ulo[3]*Bi[3]);
Bup[0]=blo;for(m=1;m<4;m++)Bup[m]=(Bi[m]+blo*u[m])/u[0];

dir="G:\\mathematica\\3DMHD\\a";atab={"0","05","07","09","098"};
read11[n_]:=(namefile=dir<>atab[[n]]<>"\\fieldline00-1.bin";
fi11=OpenRead[namefile,BinaryFormat->True];
(*Print[StringJoin[BinaryReadList[fi11,"Character8",374]]];*)
res2m=Transpose[brl=Partition[Partition[Partition[BinaryReadList[fi11,"Real32",11*32*64*256],11],256],64],{3,2,1,4}];
Close[namefile];
namefile =dir<>atab[[n]]<>"\\fieldline0000.bin";
   siz = FileByteCount[namefile]; xs = 11*256*64*32*4;
   fi11 = OpenRead[namefile, BinaryFormat -> True];
   strrr = StringJoin[BinaryReadList[fi11, "Character8", FileByteCount[namefile] - xs]];
   (*res2m = Transpose[llist = Partition[Partition[Partition[BinaryReadList[fi11, "Real32", 11*32*64*256], 11], 256], 64], {3, 2, 1, 4}];*)
   Close[namefile];);
read11[5]
stt=StringToStream[strrr];
Do[x1min=Read[stt,Number],{i,5}]
Do[dx1=Read[stt,Number],{i,3}]
Do[off=Read[stt,Number],{i,6}]
Close[stt];
x1min=-0.238785893239999`;off=0.3;
dx1=0.0216657617755842`;
rtab=off+Exp[x1min+dx1 Range[0,255]]

(*{x1min,x2min,x3min}={BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"]};
{dx1,dx2,dx3}={BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"],BinaryRead[fi11,"Real32"]};
Do[BinaryRead[fi11,"Real32"],{i,8}];BinaryRead[fi11,"Real32"];
ress=Transpose[Partition[Partition[Partition[BinaryReadList[fi11,"Real32",30*32*64*256],30],256],64],{3,2,1,4}];
Close[namefile];*)

Dimensions[avgres=((Total[res2m[[All,All,All,All]],{3}])/32)]
(*uup=Flatten[{fie[[5]],fie[[5]]fie[[6;;8]]}]*)
lpl3=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,6]]}]]];pl3=ListLogLogPlot[Transpose[{rtab,0.7(rtab^(-2.5))}],PlotStyle->Red];
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,8]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,0.3(rtab^(-1.5))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,avgres[[All,34,5]]*avgres[[All,34,7]]}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.03(rtab^(-1.5))}],PlotStyle->Red];
GraphicsRow[{Show[lpl3,pl3,PlotLabel->"(u^r)-upper index"],Show[lpl1,pl1,PlotLabel->"(u^ϕ)-upper index"],Show[lpl2,pl2,PlotLabel->"(u^θ)-upper index"]},ImageSize->1200]
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,9]]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,0.8(rtab^(-2.5))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,10]]]}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.01(rtab^(-1.5))}],PlotStyle->Red];
lpl3=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,11]]]}]]];pl3=ListLogLogPlot[Transpose[{rtab,0.8(rtab^(-2.5))}],PlotStyle->Red];
GraphicsRow[{Show[lpl1,pl1,PlotLabel->"(B^r)-upper index"],Show[lpl2,pl2,PlotLabel->"(B^θ)-upper index"],Show[lpl3,pl3,PlotLabel->"(B^ϕ)-upper index"]},ImageSize->1200]
lpl1=ListLogLogPlot[Abs[Transpose[{rtab,Abs[avgres[[All,34,1]]]}]]];pl1=ListLogLogPlot[Transpose[{rtab,6(rtab^(-0.8))}],PlotStyle->Red];
lpl2=ListLogLogPlot[Abs[Transpose[{rtab,((avgres[[All,34,2]])/(avgres[[All,34,1]]))}]]];pl2=ListLogLogPlot[Transpose[{rtab,0.25(rtab^(-0.8))}],PlotStyle->Red];
GraphicsRow[{Show[lpl1,pl1,PlotLabel->"ρ"],Show[lpl2,pl2,PlotLabel->"U"]},ImageSize->800]

Dimensions[tofile=Transpose[Table[avgres,{i,32}],{1,3,2,4}]]
Dimensions[brl]
namefile=dir<>atab[[5]]<>"\\fieldline0-10.bin";
fi22=OpenWrite[namefile,BinaryFormat->True];
BinaryWrite[fi22,tofile,"Real32"]
Close[namefile];

Dimensions[avgres]
Dimensions[tofile]
Dimensions[brl]

ruru={rg->((2 6.67 (10^(-8))4.5 (10^6)2 (10^33))/(((3 (10^10)))^2)),arcsec->((2 π 8.4 (10^3)3.08 (10^18))/(360 60 60)),pc->3.08 (10^18),year->365.35 86400,cc->3. (10^10),Msun->2. (10^33),mp->1.67 (10^(-24)),me->9.1 (10^(-28)),kb->1.38 (10^(-16)),h->6.626 (10^(-27)),ee->4.8 (10^(-10)),keV->1000 1.602 (10^(-12))};
rhonor=6 (10^5);
Bnor=(Sqrt[4.*π*rhonor*mp])*cc/.ruru
rgrav->rg/.ruru

Export[dir<>atab[[5]]<>"\\rho.dat",rhonor avgres[[All,All,1]]];
Export[dir<>atab[[5]]<>"\\ut.dat",avgres[[All,All,5]]];
Export[dir<>atab[[5]]<>"\\ur.dat",avgres[[All,All,5]]*avgres[[All,All,6]]];
Export[dir<>atab[[5]]<>"\\utheta.dat",avgres[[All,All,5]]*avgres[[All,All,7]]];
Export[dir<>atab[[5]]<>"\\uphi.dat",avgres[[All,All,5]]*avgres[[All,All,8]]];
Export[dir<>atab[[5]]<>"\\Bi1.dat",Bnor avgres[[All,All,9]]];
Export[dir<>atab[[5]]<>"\\Bi2.dat",Bnor avgres[[All,All,10]]];
Export[dir<>atab[[5]]<>"\\Bi3.dat",Bnor avgres[[All,All,11]]];

send

Δ=(r^2)-2r+(a^2);ρsq=(r^2)+(a^2)((Cos[θ])^2);Σsq=((((r^2)+(a^2)))^2)-(a^2)Δ ((Sin[θ])^2);(invgmnKS={{-1+((2r)/(ρsq)),((2r)/(ρsq)),0,-((2a r ((Sin[θ])^2))/(ρsq))},{((2r)/(ρsq)),1+((2r)/(ρsq)),0,-a(1+((2r)/(ρsq)))((Sin[θ])^2)},{0,0,ρsq,0},{-((2a r ((Sin[θ])^2))/(ρsq)),-a(1+((2r)/(ρsq)))((Sin[θ])^2),0,((Sin[θ])^2)(ρsq+(a^2)(1+((2r)/(ρsq)))((Sin[θ])^2))}})//MatrixForm;(*lower-index KS metrics*)
AMKS=DiagonalMatrix[{1,(1/(r-R0)),1/(π(h+7(1-h)(((2x2-1))^6))),1/π}];
(invnew=Inverse[AMKS].invgmnKS.Inverse[AMKS])//MatrixForm(*lower-index MKS metrics - used for the dynamical model*)

velocity  extension: (u^r) ∝ (r^(-5/2)); (u^θ) ∝ (r^(-3/2)); (u^ϕ) ∝ (r^(-3/2)); (u^t)∝1;
(Subscript[v, θ])=(u^θ)(Sqrt[(Subscript[g, θθ])]) ∝ (r^(-1/2)); (Subscript[v, ϕ])=(u^ϕ)(Sqrt[(Subscript[g, ϕϕ])]) ∝ (r^(-1/2));(Subscript[v, r])∼(u^r)(Sqrt[(Subscript[g, rr])]) ∝ (r^(-3/2));     (Subscript[u, r]) ∝(r^(-1/2)); (Subscript[u, θ])∝ (r^(1/2)); (Subscript[u, ϕ])∝ (r^(1/2)); (Subscript[u, t])∝ 1;
magnetic   field   extension: Bi1 ∝ (r^(-5/2)); Bi2 ∝ (r^(-3/2)); Bi3 ∝ (r^(-5/2)); then   blo=Sum[(Subscript[u, i]) Bi[i]]=(Subscript[u, θ])Bi2∝(r^(-1));
Bup[m]=(Bi[m]+blo*u[m])/u[0];(B^r)∝ (r^(-5/2))+(r^(-1-5/2)) ∝ (r^(-5/2));(B^θ)∝ (r^(-3/2))+(r^(-1-3/2)) ∝ (r^(-3/2));(B^ϕ)∝ (r^(-5/2))+(r^(-1-3/2)) ∝ (r^(-5/2));
(Subscript[b, r])∝ (r^(-3/2));(Subscript[b, θ])∝(r^(-1/2));(Subscript[b, ϕ])∝(r^(-3/2)) - physical fields

density  n ∝ (r^(-0.8)); temperature T∝ (r^(-1))
((B^2)/nkT)=(((Subscript[b, θ])^2)/nkT)∝ (r^0.8) - growing with radius!? => change Bi2 ∝ (r^(-5/2)); then  blo∝(r^(-2)); (Subscript[b, r])∝ (r^(-3/2));(Subscript[b, θ])∝(r^(-3/2));(Subscript[b, ϕ])∝(r^(-3/2)) 
((B^2)/nkT)∝(r^(-0.2))

(invnew/.{a->0,R0->0,h->0,x2->0,θ->(π/2),π->1})//MatrixForm

BLtoMKS=(({{1, ((2 r)/Del), 0, 0}, {0, (1/(r-R0)), 0, 0}, {0, 0, (1/x2temp), 0}, {0, (a/(Del π)), 0, (1/π)}}));uupMKS=BLtoMKS.uupBL(*multiplication from the right*);
Del=(a^2)+(-2+r) r;
x2temp=(π (h+7 (1-h) (((-1+2 x2))^6)));
h=0.65;R0=0.3;

coordinates

x1min=-0.238785893239999`;off=0.3;
dx1=0.0216657617755842`;
rtab=off+Exp[x1min+dx1 Range[0,255]]

np=32;θxx=(h (2x2-1)+(1-h)(((2x2-1))^7)+1)(π/2)/.h->0.65/.x2->(x2tab=N[Range[(1/(4np)),1-(1/(4np)),(1/(2np))]]);
Cos[θxx]
x2tab
(*ListPlot[Transpose[{x2tab,1/((π/2)(h-7(1-h)(((x2-1))^6)))/.h->0.65/.x2->x2tab}],PlotRange->All]
ListPlot[Transpose[{x2tab,1/(π+(1-h) π Cos[2 π x2])/.h->0.65/.x2->x2tab}]]*)