Sign conventions

Swanson, 2003 is chosen as the reference book:
 (Subscript[ϵ, j])=(Subscript[q, j])/|q| - sign of charge;  (Subscript[ω, cj])=|(Subscript[q, j])|(Subscript[B, 0])/(Subscript[m, j]) - positive cyclotron frequency; magnetic field (Subscript[B, 0]) is directed "up" along z axis; x & y axes are direction of polarization;  p. 24, 37 define Faraday rotation
 
 
Melrose 1991 book has the opposite sign convention!!! He has e^2=e^1 x k, whereas in all other books/papers e^2= k x e^1; k=e^1 x e^2 (Rybicki, 1979, Shcherbakov, 2008, Huang 2009, Swanson 2003)
Melrose 1971 - correct definition k=e^1 x e^2 => 
Huang 2009 - polarized transfer equation(8) coincides with that in Melrose 1991 (14.45, 14.49, 14.52); angle χ - from (a,b) to (e1,e2) clockwise rotation; radiation treated in (e1,e2)-frame, whereas (a,b) - parallel transported
                     cos((Subscript[θ, B]))=k.B;   jVc ∝ k.B

In my code (9 Apr):
shooting in equatorial plain: 
ϕ=π; costh=0; r=1000 =>   k=(1,0,0); e1=(0,1,0) =>   k is in radial direction, against the propagation direction; e1 is in "down" direction in picture plane
ϕ=0; costh=0; r=1000 =>   k=(-1,0,0); e1=(0,1,0) =>   k is in radial direction, against the propagation direction; e1 is in "down" direction in picture plane (same result for rays in polar region)
ϕ=0; costh=0; r=1000 =>   k=(-1,0,0); e1=(0,1,0) => B=(0,0,-1); kpbar>0 (B anti-aligned with direction); kbperp>0 => jVc<0; jQc>0; rQc<0; rVc>0
jVc ∝ -k.B

changes: k now looks along the ray; e2= e1 x k; jVc ∝ k.B

Under switching of viewing direction U changes sign; Q doesn't change; V switches sign;
Integration in the negative direction in affine parameter: V switches sign

Observations: V=left-right<0 (Bower 1999)
The rotation measure is defined to be positive for a magnetic field pointed towards the observer (http://iopscience.iop.org.ezp-prod1.hul.harvard.edu/0004-637X/539/1/480/50718.text.html + other links)
Negative RM<0 => magnetic field B pointed away from the observer!

Theory: Right - positive; left - negative (Melrose 1991 book; Rybicky 1978 book) - along the direction of propagation
Emissivity is positive along magnetic field B => negative for small θ => negative observed V (confirmed)
Along the ray - clockwise (right-handed) Faraday rotation, if magnetic field B is along the ray
Lower ν => more Faraday rotation
Melrose 1991, eq. 14.58(p.190):
Spin axis is aligned with magnetic field B vector (Bob) (contradiction)
Spin axis is anti-aligned with magnetic field B vector (me) (confirmed)
θ=0 corresponds to counterclockwise rotating BH, if we look at it from Earth (me) (confirmed)
Low θ configurations => positive difference ang345-ang230 - data suggests
But in C++ calculations angles change signs => should be negative difference (confirmed for θ=0.1, a=0.5)

Definitions of Stokes & rotation transformations

Primary

plp=100;myf=18;sc=10;ss=8;bc=0.13;edg=0.03;enn=0.13;
Ex=E1 Exp[I (-k t+δ)];Ey=E2 Exp[-I k t];rul={E1->1,E2->1,δ->(π/2),k->3};
Plot[{Re[Ex]/.rul,Re[Ey]/.rul},{t,0,2π},PlotStyle->{Red,Blue},PlotLabel->"Ex - red, Ey - blue",ImageSize->400];
tplrot=ParametricPlot3D[{t,u Re[Ex]/.δ->(π/2)/.rul,u Re[Ey]/.δ->(π/2)/.rul},{t,0,2π},{u,0.5,1},AxesLabel->{"t",Null,Null},Axes->{True,False,False},BaseStyle->{myf,Bold},LabelStyle->{myf+1,Bold},PlotLabel->"     clockwise rotation in fixed plane: δ>0     ",PlotStyle->Green,Mesh->False,PlotPoints->plp,BoxStyle->{AbsoluteThickness[ss]},ImageSize->400,ViewPoint->{40sc,10sc,10sc},ViewVertical->{0,1,0},Epilog->{AbsoluteThickness[4],Red,Arrowheads[{0,0.05}],Arrow[{ImageScaled[{bc,edg}],ImageScaled[{bc,enn}]}],Arrow[{ImageScaled[{bc+0.005,edg}],ImageScaled[{edg,edg}]}],
Text[Style["(e^1)",myf+4,Bold],ImageScaled[{bc+0.04,enn-0.03}]],Text[Style["(e^2)",myf+4,Bold],ImageScaled[{edg+0.04,edg+0.04}]]}]
Ex=E1 Exp[I (k z+δ)];Ey=E2 Exp[I k z];rul={E1->1,E2->1,δ->(π/2),k->3};bc=0.16;edg=0.06;enn=0.18;
Plot[{Re[Ex]/.rul,Re[Ey]/.rul},{z,0,2π},PlotStyle->{Red,Blue},PlotLabel->"Ex - red, Ey - blue",ImageSize->400];
zplrot=ParametricPlot3D[{z,u Re[Ex]/.δ->(π/2)/.rul,u Re[Ey]/.δ->(π/2)/.rul},{z,0,2π},{u,0.5,1},AxesLabel->{"z",Null,Null},Axes->{True,False,False},PlotLabel->"  right-handed rotation along the ray: δ>0  ",PlotStyle->Green,Mesh->False,PlotPoints->plp,BaseStyle->{myf,Bold},LabelStyle->{myf+1,Bold},BoxStyle->{AbsoluteThickness[ss]},ImageSize->400,ViewPoint->{40sc,10sc,10sc},ViewVertical->{0,1,0},Epilog->{AbsoluteThickness[4],Red,Arrowheads[{0,0.05}],Arrow[{ImageScaled[{bc,edg}],ImageScaled[{bc,enn}]}],Arrow[{ImageScaled[{bc+0.005,edg}],ImageScaled[{0.02,edg}]}],Arrow[{ImageScaled[{bc,edg}],ImageScaled[{1.4bc,0.01}]}],
Text[Style["(e^1)",myf+4,Bold],ImageScaled[{bc+0.04,enn-0.03}]],Text[Style["(e^2)",myf+4,Bold],ImageScaled[{edg+0.02,edg+0.04}]],Text[Style["k",myf+4,Bold],ImageScaled[{1.4bc,0.06}]]}]

Export[NotebookDirectory[]<>"tplrot1.eps",Rasterize[tplrot,RasterSize->1000]]
Export[NotebookDirectory[]<>"zplrot1.eps",Rasterize[zplrot,RasterSize->1000]]
(*Export[NotebookDirectory[]<>"tplrot2.eps",tplrot]
Export[NotebookDirectory[]<>"zplrot2.eps",zplrot]*)

Ij=ComplexExpand[Ex Ex\[Conjugate]+Ey Ey\[Conjugate]](*definitions in Rybicky, Lightman 1979*)
Qj=ComplexExpand[Ex Ex\[Conjugate]-Ey Ey\[Conjugate]]
Uj=ComplexExpand[Ex Ey\[Conjugate]+Ey Ex\[Conjugate]]
Vj=(1/I)ComplexExpand[Ex Ey\[Conjugate]-Ey Ex\[Conjugate]](*If δ>0, then Vj>0, what means that right-handed wave \[Implies] Vj>0*)

EE={Ex,Ey};
(It=Simplify[ComplexExpand[Outer[Times,EE,EE\[Conjugate]]]])//MatrixForm
(Itx=Simplify[(1/2){{Ij+Qj,Uj+I Vj},{Uj-I Vj,Ij-Qj}}])//MatrixForm
{{II+QQ,UU+I VV},{UU-I VV,II-QQ}}//MatrixForm(*standard definition of polarization tensor (w/ factor 1/2)*)

Print["Transformation for vectors has the opposite sign: mm=",(mm={{Cos[χ],-Sin[χ]},{Sin[χ],Cos[χ]}})//MatrixForm(*transformatin of vectors - with the inverse matrix*)]
Print["Electric field in new coords: ",Thread[{Eaa,Ebb}==mm.{Exx,Eyy}]]
Thread[{Ea,Eb}=mm.{Ex,Ey}];
Iaj=FullSimplify[ComplexExpand[Ea Ea\[Conjugate]+Eb Eb\[Conjugate]]](*definitions in Rybicky, Lightman 1979*)
Qaj=FullSimplify[ComplexExpand[Ea Ea\[Conjugate]-Eb Eb\[Conjugate]]]
Uaj=FullSimplify[ComplexExpand[Ea Eb\[Conjugate]+Eb Ea\[Conjugate]]]
Vaj=FullSimplify[(1/I)ComplexExpand[Ea Eb\[Conjugate]-Eb Ea\[Conjugate]]]

(RR={{1,0,0,0},{0,Cos[2χ],-Sin[2χ],0},{0,Sin[2χ],Cos[2χ],0},{0,0,0,1}})//MatrixForm
Simplify[Thread[{Iaj,Qaj,Uaj,Vaj}==RR.{Ij,Qj,Uj,Vj}]](*indeed, 2χ rotation is for (I,Q,U,V)*)

Emissivities/absorptivities etc from the response tensor (Subscript[(α^μ), ν]) etc

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->2.99792 (10^10),Msun->2. (10^33),mp->1.26 (10^(-24)),me->9.109 (10^(-28)),kb->1.38066 (10^(-16)),h->6.626176 (10^(-27)),ee->4.8032 (10^(-10)),keV->1000 1.602 (10^(-12))};γmax=5000;
ν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,γmax},{t,((νΩ)/(γ^2)),+∞}];
IjcSs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,γmax},{t,((νΩ)/(γ^2)),+∞}];
QjcSs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(2/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,γmax}];
VjcSs[θe_?NumberQ,νΩ_?NumberQ]:=(NIntegrate[(1/γ)(*((νΩ)/(γ^3))*)BesselK[(1/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,γmax},{t,((νΩ)/(γ^2)),+∞}]+NIntegrate[((νΩ)/(γ^3))BesselK[(1/3),((νΩ)/(γ^2))]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,γmax}]);
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};
quant={θ->1.,B->20,Te->6.9 (10^9),ν->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

νc=((ee B)/(2π me cc));νs=(2/9)νc ((θe)^2)Sin[θ];X=(ν/(νs));
expr=(ns=1)(((Sqrt[2])π (ee^2)νs)/(3BesselK[2,((θe)^(-1))]cc))((((X^(1/2))+(2^(11/12))(X^(1/6))))^2)Exp[-(X^(1/3))];
IjGp=expr/.θe->((kb Te)/(me (cc^2)))/.ruru/.quant

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

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;

(*FindRoot[(ρ^2)-2I ρ ξ-(ξ^2)((Sin[θ])^2)+2((Sin[θ])^2)(1-Cos[β ξ])/(β^2)==0/.{θ->0.5,β->(10^(-3)),ρ->2.5},{ξ,100-10I}]*)

acc=7;roro={ω->1,q->ee,ne->1,m->1,θ->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))}

