diff --git a/Kernel/Hyperboloidal.wl b/Kernel/Hyperboloidal.wl new file mode 100644 index 0000000..36072b4 --- /dev/null +++ b/Kernel/Hyperboloidal.wl @@ -0,0 +1,298 @@ +(* ::Package:: *) + +(* ::Chapter:: *) +(*ReggeWheelerHyperboloidal*) + + +(* ::Subsubsection:: *) +(*For details on solving the Regge-Wheeler-Zerilli equations using hyperboloidal compactification, see arXiv : gr - qc/2202.01794 and arXiv : gr - qc/2411.14976 .*) + + +(* ::Section:: *) +(*Create Package*) + + +(* ::Subsection:: *) +(*Begin Package*) + + +BeginPackage["ReggeWheeler`Hyperboloidal`", + {"ReggeWheeler`ReggeWheelerSource`", + "ReggeWheeler`ReggeWheelerRadial`", + "ReggeWheeler`ConvolveSource`", + "KerrGeodesics`KerrGeoOrbit`", + "KerrGeodesics`OrbitalFrequencies`", + "SpinWeightedSpheroidalHarmonics`"} +]; + + +(* ::Subsection:: *) +(*Begin Private Section*) + + +Begin["`Private`"]; + + +(* ::Section:: *) +(*Utility functions*) + + +(* Schwarzschild f(r) *) + f[r_,M_]:= 1-2 M/r; + +(* Fns. defined for ease of repeated use *) + L[l_] := l(l+1); + \[Lambda][l_]:=((l+2)(l-1))/2; + \[CapitalLambda][r_,l_,M_]:=\[Lambda][l] + (3M)/r; + p[r_]:=(8*\[Pi])/r^2; + q[r_,l_,M_]:= f[r,M]^2/((\[Lambda][l]+1)\[CapitalLambda][r,l,M]); + Y\[Phi][l_,m_,\[Theta]_]:= Conjugate[(((D[SphericalHarmonicY[l,m,\[Theta],\[Phi]],{\[Phi],2}])/.\[Phi]->0)+Sin[\[Theta]]Cos[\[Theta]]((D[SphericalHarmonicY[l,m,x,0],x])/.x->\[Theta]) + +L[l]/2 Sin[\[Theta]]^2 SphericalHarmonicY[l,m,\[Theta],0])]; + X\[Theta][l_,m_,\[Theta]_]:= Conjugate[Sin[\[Theta]]D[SphericalHarmonicY[l,m,x,0],x]/.x->\[Theta]]; + +(* azimuthal orbital frequency *) + \[CapitalOmega][r_,M_] := 1/r Sqrt[M/r]; + +(* angular frequency *) + \[Omega][\[CapitalOmega]_,m_]:= m \[CapitalOmega]; + +(* orbital angular momentum per unit mass *) + L0[r_,M_]:=r Sqrt[M/(r-3M)]; + +(* orbital energy per unit mass *) + \[ScriptCapitalE][r_,M_]:=(r-2M)/Sqrt[r(r-3M)]; + +(* Even/Odd parity potentials *) + VeffEven[r_,l_,M_]:=f[r,M]/(r^2 \[CapitalLambda][r,l,M]^2) (2\[Lambda][l]^2 (\[Lambda][l]+1+(3M)/r)+18 M^2/r^2 (\[Lambda][l]+M/r)); + VeffOdd[r_,l_,M_]:= f[r,M]/r^2 (L[l] - 6 M/r); + +(* RWZ gauge source terms *) + SourceTerm1Even[r_,M_,l_,m_,\[Theta]_]:=(( p[r]q[r,l,M] \[ScriptCapitalE][r,M])/(r f[r,M]\[CapitalLambda][r,l,M]) (L0[r,M]^2/\[ScriptCapitalE][r,M]^2 f[r,M]^2 \[CapitalLambda][r,l,M]-(\[Lambda][l](\[Lambda][l]+1)r^2+6\[Lambda][l]M r +15M^2)) + Conjugate[SphericalHarmonicY[l,m,\[Theta],0]]-(4p[r]L0[r,M]^2 f[r,M]^2)/(r \[ScriptCapitalE][r,M]) (l-2)!/(l+2)! Y\[Phi][l,m,\[Theta]]); + SourceTerm2Even[r_,M_,l_,m_,\[Theta]_]:=(p[r] q[r,l,M]r^2 \[ScriptCapitalE][r,M])Conjugate[SphericalHarmonicY[l,m,\[Theta],0]]; + SourceTerm1Odd[r_,M_,l_,m_,\[Theta]_]:= -((2 p[r]f[r,M]L0[r,M])/(\[Lambda][l]*L[l]))X\[Theta][l,m,\[Theta]]; + SourceTerm2Odd[r_,M_,l_,m_,\[Theta]_]:= (2 p[r] r f[r,M]^2 L0[r,M])/(\[Lambda][l]*L[l])X\[Theta][l,m,\[Theta]]; + +(* Hyperboloidal height fn. *) + H[\[Sigma]_] :=1/2 (Log[1-\[Sigma]]-1/\[Sigma]+Log[\[Sigma]]); + +(* Hyperboloidal scaling factor *) + Z[\[Sigma]_,\[Xi]_]:= 1/2 E^(\[Xi] H[\[Sigma]]); + +(* Hyperboloidal re-scaling factor *) + \[ScriptCapitalF][r_,M_,\[Xi]_]:= f[r,M] 1/(2r^2) E^(\[Xi] H[2/r]); + +(* Coefficients of hyperboloidal master eqn. operator *) + \[Alpha]2[\[Sigma]_] := \[Sigma]^2 (1-\[Sigma]); + \[Alpha]1[\[Sigma]_,\[Xi]_]:= \[Sigma](2-3\[Sigma])+\[Xi](1-2\[Sigma]^2); + (* Distinct coefficients for the even/odd parity equations *) + \[Alpha]0Odd[\[Sigma]_,l_,M_,\[Xi]_]:= -(\[Xi]^2 (1+\[Sigma])+2\[Xi] \[Sigma]+(4M^2 VeffOdd[2/\[Sigma],l,M])/((1-\[Sigma])\[Sigma]^2)); + \[Alpha]0Even[\[Sigma]_,l_,M_,\[Xi]_]:= -(\[Xi]^2 (1+\[Sigma])+2\[Xi] \[Sigma]+(4M^2 VeffEven[2/\[Sigma],l,M])/((1-\[Sigma])\[Sigma]^2)); + +(* Analytic mesh refinement (AnMR) \[Kappa] and coordinate transformation *) + \[Kappa][r0_]:=1/2 Log[r0]; + AnMRTransform[\[Chi]_,r0_]=-1(1-(2Sinh[\[Kappa][r0](1+ \[Chi])])/Sinh[2\[Kappa][r0]]); + + +(* ::Section:: *) +(*Coordinate Transformation*) + + +DomainMapping[r0_,x_,X_]:= + Module[{\[Sigma]p, \[Sigma]grid1, \[Sigma]grid2, M, AB1, AB2, map1, map2, InvMap1, InvMap2, \[Chi]ofX, \[Chi]of\[Sigma], \[Sigma]of\[Chi], a, b}, + + (* Initial setup *) + M = 1; + \[Sigma]p =(2M)/r0; + + \[Sigma]grid1 = {0,\[Sigma]p}; + \[Sigma]grid2 = {\[Sigma]p,1}; + + (* Solving for the transform coefficients *); + AB1 = Solve[{a*\[Sigma]grid1[[1]]+b==-1,a*\[Sigma]grid1[[2]]+b==1}][[1]]; + AB2 = Solve[{a*\[Sigma]grid2[[1]]+b==-1,a*\[Sigma]grid2[[2]]+b==1}][[1]]; + + (* Map to hyp. coords. *); + map1 = a*x+b/.AB1; + map2 = a*x+b/.AB2; + + (* Map to Sch. coords *); + InvMap1 = (X-b)/a/.AB1; + InvMap2 = (X-b)/a/.AB2; + + \[Chi]ofX = \[Chi]/.Solve[AnMRTransform[\[Chi],r0]==X,\[Chi]][[1]]/.C[1]->0//Quiet; + \[Chi]of\[Sigma] = \[Chi]ofX/.X->map2; + \[Sigma]of\[Chi] = InvMap2/.X->AnMRTransform[X,r0]; + + Return[{map1, map2, InvMap1, InvMap2, \[Chi]ofX, \[Chi]of\[Sigma], \[Sigma]of\[Chi]}] +]; + + +(* ::Section:: *) +(*Spectral ODE Solver*) + + +Options[HyperboloidalSolver]={"GridPoints"->32}; + + +HyperboloidalSolver[r0_, l_, m_, Xgrid_, opts:OptionsPattern[]]:=Module[ + {npts, M, \[Theta], \[Xi], \[Sigma]p, prec, map1, map2, InvMap1, InvMap2, AnMRMap\[Chi]X, AnMRMap\[Chi]\[Sigma], AnMRMap\[Sigma]\[Chi], d\[Chi]dX, d2\[Chi]dX2, x, X, \[Phi], S1, S2, ansatz, Dansatz, D2ansatz, Dmap1, Dmap2, + cs1, cs2, cs, DH, A, B, ansatz1D1, ansatz1D2, ansatz2D1, ansatz2D2, \[Alpha]21, + \[Alpha]22, \[Alpha]11, \[Alpha]12, \[Alpha]01, \[Alpha]02, BCs, BCsRHS, D\[Alpha]2, ODEs, juncs, juncs1, + juncs2, fill, juncsRHS, dom1, dom2, Mat, LARHS, + sols, sols2, csNew, sol1, sol2, map1New, map2New, y, sol1New, sol2New, poly1, poly2}, + (*Module internal number of points on grid (maybe unnecessary)*) + npts = OptionValue["GridPoints"]; + + (* Initial setup *) + M = 1; + \[Theta] = \[Pi]/2; + \[Xi] = -I \[Omega][\[CapitalOmega][r0,M],m] 4M; + (* Source radial position in hyp coords *) + \[Sigma]p = 2/r0; + (* Setting working precision *) + prec = Precision[r0]; + + (* Coordinate mappings based off source position *) + {map1, map2, InvMap1} = DomainMapping[r0,x,X][[1;;3]]; + {InvMap2[X_],AnMRMap\[Chi]X[X_], AnMRMap\[Chi]\[Sigma][x_], AnMRMap\[Sigma]\[Chi][X_]} = DomainMapping[r0,x,X][[4;;-1]]; + + (* Source terms *) + {S1,S2} = If[EvenQ[l+m], + {SourceTerm1Even[r0,M,l,m,\[Theta]],SourceTerm2Even[r0,M,l,m,\[Theta]]}, + {SourceTerm1Odd[r0,M,l,m,\[Theta]],SourceTerm2Odd[r0,M,l,m,\[Theta]]}]; + + (* Chebyshev polynomial as ansatz, using number of points on grid to determine length *) + ansatz = Table[ChebyshevT[i,X],{i,0,Length[Xgrid]+1}]; + + (* Derivatives of ansatz and mappings *) + {Dansatz, D2ansatz, Dmap1, Dmap2} = {D[ansatz,X], D[ansatz,{X,2}], D[map1,x], D[map2,x]}; + {d\[Chi]dX[X_], d2\[Chi]dX2[X_]} = {D[AnMRMap\[Chi]X[X],X], D[AnMRMap\[Chi]X[X],{X,2}]}; + + (* Initialising table of weight coefficients *) + {cs1,cs2} = {Table[Subscript[c, i, 1],{i,0,Length[Xgrid]+1}],Table[Subscript[c, i, 2],{i,0,Length[Xgrid]+1}]}; + cs = Join[cs1,cs2]; + + (*Derivative of height fn for convenience*) + DH =(D[H[x],x])/.x->\[Sigma]p; + + (* Hyperboloidal junction condition coefficients *) + {A,B} = { (2E^(-\[Xi] H[\[Sigma]p]))/(1-\[Sigma]p) (2 S1+\[Sigma]p^2/(1-\[Sigma]p) (1-\[Xi](1-\[Sigma]p)(DH))S2),-((2 \[Sigma]p^2 E^(-\[Xi] H[\[Sigma]p]))/(1-\[Sigma]p))S2}; + + (* Transforming ansatz derivatives *) + {ansatz1D1, ansatz2D1, ansatz1D2, ansatz2D2} = {Dmap1*Dansatz, Dmap2*(d\[Chi]dX[AnMRTransform[X,r0]])*Dansatz, Dmap1^2*D2ansatz, Dmap2^2*((d2\[Chi]dX2[AnMRTransform[X,r0]])*Dansatz+(d\[Chi]dX[AnMRTransform[X,r0]]^2)*D2ansatz)}; + + (* Obtaining master fn operator coefficients on Chebyshev-Gauss-Lobatto grid *) + {\[Alpha]21, \[Alpha]22, \[Alpha]11, \[Alpha]12, \[Alpha]01, \[Alpha]02} = If[EvenQ[l+m], + {\[Alpha]2[InvMap1], \[Alpha]2[InvMap2[AnMRTransform[X,r0]]], \[Alpha]1[InvMap1,\[Xi]], \[Alpha]1[InvMap2[AnMRTransform[X,r0]],\[Xi]], \[Alpha]0Even[InvMap1,l,M,\[Xi]], \[Alpha]0Even[InvMap2[AnMRTransform[X,r0]],l,M,\[Xi]]}, + {\[Alpha]2[InvMap1], \[Alpha]2[InvMap2[AnMRTransform[X,r0]]], \[Alpha]1[InvMap1,\[Xi]], \[Alpha]1[InvMap2[AnMRTransform[X,r0]],\[Xi]], \[Alpha]0Odd[InvMap1,l,M,\[Xi]], \[Alpha]0Odd[InvMap2[AnMRTransform[X,r0]],l,M,\[Xi]]}]; + + (* Defining boundary conditions, junction conditions and the ODE at every other grid point *) + {BCs, ODEs, juncs1} = {{((\[Alpha]11 ansatz1D1 + \[Alpha]01 ansatz)/.X -> -1),((\[Alpha]12 ansatz2D1 + \[Alpha]02 ansatz)/.X -> 1)}, + {((\[Alpha]21 ansatz1D2 + \[Alpha]11 ansatz1D1 + \[Alpha]01 ansatz)/.X -> Xgrid),((\[Alpha]22 ansatz2D2 + \[Alpha]12 ansatz2D1 + \[Alpha]02 ansatz)/.X -> Xgrid)}, + {((cs2 . ansatz/.X -> -1) - (cs1 . ansatz/.X -> 1)),((cs2 . ansatz2D1/.X -> -1) - (cs1 . ansatz1D1/.X -> 1))}}; + + (* Obtaining values for junction conditions *) + juncs2 = Table[Coefficient[juncs1[[j]],cs[[i]]], {j, Length[juncs1]}, {i, Length[cs]}]; + + (* Placeholder column of zeroes for matrix construction (replaced by junction condition vectors in full matrix ) *) + fill = ConstantArray[0, Length[ansatz]]; + + (* Constructing matrices to be inverted for each domain *); + dom1 = MapThread[Append, {MapThread[Prepend,{ODEs[[1]],BCs[[1]]}],fill}]; + dom2 = MapThread[Append,{MapThread[Prepend,{ODEs[[2]],fill}],BCs[[2]]}]; + + (* Constructing overall block diagonal matrix to be inverted *) + Mat = ArrayFlatten[{{dom1,0},{0,dom2}}]; + (*Setting junction conditions*) + {Mat[[;;,Length[Mat]/2]], Mat[[;;,Length[Mat]/2+1]]} = juncs2; + + (*Initialising values for RHS of OVERALL equation Mx = b.*); + BCsRHS = {0,0}; + D\[Alpha]2 = (D[\[Alpha]2[x],x])/.x->\[Sigma]p; + juncsRHS = {B/\[Alpha]2[\[Sigma]p ],(A/\[Alpha]2[\[Sigma]p ]+B (D\[Alpha]2 -\[Alpha]1[\[Sigma]p,\[Xi] ])/\[Alpha]2[\[Sigma]p ]^2)}; + + (* Constructing RHS vector *); + LARHS = Join[{BCsRHS[[1]]},ConstantArray[0,Length[Xgrid]],juncsRHS,ConstantArray[0,Length[Xgrid]],{BCsRHS[[2]]}]; + + (*Inverting matrix to obtain weight coefficients *); + sols = LARHS . Inverse[Mat]; + + + (* Creating table of replacement rules for weight coefficients *); + csNew = Table[cs[[i]]->sols[[i]],{i,1,Length[cs]}]; + + (* Creating final polynomial *); + sol1[x_]:= (cs1 . ansatz/.csNew)/.X->x; + sol2[x_] := (cs2 . ansatz/.csNew)/.X->x; + map1New[y_]:= map1/.x->y; + map2New[y_]:= AnMRMap\[Chi]\[Sigma][y]; + sol1New[x_]:= sol1[map1New[x]]; + sol2New[x_]:= sol2[map2New[x]]; + poly1 = Function[\[Sigma],Which[\[Sigma]<\[Sigma]p,sol1New[\[Sigma]],\[Sigma]>\[Sigma]p,sol2New[\[Sigma]]]]; + poly2 = Function[r, Which[rr0,Z[2/r,\[Xi]]poly1[2/r]]]; + Return[{poly1,poly2}] +] + + +(* ::Section:: *) +(*Overall Module*) + + +Options[ReggeWheelerHyperboloidal]={"GridPoints" -> 32}; + + +ReggeWheelerHyperboloidal[s_Integer, l_Integer, m_Integer, n_Integer, orbit_KerrGeoOrbitFunction, opts:OptionsPattern[]]:= + Module[{r0, M, w, grid,Xgrid,prec,npts, S, R, R\[Sigma], Z}, + + (* Initial setup *) + M = 1; + r0 = orbit["p"]; + + (* Initialising number of grid points *) + npts = OptionValue["GridPoints"]; + + (* Setting working precision *) + prec = Precision[r0]; + + (* Initialising Chebyshev-Gauss-Lobatto grid *) + Xgrid = N[Cos[(Range[0,npts]\[Pi])/npts][[2;;-2]],prec]; + + (* Output *) + R = HyperboloidalSolver[r0, l, m, Xgrid, + "GridPoints" -> npts + ][[2]]; + + R\[Sigma] = HyperboloidalSolver[r0, l, m, Xgrid, + "GridPoints" -> npts + ][[1]]; + + S = SpinWeightedSpheroidalHarmonicS[s, l, m, 0]; + Z = <| "\[ScriptCapitalI]" -> R\[Sigma][0], "\[ScriptCapitalH]" -> R\[Sigma][1] |>; + w = \[Omega][\[CapitalOmega][r0,M],m]; + + assoc = <| "s" -> 2, + "l" -> l, + "m" -> m, + "\[Omega]" -> w, + "Type" -> {"PointParticleCircular","Orbital Radius"->ToString[r0] <>"M"}, + "RadialFunction" -> R, + "AngularFunction" -> S, + "Amplitudes" -> Z, + "Method" -> {"Hyperboloidal", "GridPoints"->npts} + |>; + + ReggeWheeler`ReggeWheelerMode`ReggeWheelerMode[assoc] +] + + +(* ::Section:: *) +(*End Package*) + + +(* ::Subsection:: *) +(*End*) + + +End[]; +EndPackage[]; diff --git a/Kernel/ReggeWheeler.m b/Kernel/ReggeWheeler.m index 6eb4f25..54068be 100644 --- a/Kernel/ReggeWheeler.m +++ b/Kernel/ReggeWheeler.m @@ -10,6 +10,7 @@ ]; Get["ReggeWheeler`NumericalIntegration`"]; +Get["ReggeWheeler`Hyperboloidal`"]; Get["ReggeWheeler`ReggeWheelerRadial`"]; Get["ReggeWheeler`ReggeWheelerSource`"]; Get["ReggeWheeler`ReggeWheelerMode`"]; diff --git a/Kernel/ReggeWheelerMode.m b/Kernel/ReggeWheelerMode.m index cd5bfd0..8a38203 100644 --- a/Kernel/ReggeWheelerMode.m +++ b/Kernel/ReggeWheelerMode.m @@ -4,11 +4,11 @@ (*ReggeWheelerMode*) -(* ::Section::Closed:: *) +(* ::Section:: *) (*Create Package*) -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*BeginPackage*) @@ -38,11 +38,15 @@ "ReggeWheelerMode representing a solution to the Regge-Wheeler equation with a point particle source."; -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Error Messages*) ReggeWheelerPointParticleMode::nospin = "Regge-Wheeler perturbations are only for Schwarzschild black holes but spin `1` is not zero." +ReggeWheelerPointParticleMode::eccentricity = "The hyperboloidal package does not currently accept eccentric orbits. Please set eccentricity ('e') to zero." +ReggeWheelerPointParticleMode::spin2field = "The hyperboloidal package currently only works for spin = 2 fields, but the fluxes and radial fns. are correct for spin = -2. Please set field spin ('s') to two." +ReggeWheelerPointParticleMode::inclination = "The hyperboloidal package currently only works for orbits in the equatorial plane. Please set orbital inclination ('x') to one." +ReggeWheelerPointParticleMode::eccentricitymode = "The hyperboloidal package currently only works for circular orbit modes ('m'). Please set the eccentricity mode ('n') to zero." (* ::Subsection::Closed:: *) @@ -52,7 +56,7 @@ Begin["`Private`"]; -(* ::Section::Closed:: *) +(* ::Section:: *) (*ReggeWheelerPointParticleMode*) @@ -60,15 +64,64 @@ {"ArgumentsPattern" -> {_, _, _, _, _, OptionsPattern[]}}; -Options[ReggeWheelerPointParticleMode] = {}; +Options[ReggeWheelerPointParticleMode] = {Method->Automatic}; ReggeWheelerPointParticleMode[s_Integer, l_Integer, m_Integer, n_Integer, orbit_KerrGeoOrbitFunction, opts:OptionsPattern[]] /; AllTrue[orbit["Frequencies"], InexactNumberQ] := - Module[{source, assoc, radialopts, R, S, \[Omega], \[CapitalOmega]r, \[CapitalOmega]\[Phi], \[CapitalOmega]\[Theta], Z}, + Module[{source, assoc, radialopts, R, S, \[Omega], \[CapitalOmega]r, \[CapitalOmega]\[Phi], \[CapitalOmega]\[Theta], Z,hypopts}, If[orbit["a"] != 0, Message[ReggeWheelerPointParticleMode::nospin, orbit["a"]]; Return[$Failed]; ]; + + + Switch[OptionValue["Method"], + "Hyperboloidal", + If[orbit["e"] != 0, + Message[ReggeWheelerPointParticleMode::eccentricity]; + Return[$Failed]; + ]; + + If[s != 2, + Message[ReggeWheelerPointParticleMode::spin2field]; + Return[$Failed]; + ]; + + If[orbit["Inclination"] != 1, + Message[ReggeWheelerPointParticleMode::inclination]; + Return[$Failed]; + ]; + + If[n != 0, + Message[ReggeWheelerPointParticleMode::eccentricitymode]; + Return[$Failed]; + ]; + Return[ReggeWheeler`Hyperboloidal`Private`ReggeWheelerHyperboloidal[s,l,m,n,orbit]], + + {"Hyperboloidal",OptionsPattern[ReggeWheeler`Hyperboloidal`Private`ReggeWheelerHyperboloidal]}, + If[orbit["e"] != 0, + Message[ReggeWheelerPointParticleMode::eccentricity]; + Return[$Failed]; + ]; + + If[s != 2, + Message[ReggeWheelerPointParticleMode::spin2field]; + Return[$Failed]; + ]; + + If[orbit["Inclination"] != 1, + Message[ReggeWheelerPointParticleMode::inclination]; + Return[$Failed]; + ]; + + If[n != 0, + Message[ReggeWheelerPointParticleMode::eccentricitymode]; + Return[$Failed]; + ]; + hypopts = Sequence@@FilterRules[{OptionValue["Method"][[2;;]]}, Options[ReggeWheeler`Hyperboloidal`Private`ReggeWheelerHyperboloidal]]; + Return[ReggeWheeler`Hyperboloidal`Private`ReggeWheelerHyperboloidal[s,l,m,n,orbit,hypopts]] + ]; + (*{\[CapitalOmega]r, \[CapitalOmega]\[Theta], \[CapitalOmega]\[Phi]} = orbit["Frequencies"];*) (*This gives Mino frequencies, need BL frequencies*) {\[CapitalOmega]r, \[CapitalOmega]\[Theta], \[CapitalOmega]\[Phi]} = Values[KerrGeoFrequencies[orbit["a"], orbit["p"], orbit["e"], orbit["Inclination"]]]; @@ -93,14 +146,15 @@ "Type" -> {"PointParticleCircular", "Radius" -> orbit["p"]}, "RadialFunctions" -> R, "AngularFunction" -> S, - "Amplitudes" -> Z + "Amplitudes" -> Z, + "Method" -> R["In"]["Method"] |>; ReggeWheelerMode[assoc] ] -(* ::Section::Closed:: *) +(* ::Section:: *) (*ReggeWheelerMode*) @@ -147,7 +201,10 @@ ReggeWheelerMode[assoc_][string_] := assoc[string]; -(* ::Section::Closed:: *) +Keys[m_ReggeWheelerMode] ^:= Keys[m[[1]]] + + +(* ::Section:: *) (*Fluxes*) @@ -156,14 +213,19 @@ EnergyFlux[mode_ReggeWheelerMode] := - Module[{l, m, \[Omega], Z, FluxInf, FluxH}, + Module[{l, m, \[Omega], Z, FluxInf, FluxH,\[Xi]}, l = mode["l"]; m = mode["m"]; \[Omega] = mode["\[Omega]"]; Z = mode["Amplitudes"]; - - FluxInf = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1))*Abs[\[Omega]*Z["\[ScriptCapitalI]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2))*Abs[\[Omega]*Z["\[ScriptCapitalI]"]]^2/(16*Pi)]; - FluxH = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1))*Abs[\[Omega]*Z["\[ScriptCapitalH]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2))*Abs[\[Omega]*Z["\[ScriptCapitalH]"]]^2/(16*Pi)]; + + If[mode["Method"][[1]] == "Hyperboloidal", + FluxInf = ((l+2)!/(l-2)!)1/(256\[Pi]*16)Abs[-I \[Omega] 4 Z[[1]]]^2; + FluxH = ((l+2)!/(l-2)!)1/(256\[Pi]*16)Abs[-I \[Omega] 4 Z[[2]]]^2; + , + FluxInf = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1))*Abs[\[Omega]*Z["\[ScriptCapitalI]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2))*Abs[\[Omega]*Z["\[ScriptCapitalI]"]]^2/(16*Pi)]; + FluxH = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1))*Abs[\[Omega]*Z["\[ScriptCapitalH]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2))*Abs[\[Omega]*Z["\[ScriptCapitalH]"]]^2/(16*Pi)]; + ]; <| "\[ScriptCapitalI]" -> FluxInf, "\[ScriptCapitalH]" -> FluxH |> ]; @@ -174,14 +236,19 @@ AngularMomentumFlux[mode_ReggeWheelerMode] := - Module[{l, m, \[Omega], Z, FluxInf, FluxH}, + Module[{l, m, \[Omega], Z, FluxInf, FluxH,\[Xi]}, l = mode["l"]; m = mode["m"]; \[Omega] = mode["\[Omega]"]; Z = mode["Amplitudes"]; - FluxInf = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1)) m \[Omega] Abs[Z["\[ScriptCapitalI]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2)) m \[Omega] Abs[Z["\[ScriptCapitalI]"]]^2/(16*Pi)]; - FluxH = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1)) m \[Omega] Abs[Z["\[ScriptCapitalH]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2)) m \[Omega] Abs[Z["\[ScriptCapitalH]"]]^2/(16*Pi)]; + If[mode["Method"][[1]] == "Hyperboloidal", + FluxInf = I*m*(-I \[Omega] 4)((l+2)!/(l-2)!)1/(64\[Pi]*16)Abs[Z[[1]]]^2; + FluxH = I*m*(-I \[Omega] 4)((l+2)!/(l-2)!)1/(64\[Pi]*16)Abs[Z[[2]]]^2; + , + FluxInf = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1)) m \[Omega] Abs[Z["\[ScriptCapitalI]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2)) m \[Omega] Abs[Z["\[ScriptCapitalI]"]]^2/(16*Pi)]; + FluxH = If[EvenQ[l+m], (l-1)*(l+2)/(l*(l+1)) m \[Omega] Abs[Z["\[ScriptCapitalH]"]]^2/(4*Pi), (l*(l+1))/((l-1)*(l+2)) m \[Omega] Abs[Z["\[ScriptCapitalH]"]]^2/(16*Pi)]; + ]; <| "\[ScriptCapitalI]" -> FluxInf, "\[ScriptCapitalH]" -> FluxH |> ]; diff --git a/Kernel/ReggeWheelerRadial.m b/Kernel/ReggeWheelerRadial.m index 15cba33..bb24b28 100644 --- a/Kernel/ReggeWheelerRadial.m +++ b/Kernel/ReggeWheelerRadial.m @@ -59,7 +59,7 @@ Begin["`Private`"]; -(* ::Section::Closed:: *) +(* ::Section:: *) (*ReggeWheelerRadial*) @@ -269,7 +269,7 @@ ]; -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*ReggeWheelerRadial*) @@ -287,7 +287,7 @@ }; -(* ::Subsubsection:: *) +(* ::Subsubsection::Closed:: *) (*Static modes*) @@ -342,6 +342,11 @@ (* Potential *) pot = OptionValue["Potential"]; + (* Don't use ReggeWheelerRadial if Method = Hyperboloidal *) + If[OptionValue["Method"] == "Hyperboloidal", + Return[Null]; + ]; + (*Zerilli Potential only valid for Abs[s]=2*) If[MatchQ[pot, "Zerilli"]&&Abs[s]!=2, Message[ReggeWheelerRadial::szopt]; diff --git a/PacletInfo.wl b/PacletInfo.wl index 8001efe..1129d10 100644 --- a/PacletInfo.wl +++ b/PacletInfo.wl @@ -1,3 +1,5 @@ +(* ::Package:: *) + Paclet[ "Name" -> "ReggeWheeler", "Version" -> "0.3.0",