1

η=.;(rmurnu=FullSimplify[((((β/ω)))^2)Outer[Times,((ω Sin[θ])/β){η(1-Cos[β ξ]),Cos[θ](Sin[β ξ]-β ξ)},((ω Sin[θ])/β){-η(1-Cos[β ξ]),Cos[θ](Sin[β ξ]-β ξ)}]/.(η^2)->1])//MatrixForm

Off[NIntegrate::"eincr",NIntegrate::"slwcon"];
acc=7;roro={ω->1,q->ee,ne->1,m->1,θ->(π/4),ρ->((me (cc^2))/(kb Te))}/.ruru/.quant;
resx=NIntegrate[(Cos[α]+I Sin[α])magn/.roro/.ruru/.quant/.ξ->(Cos[α]+I Sin[α])x/.α->(π/4),{x,0,+∞},MaxRecursion->40,PrecisionGoal->acc]

η->1

η->-1

prec=50;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))};
Do[quant={θ->0.1,B->20,Te->3.8 (10^9),ν->2. (10^10)(1.6^(-i)),ne->1};AppendTo[quant,β->((ee B)/(me cc 2π ν))/.quant/.ruru];If[(β/.quant)>0.7,Break[]];
acc=6;roro={ω->1,q->ee,ne->1,m->1,θ->(π/4),ρ->((me (cc^2))/(kb Te))}/.ruru/.quant;
resx=NIntegrate[(Cos[α]+I Sin[α])(10^16)magn/.roro/.ruru/.quant/.ξ->(Cos[α]+I Sin[α])x/.α->(π/20.),{x,0,300},MaxRecursion->40,PrecisionGoal->acc];
res=(10^(-16))(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]]];
Print[ν/.quant,"  ", {αIT,αQT,αVT,ρQT,ρVT}];
{((αIS)/(αIT)),((αQS)/(αQT)),((αVS)/(αVT)),((ρQS)/(ρQT)),((ρVS)/(ρVT))},{i,17}]

((ee B)/(2π me cc))/.quant/.ruru

Negative RM < 0 => magnetic field B pointed away from the observer!

Theory : Right - positive; left - negative (Melrose 1991 book; Rybicky 1978 book) - along the direction of propagation
Emissivity is positive along magnetic field B => negative for small θ => negative observed V (confirmed)

ttx=SessionTime[];Clear[inum,jnum,θetab,νΩtab,θe,νΩ,partab];γmax=5000;
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,γmax];
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,γmax},{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,γmax}]],{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,γmax}]+NIntegrate[(1/γ)(*((νΩ)/(γ^3))*)BesselK[(1/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])/.{θe->θetab[[i]],νΩ->νΩtab[[j]]},{γ,1,γmax},{t,((νΩ)/(γ^2))/.{νΩ->νΩtab[[j]]},+∞}]],{i,Length[θetab]},{j,Length[νΩtab]}]];
Print[SessionTime[]-ttx]

dir="G:\\mathematica\\3DMHD\\";
Export[dir<>"lookupjIy.dat",N[lookupjI],"Table"]
Export[dir<>"lookupjQy.dat",lookupjQ,"Table"]
Export[dir<>"lookupjVy.dat",lookupjV,"Table"]
Export[dir<>"lookuprVy.dat",lookuprV,"Table"]
Export[dir<>"lookuprQy.dat",lookuprQ,"Table"]

dir="G:\\mathematica\\3DMHD\\";
Export[dir<>"lookupjIz.dat",N[lookupjI],"Table"]
Export[dir<>"lookupjQz.dat",lookupjQ,"Table"]
Export[dir<>"lookupjVz.dat",lookupjV,"Table"]
Export[dir<>"lookuprVz.dat",lookuprV,"Table"]
Export[dir<>"lookuprQz.dat",lookuprQ,"Table"]

lookupjI

ComplexExpand[((2cc)/ν)Im[Integrate[Simplify[Normal[Series[Normal[Series[((I (ee^2)ne ω (ρ^2))/((cc^2)me BesselK[2,ρ]))(rot11 ((BesselK[2,Rex])/(Rex^2))-((((ω/(Ω0))))^2)rot22 ((BesselK[3,Rex])/(Rex^3))),{Ω0,0,1}]],{ρ,+∞,0}]]],{ξ,0,+∞},GenerateConditions->False][[1,2]]]/.{ω->2π ν,Ω0->((ee B)/(me cc))}](*Faraday rotation coefficient ρV; dψ/ds=(1/2)ρV; ρV has coefficient 1 (Melrose, Huang have coefficients (1/2))*)

Integrate[((n ρ )/(BesselK[2,ρ]))(En+1)En Exp[-ρ (En+1)],{En,0,+∞},GenerateConditions->False]
(*Integrate[((n ρ )/(BesselK[2,ρ]))(En+1)(Sqrt[(En^2)+2En]) Exp[-ρ (En+1)],{En,0,+∞},GenerateConditions->False]*)
Integrate[((n ρ )/(BesselK[2,ρ]))γ (Sqrt[(γ^2)-1]) Exp[-ρ γ],{γ,1,+∞},GenerateConditions->False]

Rotation of response tensor (Subscript[(α^μ), ν])

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

1

T11=(({{Cos[α β], - Sin[α β], 0}, {Sin[α β], Cos[α β], 0}, {0, 0, 1}}));Clear[η];
T22=(({{((Sin[θ])^2)((Sin[α β])^2), - ((Sin[θ])^2)Sin[α β](1-Cos[α β]), Sin[θ]Cos[θ]α β Sin[α β]}, {((Sin[θ])^2)Sin[α β](1-Cos[α β]), -((Sin[θ])^2)(((1-Cos[α β]))^2), Sin[θ]Cos[θ]α β(1-Cos[α β])}, {Sin[θ]Cos[θ]α β Sin[α β], - Sin[θ]Cos[θ]α β(1-Cos[α β]), ((Cos[θ])^2)(α^2)(β^2)}}));
(rot11=FullSimplify[invmat.invm90.T11.m90.mat][[1;;2,1;;2]]/.α->ξ )//MatrixForm
(rot22=FullSimplify[invmat.invm90.T22.m90.mat][[1;;2,1;;2]]/.α->ξ)//MatrixForm
T11=(({{Cos[β ξ], η Sin[β ξ], 0}, {-η Sin[β ξ], Cos[β ξ], 0}, {0, 0, 1}}));
T22=(({{-(((1-Cos[β ξ]))^2)((Sin[θ])^2), η (1-Cos[β ξ])Sin[β ξ]((Sin[θ])^2), η β ξ(1-Cos[β ξ])Sin[θ]Cos[θ]}, {-η (1-Cos[β ξ])Sin[β ξ]((Sin[θ])^2), ((Sin[β ξ])^2)((Sin[θ])^2), β ξ Sin[β ξ]Sin[θ]Cos[θ]}, {-η β ξ(1-Cos[β ξ])Sin[θ]Cos[θ], β ξ Sin[β ξ]Sin[θ]Cos[θ], (((β ξ))^2)((Cos[θ])^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;

Clear[η];rmurnu=FullSimplify[(((((Ω0)/ω)))^2)Outer[Times,((ω Sin[θ])/(Ω0)){η(1-Cos[β ξ]),Cos[θ](Sin[β ξ]-β ξ)},((ω Sin[θ])/(Ω0)){-η(1-Cos[β ξ]),Cos[θ](Sin[β ξ]-β ξ)}]/.(η^2)->1]//MatrixForm

invm90.T11.m90//MatrixForm
invm90.T22.m90//MatrixForm

d{E1,E2}/ds=={{(α^11),(α^12)},{-(α^12),α22}}.{E1,E2}
((ΔK)^(αβ))=((4π c)/(ω^2))(Subscript[(α^μ), ν]) => (dE^μ)/ds=i /ν (Subscript[(α^μ), ν]) (E^ν)

(E^α)(dE^( *β))/ds==-i  c/ν (Subscript[(α^( *β)), ν]) (E^( *ν))(E^α)
(E^( *β))(dE^α)/ds==i  c/ν (Subscript[(α^α), ν]) (E^ν)(E^( *β))
<(E^α)(E^( *β))>/ds==i  c/ν ((Subscript[(α^α), ν]) <(E^ν)(E^( *β))>- (Subscript[(α^( *β)), ν])<(E^( *ν))(E^α)>)
<(E^1)(E^( *1))>/ds==i  c/ν ((Subscript[(α^1), ν]) <(E^ν)(E^( *1))>- (Subscript[(α^( *1)), ν])<(E^( *ν))(E^1)>)
<(E^1)(E^( *1))>/ds==i  c/ν (((Subscript[(α^1), 1])-(Subscript[(α^( *1)), 1])) <(E^1)(E^( *1))>+(Subscript[(α^1), 2])<(E^2)(E^( *1))>-(Subscript[(α^( *1)), 2])<(E^1)(E^( *2))>)

rux={A[n_,m_]->ReA[n,m]+I ImA[n,m]};
α={{A[1,1],A[1,2]},{-A[1,2],A[2,2]}}
α-α\[Conjugate]
ComplexExpand[α[[1,1]]-α[[1,1]]\[Conjugate]/.rux]

(*old definition*)Clear[α];rux={A[n_,m_]->ReA[n,m]+I ImA[n,m]};
(*Table[Sum[I  c/ν ((Subscript[(α^α), μ]) ((E^μ)(E^( *β)))- (Subscript[(α^( *β)), μ])((E^α)(E^( *μ)))),{μ,2}],{α,2},{β,2}]*)
sisi=Simplify[ComplexExpand[Table[Sum[I  c/ν (A[α,μ]EE[μ,β]-A[β,μ]\[Conjugate]EE[α,μ]),{μ,2}],{α,2},{β,2}]/.{A[2,1]->-A[1,2]}/.rux]];
Esol=Solve[{Ii==EE[1,1]+EE[2,2],Qi==EE[1,1]-EE[2,2],Ui==EE[1,2]+EE[2,1],Vi==(1/I)(EE[1,2]-EE[2,1])},Flatten[Table[EE[n,m],{n,2},{m,2}]]][[1]]
norm=(ν/c);mat=FullSimplify[norm sisi/.Esol]
ans=FullSimplify[{(*dIds==*)mat[[1,1]]+mat[[2,2]],(*dQds==*)mat[[1,1]]-mat[[2,2]],(*dUds==*)mat[[1,2]]+mat[[2,1]],(*dVds==*)(1/I)(mat[[1,2]]-mat[[2,1]])}];
Xi={Ii,Qi,Ui,Vi};
-Table[D[ans[[k]],Xi[[l]]],{k,4},{l,4}]//MatrixForm

(*new definition*)Clear[α];rux={A[n_,m_]->ReA[n,m]+I ImA[n,m]};
(*Table[Sum[I  c/ν ((Subscript[(α^α), μ]) ((E^μ)(E^( *β)))- (Subscript[(α^( *β)), μ])((E^α)(E^( *μ)))),{μ,2}],{α,2},{β,2}]*)
sisi=Simplify[ComplexExpand[Table[Sum[I /ν (A[α,μ]EE[μ,β]-A[β,μ]\[Conjugate]EE[α,μ]),{μ,2}],{α,2},{β,2}]/.{A[2,1]->-A[1,2]}/.rux]];
Esol=Solve[{Ii==EE[1,1]+EE[2,2],Qi==EE[1,1]-EE[2,2],Ui==EE[1,2]+EE[2,1],Vi==I(EE[1,2]-EE[2,1])},Flatten[Table[EE[n,m],{n,2},{m,2}]]][[1]]
norm=ν;mat=FullSimplify[ν sisi/.Esol]
ans=FullSimplify[{(*dIds==*)mat[[1,1]]+mat[[2,2]],(*dQds==*)mat[[1,1]]-mat[[2,2]],(*dUds==*)mat[[1,2]]+mat[[2,1]],(*dVds==*)I(mat[[1,2]]-mat[[2,1]])}];
Xi={Ii,Qi,Ui,Vi};
-Table[D[ans[[k]],Xi[[l]]],{k,4},{l,4}]//MatrixForm

Table[α+2β,{α,2},{β,2}]//MatrixForm

(*new definition*)Clear[α];rux={A[n_,m_]->ReA[n,m]+I ImA[n,m]};
(*Table[Sum[I  c/ν ((Subscript[(α^α), μ]) ((E^μ)(E^( *β)))- (Subscript[(α^( *β)), μ])((E^α)(E^( *μ)))),{μ,2}],{α,2},{β,2}]*)
sisi=Simplify[ComplexExpand[Table[Sum[I /ν (A[α,μ]EE[μ,β]-A[β,μ]\[Conjugate]EE[α,μ]),{μ,2}],{α,2},{β,2}]/.{A[2,1]->-A[1,2]}/.rux]];
Esol=Solve[{Ii==EE[1,1]+EE[2,2],Qi==EE[1,1]-EE[2,2],Ui==EE[1,2]+EE[2,1],Vi==I(EE[1,2]-EE[2,1])},Flatten[Table[EE[n,m],{n,2},{m,2}]]][[1]]
norm=ν;mat=FullSimplify[ν sisi/.Esol]
ans=FullSimplify[{(*dIds==*)mat[[1,1]]+mat[[2,2]],(*dQds==*)mat[[1,1]]-mat[[2,2]],(*dUds==*)mat[[1,2]]+mat[[2,1]],(*dVds==*)I(mat[[1,2]]-mat[[2,1]])}];
Xi={Ii,Qi,Ui,Vi};
-Table[D[ans[[k]],Xi[[l]]],{k,4},{l,4}]//MatrixForm

-D[ans[[1]],Ii]
-D[ans[[4]],Vi]
-D[ans[[2]],Ui]
-D[ans[[3]],Vi]

α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]]];

Derivation of response tensor

((∂f)/(∂t))+v ∇f+η e(ΔE+(v/c)×(B+ΔB))(((Subscript[∇, p])f)/(m c))==0;   q=η e (η - sign of charge)                   p==(Sqrt[(γ^2)-1])==(Sqrt[(ϵ^2)/(((m (c^2)))^2)-1]);           (Subscript[f, 0])=(ne/(4π))((Exp[-(Sqrt[1+(p^2)])/θe])/(θe BesselK[2,((θe)^(-1))]));  (∫_0^(+∞))(Subscript[f, 0])4π (p^2)dp=ne
       (Overscript[B0, ->])==B0 (Overscript[(Subscript[e, 3]), ->]);  (Subscript[k, ⊥])=(Subscript[k, y]);  (Subscript[k, x])=0;   px=(Subscript[p, ⊥])Cos[ϕ];  py=(Subscript[p, ⊥])Sin[ϕ]
(Subscript[f, 1])(r,p,t)==Exp[i((Overscript[k, ->]) (Overscript[r, ->])-ω t)](Subscript[f, 0])Φ(p);    v=p c/γ;  (Overscript[k, ->]) (Overscript[r, ->])==(Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ];  ΔE==Exp[i((Overscript[k, ->]) (Overscript[r, ->])-ω t)]ΔE
-i (Subscript[ωf, 1])+((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)(Subscript[f, 1])+η e(ΔE+(v/c)×ΔB)(((Subscript[∇, p])(Subscript[f, 0]))/(m c))+((η e)/γ)(p×(Subscript[B, 0]))(((Subscript[∇, p])(Subscript[f, 1]))/(m c))==0,    since    (p×(Subscript[B, 0]))(Subscript[∇, p])(Subscript[f, 0])==0 – only for isotropic particle distributions
(Subscript[∇, p])(Subscript[f, 0])==-ne((Exp[-(Sqrt[1+(p^2)])/θe])/(((θe)^2)BesselK[2,((θe)^(-1))](Sqrt[1+(p^2)])))p==-((Subscript[f, 0])/(θe(Sqrt[1+(p^2)])))p    =>   ((v/c)×ΔB)(((Subscript[∇, p])(Subscript[f, 0]))/(m c))==0 – only for isotropic particle distributions;   η e ΔE(((Subscript[∇, p])(Subscript[f, 0]))/(m c))==-((Subscript[f, 0])/(θe(Sqrt[1+(p^2)])))((η e)/(m c))(p ΔE)==-((Subscript[f, 0])/(θeγ))((η e)/(m c))(p ΔE)
(Subscript[∇, p])(Subscript[f, 1])==-((ne (Subscript[f, 0]))/(θe(Sqrt[1+(p^2)])))p+(Subscript[f, 1])(((Subscript[∇, p])Φ)/Φ);     (Subscript[∇, p])Φ==(OverHat[(Subscript[p, ⊥])])((∂Φ)/(∂(Subscript[p, ⊥])))+((OverHat[ϕ])/(Subscript[p, ⊥]))((∂Φ)/(∂ϕ))+(OverHat[(Subscript[p, ||])])((∂Φ)/(∂(Subscript[p, ||])))                  
(*1*) (η e)/γ (p×(Subscript[B, 0])) (Subscript[∇, p])(Subscript[f, 1])/(m c) == (η e (Subscript[B, 0]))/(γ m c) (Subscript[f, 1])/Φ (p×(OverHat[(Subscript[p, ||])])) (Subscript[∇, p])Φ == ((η e (Subscript[B, 0]))/(γ m c))(Subscript[f, 1])/Φ ((Subscript[p, ⊥]), 0, (Subscript[p, ||]))×(0, 0, (OverHat[(Subscript[p, ||])])).(Subscript[∇, p])Φ == -((η e (Subscript[B, 0]))/(γ m c)) ((Subscript[f, 1])/Φ) ((Subscript[p, ⊥])/(Subscript[p, ⊥])) ∂Φ/∂ϕ == -((η (Subscript[ω, c]))/γ)(Subscript[f, 0])Exp[i(k r-ω t)]((∂Φ)/(∂ϕ)) =>
-i ω Φ+((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)Φ-(1/(θeγ))((η e)/(m c))(p ΔE)-((η (Subscript[ω, c]))/γ)((∂Φ)/(∂ϕ))==0;      (Subscript[n, ⊥])=(Subscript[k, ⊥])c/ω;   (Subscript[n, ||])=(Subscript[k, ||])c/ω
i ω Φ-((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)Φ+(1/(θeγ))((η e)/(m c))(p ΔE)+((η (Subscript[ω, c]))/γ)((∂Φ)/(∂ϕ))==0;    ((Subscript[ω, c])/ω)==((Subscript[ν, B])/ν)==β
i (ω/(η (Subscript[ω, c])))( γ-(Subscript[n, ||])(Subscript[p, ||])-(Subscript[n, ⊥])(Subscript[p, ⊥])Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==-(1/(θe (Subscript[B, 0])))(p ΔE)
(i/(η β))( γ-(Subscript[n, ||])(Subscript[p, ||])-(Subscript[n, ⊥])(Subscript[p, ⊥])Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==-(1/(θe (Subscript[B, 0])))(p ΔE);   i( a-b Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==-F;
Φ[ϕ]->-(e^(-i (a ϕ+ b Cos[ϕ]))) (∫_(η ∞)^ϕ)(e^(i(a ψ+b Cos[ψ]))) F[ψ]dψ ;        ψ->ϕ-ξ;      (*need to integrate from -∞ to catch all wave history;  special solution dies out with time *)
Φ[ϕ]->-(∫_(-η ∞)^0)(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) F[ξ]dψ     (*simultaneously switch lower-upper boundary & dψ -> –dξ => sign is unchanged*)
Φ[ϕ]->-(∫_0^(-η ∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) F[ξ]dξ;   a=((γ-(Subscript[n, ||])(Subscript[p, ||]))/(η β));    b=(((Subscript[n, ⊥])(Subscript[p, ⊥]))/(η β));  F=(1/((Subscript[θ, e])(Subscript[B, 0])))(p ΔE)==(1/((Subscript[θ, e])(Subscript[B, 0])))((Subscript[E, x])(Subscript[p, ⊥])Cos[ξ-ϕ]-(Subscript[E, y])(Subscript[p, ⊥])Sin[ξ-ϕ]+(Subscript[p, ||])(Subscript[E, z])) 
                 Im[a]∝ η Im[ω]<0 for η= –1;        Im[a]>0 for η=1;   (e^(-i a ϕ-i b Cos[ϕ]))∝(e^(Im[a] ϕ)) – need to die;  ϕ -> +∞ for η= –1;   ϕ->-∞ for η=1 =>  integration from -∞ to ϕ for η=-1 and from +∞ to ϕ for η=1;
(Overscript[j, ->])=η e∫(Subscript[f, 1])(((Overscript[p, ->])c)/γ) d3p;     (Subscript[p, α])=((∂/(∂(Subscript[s, α])))Exp[(Overscript[s, ->])(Overscript[p, ->])])   for    (Overscript[s, ->])->0 ;   (Subscript[p, α])((Overscript[p, ->])(Overscript[E, ->]))=(Subscript[E, β])(((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[(Overscript[s, ->])(Overscript[p, ->])+(Subscript[s', x])(Subscript[p, ⊥])Cos[ξ-ϕ]-(Subscript[s', y])(Subscript[p, ⊥])Sin[ξ-ϕ]+(Subscript[s', z])(Subscript[p, ||])])   for    (Overscript[s, ->])->0
(Subscript[j, α])(k,ω)=-i ω ((((Subscript[ε, αβ])-(Subscript[δ, αβ])))/(4π))(Subscript[E, β])(k,ω)     =>  (Subscript[j, α]) =(Subscript[σ, αβ])(Subscript[E, β])=(((i ω)/c)(Subscript[σ, αβ]))(Subscript[A, β]);
(Subscript[α, αβ])=((η e ω)/c)∫(∫_0^(-η ∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])(c/γ)((Subscript[f, 0])/((Subscript[θ, e])(Subscript[B, 0])))dξ d3p
(Subscript[α, αβ]) = η ((e ω )/(Subscript[B, 0])) (c/c) ne (1/(((θe)^2) BesselK[2,((θe)^(-1))])) ∫(∫_0^(-η ∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])((Exp[-(Sqrt[1 + (p^2)])/θe])/(4π(Sqrt[1+(p^2)])))dξ d3p;    ((e ω )/(Subscript[B, 0])) (c/c) ne= (((e^2) ω me c  )/(me c e (Subscript[B, 0])))ne=(((e^2)ne)/(me c))(ν/(Subscript[ν, B]))=(((Subscript[ω, p])^2)/(4π c))(ν/(Subscript[ν, B]))
(Subscript[α, αβ]) =-η (((e^2)ne  )/(me c))(ν/(Subscript[ν, B]))(1/(((θe)^2) BesselK[2,((θe)^(-1))])) ∫(∫_0^(-η ∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])((Exp[-(Sqrt[1 + (p^2)])/θe])/(4π(Sqrt[1+(p^2)])))dξ d3p
(*3*)Exp[-((Sqrt[1 + (p^2)])/(θe))](e^(-i ((γ-(Subscript[n, ||])(Subscript[p, ||]))/(η β)) ξ+i (((Subscript[n, ⊥])(Subscript[p, ⊥]))/(η β))(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[p, α])(Subscript[((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])==((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[-A(Sqrt[1+(p^2)])-i (Overscript[a, ->])(Overscript[p, ->])]
change variable ξ->β ξ
(Subscript[α, αβ]) = -η (((e^2)ne  )/(me c))(1/(((θe)^2) BesselK[2,((θe)^(-1))])) ∫(∫_0^(-η ∞))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))((Exp[-A(Sqrt[1 + (p^2)])-i (Overscript[a, ->])(Overscript[p, ->])])/(4π(Sqrt[1+(p^2)])))dξ d3p
A->(1/(θe))+((i ξ)/η),ax->((Subscript[n, p])/(β η))-((Cos[β ξ] (Subscript[n, p]))/(β η))+i (Subscript[s, x])+i Cos[β ξ] (Subscript[(s'), x])-i Sin[β ξ] (Subscript[(s'), y]),ay->-((Sin[β ξ] (Subscript[n, p]))/(β η))+i (Subscript[s, y])+i Sin[β ξ] (Subscript[(s'), x])+i Cos[β ξ] (Subscript[(s'), y]),az->-((ξ (Subscript[n, z]))/η)+i (Subscript[s, z])+i (Subscript[(s'), z])
(*4*)(Subscript[α, αβ]) =- η (((e^2)ne  )/(me c))(1/(((θe)^2) BesselK[2,((θe)^(-1))])) ∫(∫_0^(-η ∞))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))((Exp[-A(Sqrt[1 + (p^2)])-i (Overscript[a, ->])(Overscript[p, ->])])/(4π (Sqrt[1+(p^2)])))dξ Sin[χ]dχdϕ'(p^2)dp== η (((e^2)ne  )/(me c))((-1)/(((θe)^2) BesselK[2,((θe)^(-1))])) ((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))(∫_0^(+∞))(∫_0^(+∞))((Exp[-A(Sqrt[1 + (p^2)])])/(Sqrt[1+(p^2)]))((Sin[a p])/(a p))dξ (p^2)dp
(*5*)(∫_0^(+∞))((Exp[-A(Sqrt[1 + (p^2)])])/(Sqrt[1+(p^2)]))((Sin[a p])/(a p))(p^2)dp==(∫_0^(+∞))((Exp[-A(Sqrt[1 + (p^2)])])/(Sqrt[1+(p^2)]))((Sin[a p])/a)((dp^2)/2)==(∫_0^(+∞))Exp[-A(Sqrt[1 + (p^2)])]((Sin[a p])/a)d (Sqrt[1+(p^2)])==(∫_0^(+∞))-((Sin[a p])/(A a))dExp[-A(Sqrt[1 + (p^2)])]=(∫_0^(+∞))((Cos[a p])/A)Exp[-A(Sqrt[1 + (p^2)])]dp=(1/(2A))(∫_(-∞)^(+∞))Exp[-A(Sqrt[1 + (p^2)])+i a p]dp==(1/(2A))(∫_(-∞)^(+∞))Exp[-A Cosh[z]+i a Sinh[z]]Cosh[z]dz==-(1/(2A))(∂/(∂A))(∫_(-∞)^(+∞))Exp[-A Cosh[z]+i a Sinh[z]]dz==
-(1/A)(∂/(∂A))BesselK[0,(Sqrt[(A^2)+(a^2)])]==((BesselK[1,(Sqrt[(a^2)+(A^2)])])/(Sqrt[(a^2)+(A^2)])) (*Gradstein, Ryzhik, p. 372, eq. 3.547.4*)
(Subscript[α, αβ]) = -η (((e^2)ne i  )/(me c))(1/(((θe)^2) BesselK[2,((θe)^(-1))]))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β]))) (∫_0^(-η ∞))((BesselK[1,(Sqrt[(a^2)+(A^2)])])/(Sqrt[(a^2)+(A^2)]))dξ
(Subscript[α, αβ]) = (((e^2)ne i  )/(me c))(1/(((θe)^2) BesselK[2,((θe)^(-1))]))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β]))) (∫_0^∞)((BesselK[1,(Sqrt[(a^2)+(A^2)])])/(Sqrt[(a^2)+(A^2)]))dξ  (*ξ-> –η ξ*);   (1/((θe)^2))-((2 i ξ)/(θe))-((ξ^2)+2 (((-1+Cos[β ξ]))/(β^2)))((Sin[θ])^2))
                                                      δ(γ-(Subscript[γ, 0]))=∫Exp[-(γ/(θe))]γ(Sqrt[(γ^2)-1])K(γ,θe)dθe;  f(γ')=∫δ(γ-(Subscript[γ, 0]))f((Subscript[γ, 0]))Exp[-(γ/(θe))]γ(Sqrt[(γ^2)-1])K(γ,θe)dθe (Subscript[dγ, 0])(*dealing with the tails of small T calculations*) Exp[-A((γ/(θe))+(Sqrt[(γ/(θe))]))]
             
=>   rot(rot((A^ν)))=-(1/(c^2))(((∂^2)(A^ν))/(∂(t^2)))+((4π)/c)(Subscript[α, μν])(A^ν)=-ΔA;     ((k^2)-((ω^2)/(c^2))-((4π)/c)(Subscript[α, μν]))A==0  => ((((k^2)(c^2))/(ω^2))-1-((4π c)/(ω^2))(Subscript[α, μν]))A==0

(*2*)DSolve[{I(a-b Sin[ϕ])Φ[ϕ]+D[Φ[ϕ],ϕ]==-F[ϕ]},Φ[ϕ],ϕ][[1]]/.K[1]->ψ

(mm={{1,0,0},{0,Cos[θ],Sin[θ]},{0,-Sin[θ],Cos[θ]}})//MatrixForm
mm[[3,2]]

DSolve[{I(a-b Sin[ϕ])Φ[ϕ]+D[Φ[ϕ],ϕ]==F[ϕ],Φ[1]==0},Φ[ϕ],ϕ][[1]]/.K[1]->ψ
 {px->(Subscript[p, ⊥])Cos[ϕ],py->(Subscript[p, ⊥])Sin[ϕ]}

(*(*1*){pper,0,ppar}\[Cross]{0,0,ppar}
(*2*)DSolve[{I(a-b Sin[ϕ])Φ[ϕ]+D[Φ[ϕ],ϕ]==-F[ϕ],Φ[1]==0},Φ[ϕ],ϕ][[1]]/.K[1]->ψ
 {px->(Subscript[p, ⊥])Cos[ϕ],py->(Subscript[p, ⊥])Sin[ϕ]}
(*3*) - calculating the modulus*)
xxx=-A γ-I (ax (Subscript[p, x]) +ay (Subscript[p, y])+az (Subscript[p, z]))-TrigExpand[-(γ/(θe))-I ((γ-(Subscript[n, z])(Subscript[p, z]))/(η β)) ξ+I (((Subscript[n, p])(Subscript[p, p]))/(η β))(Cos[ϕ-ξ]-Cos[ϕ])+(Subscript[s, x])(Subscript[p, x])+(Subscript[s, y])(Subscript[p, y])+(Subscript[s, z])(Subscript[p, z])+(Subscript[s', x])(Subscript[p, p])Cos[ξ-ϕ]-(Subscript[s', y])(Subscript[p, p])Sin[ξ-ϕ]+(Subscript[s', z])(Subscript[p, z])]/.{Cos[ϕ]->((Subscript[p, x])/(Subscript[p, p])),Sin[ϕ]->((Subscript[p, y])/(Subscript[p, p]))};vars={γ,(Subscript[p, x]),(Subscript[p, y]),(Subscript[p, z])};
solvA=Expand[Solve[Thread[Table[D[xxx,vars[[i]]],{i,Length[vars]}]==0],{A,ax,ay,az}][[1]]]/.ξ->-η β ξ/.{Cos[β η ξ]->Cos[β ξ],Sin[β η ξ]->η Sin[β ξ]}
Simplify[Collect[Simplify[Expand[(A^2)+(ax^2)+(ay^2)+(az^2)/.solvA]/.{(Subscript[s, x])->0,(Subscript[s, y])->0,(Subscript[s, z])->0,(Subscript[(s'), z])->0,(Subscript[(s'), y])->0,(Subscript[(s'), x])->0}/.{(Subscript[n, z])->Cos[θ],(Subscript[n, p])->Sin[θ]}/.(η^(-2))->1],{β,Sin[θ],Cos[θ]}]==(modul=(1/((θe)^2))-((2 I ξ)/(θe))-((ξ^2)+2 (((-1+Cos[β ξ]))/(β^2)))((Sin[θ])^2))]
(*(*4*)Integrate[(1/(4π))Exp[-I a p Cos[χ]]Sin[χ],{χ,0,π},{ϕ',0,2π}]*)
(*5*)AA=0.2-0.9I;aa=3.+0.1I;
ff[A_,a_]:=Quiet[NIntegrate[Exp[-A Cosh[z]+I a Sinh[z]],{z,-∞,+∞},MaxRecursion->20]];{ff[AA,aa],2BesselK[0,(Sqrt[(AA^2)+(aa^2)])]}
-(1/A)D[BesselK[0,(Sqrt[(A^2)+(a^2)])],A]
stab={(Subscript[s, x]),(Subscript[s, y]),(Subscript[s, z])};sptab={(Subscript[s', x]),(Subscript[s', y]),(Subscript[s', z])};
ann=Simplify[Table[D[((BesselK[1,(Sqrt[(ax^2)+(ay^2)+(az^2)+(A^2)])])/(Sqrt[(ax^2)+(ay^2)+(az^2)+(A^2)]))/.solvA,stab[[i]],sptab[[j]]]/.{(Subscript[s, x])->0,(Subscript[s, y])->0,(Subscript[s, z])->0,(Subscript[(s'), z])->0,(Subscript[(s'), y])->0,(Subscript[(s'), x])->0},{i,3},{j,3}]/.((((1/(θe))-I ξ))^2)+((((Sin[β ξ])^2) (Subscript[n, p]^2))/(β^2))+(((((Subscript[n, p])/(β η))-((Cos[β ξ] (Subscript[n, p]))/(β η))))^2)+(ξ^2) (Subscript[n, z]^2)->(Rex^2),Rex>0]

rel1=2 BesselK[1,x]==x BesselK[2,x]- x BesselK[0,x];
rel2=4 BesselK[2,x]==x BesselK[3,x]- x BesselK[1,x];FullSimplify[rel1]&&FullSimplify[rel2]
sol2=Solve[rel2,BesselK[1,x]][[1]]
sol1=Simplify[Solve[rel1,BesselK[0,x]][[1]]/.sol2]
final=ann/.{BesselK[1,x_]->((-4 BesselK[2,x]+x BesselK[3,x])/x),BesselK[0,x_]->(((8+(x^2)) BesselK[2,x]-2 x BesselK[3,x])/(x^2))};
Simplify[(Rex^2)D[final,BesselK[2,Rex]]]//MatrixForm
Simplify[-(Rex^3)(β^2)D[final,BesselK[3,Rex]]/.{(Subscript[n, z])->Cos[θ],(Subscript[n, p])->Sin[θ],(η^2)->1,(η^(-2))->1,(η^(-1))->η}]/.{((Sin[((β ξ)/2)])^2)->((1-Cos[β ξ])/2),((Sin[((β ξ)/2)])^4)->(((((1-Cos[β ξ])/2)))^2),Sin[2 θ]->2Sin[θ]Cos[θ]}//MatrixForm
Rex->modul

(*for arbitrary isotropic particle distribution*)
((∂f)/(∂t))+v ∇f+η e(ΔE+(v/c)×(B+ΔB))(((Subscript[∇, p])f)/(m c))==0;   q=η e (η - sign of charge)                   p==(Sqrt[(γ^2)-1])==(Sqrt[(ϵ^2)/(((m (c^2)))^2)-1]);           (Subscript[f, 0])=(ne/(4π))((Exp[-(Sqrt[1+(p^2)])/θe])/(θe BesselK[2,((θe)^(-1))]));  (∫_0^(+∞))(Subscript[f, 0])4π (p^2)dp=ne
       (Overscript[B0, ->])==B0 (Overscript[(Subscript[e, 3]), ->]);  (Subscript[k, ⊥])=(Subscript[k, y]);  (Subscript[k, x])=0;   px=(Subscript[p, ⊥])Cos[ϕ];  py=(Subscript[p, ⊥])Sin[ϕ]
(Subscript[f, 1])(r,p,t)==Exp[i((Overscript[k, ->]) (Overscript[r, ->])-ω t)](Subscript[f, 0])Φ(p);    v=p c/γ;  (Overscript[k, ->]) (Overscript[r, ->])==(Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ];  ΔE==Exp[i((Overscript[k, ->]) (Overscript[r, ->])-ω t)]ΔE
-i (Subscript[ωf, 1])+((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)(Subscript[f, 1])+η e(ΔE+(v/c)×ΔB)(((Subscript[∇, p])(Subscript[f, 0]))/(m c))+((η e)/γ)(p×(Subscript[B, 0]))(((Subscript[∇, p])(Subscript[f, 1]))/(m c))==0,    since    (p×(Subscript[B, 0]))(Subscript[∇, p])(Subscript[f, 0])==0 – only for isotropic particle distributions
(Subscript[∇, p])(Subscript[f, 0])==((∂(Subscript[f, 0]))/(∂p))p/p    =>   ((v/c)×ΔB)(((Subscript[∇, p])(Subscript[f, 0]))/(m c))==0 – only for isotropic particle distributions;   η e ΔE(((Subscript[∇, p])(Subscript[f, 0]))/(m c))==∂(Subscript[f, 0])/∂p((η e)/(p m c))(p ΔE)
(Subscript[∇, p])(Subscript[f, 1])==Exp[i((Overscript[k, ->]) (Overscript[r, ->])-ω t)]Φ((∂(Subscript[f, 0]))/(p∂p))p+(Subscript[f, 1])(((Subscript[∇, p])Φ)/Φ);     (Subscript[∇, p])Φ==(OverHat[(Subscript[p, ⊥])])((∂Φ)/(∂(Subscript[p, ⊥])))+((OverHat[ϕ])/(Subscript[p, ⊥]))((∂Φ)/(∂ϕ))+(OverHat[(Subscript[p, ||])])((∂Φ)/(∂(Subscript[p, ||])))                  
(*1*) (η e)/γ (p×(Subscript[B, 0])) (Subscript[∇, p])(Subscript[f, 1])/(m c) == (η e (Subscript[B, 0]))/(γ m c) (Subscript[f, 1])/Φ (p×(OverHat[(Subscript[p, ||])])) (Subscript[∇, p])Φ == ((η e (Subscript[B, 0]))/(γ m c))(Subscript[f, 1])/Φ ((Subscript[p, ⊥]), 0, (Subscript[p, ||]))×(0, 0, (OverHat[(Subscript[p, ||])])).(Subscript[∇, p])Φ == -((η e (Subscript[B, 0]))/(γ m c)) ((Subscript[f, 1])/Φ) ((Subscript[p, ⊥])/(Subscript[p, ⊥])) ∂Φ/∂ϕ == -((η (Subscript[ω, c]))/γ)(Subscript[f, 0])Exp[i(k r-ω t)]((∂Φ)/(∂ϕ)) =>
-i ω Φ+((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)Φ+((∂log((Subscript[f, 0])))/(∂p))((η e)/(p m c))(p ΔE)-((η (Subscript[ω, c]))/γ)((∂Φ)/(∂ϕ))==0;      (Subscript[n, ⊥])=(Subscript[k, ⊥])c/ω;   (Subscript[n, ||])=(Subscript[k, ||])c/ω
i ω Φ-((i((Subscript[k, ||])(Subscript[p, ||])+(Subscript[k, ⊥])(Subscript[p, ⊥])Sin[ϕ])c)/γ)Φ-((∂log((Subscript[f, 0])))/(∂p))((η e)/(p m c))(p ΔE)+((η (Subscript[ω, c]))/γ)((∂Φ)/(∂ϕ))==0;    ((Subscript[ω, c])/ω)==((Subscript[ν, B])/ν)==β    (+(1/(θeγ))((η e)/(m c))(p ΔE))
i (ω/(η (Subscript[ω, c])))( γ-(Subscript[n, ||])(Subscript[p, ||])-(Subscript[n, ⊥])(Subscript[p, ⊥])Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==((∂log((Subscript[f, 0])))/(∂γ))(1/(Subscript[B, 0]))(p ΔE)
(i/(η β))( γ-(Subscript[n, ||])(Subscript[p, ||])-(Subscript[n, ⊥])(Subscript[p, ⊥])Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==((∂log((Subscript[f, 0])))/(∂γ))(1/(Subscript[B, 0]))(p ΔE);   i( a-b Sin[ϕ])Φ+((∂Φ)/(∂ϕ))==-F;
Φ[ϕ]->-(e^(-i (a ϕ+ b Cos[ϕ]))) (∫_(-∞)^ϕ)(e^(i(a ψ+b Cos[ψ]))) F[ψ]dψ ;        ψ->ϕ-ξ;      (*need to integrate from -∞ to catch all wave history;  special solution dies out with time *)
Φ[ϕ]->-(∫_(-∞)^0)(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) F[ψ]dψ
Φ[ϕ]->(∫_0^(+∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) F[ξ]dξ;   a=((γ-(Subscript[n, ||])(Subscript[p, ||]))/(η β));    b=(((Subscript[n, ⊥])(Subscript[p, ⊥]))/(η β));  F=-((∂log((Subscript[f, 0])))/(∂γ))(1/(Subscript[B, 0]))(p ΔE)==((∂log((Subscript[f, 0])))/(∂γ))(1/(Subscript[B, 0]))((Subscript[E, x])(Subscript[p, ⊥])Cos[ξ-ϕ]-(Subscript[E, y])(Subscript[p, ⊥])Sin[ξ-ϕ]+(Subscript[p, ||])(Subscript[E, z]))         
(Overscript[j, ->])=η e∫(Subscript[f, 1])(((Overscript[p, ->])c)/γ) d3p;     (Subscript[p, α])=((∂/(∂(Subscript[s, α])))Exp[(Overscript[s, ->])(Overscript[p, ->])])   for    (Overscript[s, ->])->0 ;   (Subscript[p, α])((Overscript[p, ->])(Overscript[E, ->]))=(Subscript[E, β])(((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[(Overscript[s, ->])(Overscript[p, ->])+(Subscript[s', x])(Subscript[p, ⊥])Cos[ξ-ϕ]-(Subscript[s', y])(Subscript[p, ⊥])Sin[ξ-ϕ]+(Subscript[s', z])(Subscript[p, ||])])   for    (Overscript[s, ->])->0 
(Subscript[j, α])(k,ω)=-i ω ((((Subscript[ε, αβ])-(Subscript[δ, αβ])))/(4π))(Subscript[E, β])(k,ω)     =>  (Subscript[j, α]) =(Subscript[σ, αβ])(Subscript[E, β])=(((i ω)/c)(Subscript[σ, αβ]))(Subscript[A, β]);
(Subscript[α, αβ])=i((η e ω)/c)∫(∫_0^(+∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])(c/γ)(1/(Subscript[B, 0]))(-((∂(Subscript[f, 0]))/(∂γ)))dξ d3p
(Subscript[α, αβ]) = η ((e ω )/(Subscript[B, 0])) (c/c) ∫(∫_0^(+∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])(-((∂(Subscript[f, 0]))/(p ∂p)))dξ d3p;    ((e ω )/(Subscript[B, 0])) (c/c) ne= (((e^2) ω me c  )/(me c e (Subscript[B, 0])))ne=(((e^2)ne)/(me c))(ν/(Subscript[ν, B]))=(((Subscript[ω, p])^2)/(4π c))(ν/(Subscript[ν, B]))
(Subscript[α, αβ]) = η (((e^2)ne  )/(me c))(ν/(Subscript[ν, B]))∫(∫_0^(+∞))(e^(-i a ξ+i b(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[(Subscript[p, α])((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])(-((∂(Subscript[f, 0]))/(p ∂p)))dξ d3p
(*3*)(-((∂(Subscript[f, 0]))/(p ∂p)))(e^(-i ((γ-(Subscript[n, ||])(Subscript[p, ||]))/(η β)) ξ+i (((Subscript[n, ⊥])(Subscript[p, ⊥]))/(η β))(Cos[ϕ-ξ]-Cos[ϕ]))) (Subscript[p, α])(Subscript[((Subscript[p, β]) ((Subscript[ΔE, β])/(ΔE))), ϕ->ϕ-ξ])==(-((∂(Subscript[f, 0]))/(p ∂p)))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[-A(Sqrt[1+(p^2)])-i (Overscript[a, ->])(Overscript[p, ->])]
change variable ξ->β ξ
(Subscript[α, αβ]) = (((e^2)ne  )/(me c)) ∫(∫_0^(+∞))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[-A(Sqrt[1 + (p^2)])-i (Overscript[a, ->])(Overscript[p, ->])](-((∂(Subscript[f, 0]))/(p ∂p)))dξ d3p
A->((i ξ)/η),ax->((Subscript[n, p])/(β η))-((Cos[β ξ] (Subscript[n, p]))/(β η))+i (Subscript[s, x])+i Cos[β ξ] (Subscript[(s'), x])-i Sin[β ξ] (Subscript[(s'), y]),ay->-((Sin[β ξ] (Subscript[n, p]))/(β η))+i (Subscript[s, y])+i Sin[β ξ] (Subscript[(s'), x])+i Cos[β ξ] (Subscript[(s'), y]),az->-((ξ (Subscript[n, z]))/η)+i (Subscript[s, z])+i (Subscript[(s'), z])
(*4*)(Subscript[α, αβ]) = (((e^2)ne  )/(me c)) ∫(∫_0^(+∞))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))Exp[-A(Sqrt[1 + (p^2)])-i (Overscript[a, ->])(Overscript[p, ->])](-((∂(Subscript[f, 0]))/(p ∂p)))dξ Sin[χ]dχdϕ'(p^2)dp== η (((e^2)ne  )/(me c)) ((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β])))(∫_0^(+∞))(∫_0^(+∞))Exp[-A(Sqrt[1 + (p^2)])]((Sin[a p])/(a p))(-((∂(Subscript[f, 0]))/(p ∂p)))dξ (p^2)dp
(*5*)(∫_0^(+∞))Exp[-A(Sqrt[1 + (p^2)])](-((∂(Subscript[f, 0]))/(p ∂p)))((Sin[a p])/(a p))(p^2)dp==(Underoverscript[(1/a)∫, 0, (+∞)])Exp[-A(Sqrt[1 + (p^2)])](-((∂(Subscript[f, 0]))/(∂p)))Sin[a p]dp==(Underoverscript[(1/a)∫, 0, (+∞)])(Subscript[f, 0])d(Exp[-A(Sqrt[1 + (p^2)])]Sin[a p])=
(∫_1^(+∞))(Subscript[f, 0]) Exp[-A γ](-((A p)/(a γ))Sin[a p]+Cos[a p])(γ/P)dγ= Exp[-A (Subscript[γ, 0])](-(A/a)Sin[a (Subscript[p, 0])]+((Subscript[γ, 0])/(Subscript[p, 0]))Cos[a (Subscript[p, 0])]) (*for (Subscript[f, 0])=δ(γ-(Subscript[γ, 0]))*)
 (*Gradstein, Ryzhik, p. 372, eq. 3.547.4*)
(Subscript[α, αβ]) = (((e^2)ne i  )/(me c))(1/(((θe)^2) BesselK[2,((θe)^(-1))]))((∂^2)/(∂(Subscript[s, α])∂(Subscript[s', β]))) (∫_0^(+∞))((BesselK[1,(Sqrt[(a^2)+(A^2)])])/(Sqrt[(a^2)+(A^2)]))dξ
                                                      δ(γ-(Subscript[γ, 0]))=∫Exp[-(γ/(θe))]γ(Sqrt[(γ^2)-1])K(γ,θe)dθe;  f(γ')=∫δ(γ-(Subscript[γ, 0]))f((Subscript[γ, 0]))Exp[-(γ/(θe))]γ(Sqrt[(γ^2)-1])K(γ,θe)dθe (Subscript[dγ, 0])(*dealing with the tails of small T calculations*) Exp[-A((γ/(θe))+(Sqrt[(γ/(θe))]))]
             
=>   rot(rot((A^ν)))=-(1/(c^2))(((∂^2)(A^ν))/(∂(t^2)))+((4π)/c)(Subscript[α, μν])(A^ν)=-ΔA;     ((k^2)-((ω^2)/(c^2))-((4π)/c)(Subscript[α, μν]))A==0  => ((((k^2)(c^2))/(ω^2))-1-((4π c)/(ω^2))(Subscript[α, μν]))A==0

Integrate[((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,+∞},GenerateConditions->False]
NN[θe_]:=NIntegrate[((Exp[-((Sqrt[1+(p^2)])/(θe))])/(θe BesselK[2,((θe)^(-1))]))(p^2),{p,0,+∞}];{NN[1.],NN[2.4]}

Solve[D[((γ/(4TT))-(Sqrt[(γ/(4TT))])),γ]==0][[1]]
TT=.;(γ/(4TT))-(Sqrt[(γ/(4TT))])/.γ->TT
TT=5.4;Plot[Exp[-A((γ/(4TT))-(Sqrt[(γ/(4TT))])+(1/4))]/.A->(10^3.5),{γ,0.9TT,1.1TT},PlotRange->All]

Emissivities Huang & me

signs
Sazonov(1969):  positive B => (Subscript[V, 0])<0; f(RM)>0; h<0 
Melrose (1971):B>0 => (Subscript[ϵ, V])>0;f(RM)>0 (same) 
Shcherbakov (2008): B>0 =>(Subscript[ρ, V])(RM)>0 (same); (Subscript[ρ, Q])<0, but U & V switched =>coincides with Sazonov 
Melrose 2010:

Huang and definitions

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))};
(*Huang*)
Off[NIntegrate::"inumr"];
C1=((3 ee)/(4π (me^3)(cc^5)));C3=((Sqrt[3])/(4π))((ee^3)/(me (cc^2)));n[en_]:=N0 (en^2)Exp[-en/(kb Te)];N0=(1/(2(((kb Te))^3)));
IjcH[x_?NumberQ]:=((1/x)NIntegrate[(z^2)Exp[-z](x/(z^2)) BesselK[(5/3),y],{z,0,+∞},{y,(x/(z^2)),+∞}]);
xm=(ν/(C1 B Sin[θ](((kb Te))^2)));num=5;freq=Table[(10^9.5)(1000^(i/num)),{i,0,num}];Ttab=Table[5 (10^9)(30^(i/num)),{i,0,num}];
expr=(((ee^2)(cc^3)  (me^2))/(Sqrt[3]))N0 kb Te ν IjcH[xm];IjH1=expr/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;
IjH2=expr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
IjH=Flatten[{IjH1,IjH2}]

Sazonov 1969

Clear[en];ν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[θ]))];
IjSs1=Iexpr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;IjSs2=Iexpr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
QjSs1=Qexpr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;QjSs2=Qexpr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
VjSs1=Vexpr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;VjSs2=Vexpr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
IjSs=Flatten[{IjSs1,IjSs2}]
QjSs=Flatten[{QjSs1,QjSs2}]
VjSs=Flatten[{VjSs1,VjSs2}](*negative (Subscript[ϵ, V]) for positive B*)

(QjSs/(IjSs+QjSs))
(VjSs/IjSs)

Comparison

(IjRs/IjGp)
(IjSs/IjGp)
(IjH/IjGp)

Good approximation - surrogate from Melrose, 1971

IjcRs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}];expr=(ns=1)(((Sqrt[3]) B (ee^3) Sin[θ])/(4π (cc^2) me)) IjcRs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];
(*IjcRs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t] Exp[-((γ-1)/(θe))]((((γ-1))^2)/(2((θe)^3))),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}];expr=(ns=1)(((Sqrt[3]) B (ee^3) Sin[θ])/(4π (cc^2) me)) IjcRs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];*)
(*IjcRs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[νΩ BesselK[(5/3),t]Exp[-(γ/(θe))](1/(2((θe)^3))),{γ,0,3000},{t,((νΩ)/(γ^2)),+∞}];expr=(ns=1)(((Sqrt[3]) B (ee^3) Sin[θ])/(4π (cc^2) me)) IjcRs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];*)
(*IjcRs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[BesselK[(5/3),t]Exp[-(γ/(θe))](1/(2((θe)^3))),{γ,0,3000},{t,((νΩ)/(γ^2)),+∞}];expr=(((ee^2) ν)/((Sqrt[3]) cc))IjcRs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];*)
(*IjcRs[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[BesselK[(5/3),t]Exp[-γ](1/(2((θe)^2))),{γ,0,3000θe},{t,((νΩ)/(((γ θe))^2)),+∞}];expr=(((ee^2) ν)/((Sqrt[3]) cc))IjcRs[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];*)
IjRs1=expr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;IjRs2=expr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
IjRs=Flatten[{IjRs1,IjRs2}]

Melrose 1971 exact

IjcRp[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/γ)BesselK[(5/3),t](((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])),{γ,1,3000},{t,((νΩ)/γ),+∞}];
IjcRp[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t](((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1])),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}];
(*IjcRp[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]γ Exp[-((γ-1)/(θe))]((γ^2)/(2((θe)^3))),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}]*)
expr=(ns=1)(((Sqrt[3])(ee^3))/(16(π^2)(cc^2)me))B Sin[θ]IjcRp[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))];
IjRp1=expr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;IjRp2=expr/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
IjRp=Flatten[{IjRp1,IjRp2}]

Gammie - very good approximation!

(*Gammie*)νc=((ee B)/(2π me cc));νs=(2/9)νc ((θe)^2)Sin[θ];X=(ν/(νs));
IjG1=(ns=1)(((Sqrt[2])π (ee^2)ν)/(6((θe)^2)cc))Exp[-(((ν/νs))^(1/3))]/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;
IjG2=(ns=1)(((Sqrt[2])π (ee^2)ν)/(6((θe)^2)cc))Exp[-(((ν/νs))^(1/3))]/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;IjG=Flatten[{IjG1,IjG2}];
expr=(ns=1)(((Sqrt[2])π (ee^2)νs)/(3BesselK[2,((θe)^(-1))]cc))((((X^(1/2))+(2^(11/12))(X^(1/6))))^2)Exp[-(X^(1/3))];IjGp1=expr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq;
IjGp2=expr/.θe->((kb Te)/(me (cc^2)))/.ruru/.{θ->1.,B->20,Te->#,ν->3 (10^11)}&/@Ttab;
IjGp=Flatten[{IjGp1,IjGp2}]

IjcRp[θe_?NumberQ,νΩ_?NumberQ]:=NIntegrate[((νΩ)/(γ^2))BesselK[(5/3),t]γ Exp[-((γ-1)/(θe))]((γ^2)/(2((θe)^3))),{γ,1,3000},{t,((νΩ)/(γ^2)),+∞}]
IjRp=(*2 π*)(ns=1)((Sqrt[3])/(16(π^2)(cc^2)))((ee^3)/me)2B Sin[θ]IjcRp[((kb Te)/(me (cc^2))),((2me cc 2π ν)/(3ee B Sin[θ]))]/.ruru/.{θ->1.,B->20,Te->2. (10^10),ν->#}&/@freq

Distributions

Integrate[((Exp[-(γ/(θe))])/(θe BesselK[2,((θe)^(-1))]))γ(Sqrt[(γ^2)-1]),{γ,1,+∞},GenerateConditions->False]
Integrate[Exp[-((γ-1)/(θe))]((((γ-1))^2)/(2((θe)^3))),{γ,1,+∞},GenerateConditions->False]

Transfer equations

(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

χexp=2π s; χrul={χ[s]->χexp,χ'[s]->D[χexp,s]};acc=12;
S={Ij[s],Qj[s],Uj[s],Vj[s]};ϵ={4,2,0,1};
sys=Flatten[{Thread[R.D[S,s]==R.ϵ]/.χrul,Ij[0]==0,Qj[0]==0,Uj[0]==0,Vj[0]==0}];
newsys=Flatten[{Thread[D[R.S,s]==R.ϵ]/.χrul,Ij[0]==0,Qj[0]==0,Uj[0]==0,Vj[0]==0}];

sol=NDSolve[sys,S,{s,0,1}][[1]];
newsol=NDSolve[newsys,S,{s,0,1},AccuracyGoal->acc,PrecisionGoal->acc][[1]];
S/.sol/.s->1.
Chop[R.S/.sol/.χrul/.s->1.]
Chop[S/.newsol/.s->1.]

(mat={{αI,αQ,0,αV},{αQ,αI,ρV,0},{0,-ρV,αI,ρQ},{αV,0,-ρQ,αI}})//MatrixForm

FullSimplify[Transpose[R].mat.R]//MatrixForm

Full compilation of transport 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[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]]

CPU hours

(*π simulations of certain spin etc.*)
((70(*hr*))/(8600(*M*)))((20000.(*M*))/(24(*hr/day*)))

Exp[0.4]
0.60 (0.96^(#))&/@Range[0,37]
(*Stage 1 - searching for the line in (heat, rhonor) plane for spectrum fitting*)Print[tdes=31(*angles*) 38.(*heat values*)200(*s per pair (θ,heat)*) (1/(3600(*h per s*)))," h = desktop time"]
ratio=((8 2.2)/(4 2.83))(*1 node speed Odyssey/desktop*)
Print[5(*spins*)(tdes/ratio)(1/(31(*nodes*))),"h = full search of (heat, rhonor) pairs for all spins for (Subscript[F, ν])"]

hr=(1/(3600(*h/s*)))40(*s*) 11(*iterations per value*) 38.(*heat values*)(*56 cores*)
CPUhr=56hr
tothr=(1/(32(*8-core nodes*)3600(*s/hr*)24(*hr/day*)))40(*s*) 11(*iterations per value*) 38.(*heat values*)17(*θ values*)5(*values of spin*)10(*timeshots per spin*)

Times

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))};
3000 (rg/(2cc))(1/(60 60))/.ruru (*hours*)
Print["Period = ",150 (rg/(2cc))(1/(60 ))/.ruru ," min"](*hours*)
(*2π (Sqrt[((r^3)/(G M))])==2π ((((r/M)))^(3/2))(rg/(2cc))*)
rr=25;
Print["Period at r=",rr,"M:  T=",Round[2.π (((rr))^(3/2))/.ruru],"M"]
Print[dd=(rg/(2cc 60))(*M*)/.ruru, " minutes per M"]
17/dd