Get the flux working for the 22-mode

ConvolveSource.m

@@ -6,7 +6,7 @@
 
 
 ConvolveSource[RF_, SO_] :=
-	If[SO["type"]=="PointParticleCircular",ConvolvePointParticleSourceCircular[RF,SO],Nothing];
+	If[SO["type"]=="PointParticleCircular", Return[ConvolvePointParticleSourceCircular[RF,SO]], Nothing];
 
 
 ConvolvePointParticleSourceCircular[RF_,SO_]:=
@@ -17,10 +17,10 @@
 		\[Omega]=m*Sqrt[1/r0^3];
 		rm2M=r0-2;
 		If[EvenQ[l+m],
-			PsiOddIn=RF["In"]["Psi"][r0];
-			dPsiOddIndr=RF["In"]["dPsidr"][r0];
-			PsiOddUp=RF["Up"]["Psi"][r0];
-			dPsiOddUpdr=RF["Up"]["dPsidr"][r0];
+			PsiOddIn=RF["In"][r0];
+			dPsiOddIndr=RF["In"]'[r0];
+			PsiOddUp=RF["Up"][r0];
+			dPsiOddUpdr=RF["Up"]'[r0];
 			n=(l-1)*(l+2);
 			np6M=n*r0+6;
 			denom=n*(n+2)+12*I*\[Omega];
@@ -32,17 +32,18 @@
 			dPsiIn=(-12*\[Omega]^2/(1-2/r0)*PsiOddIn+c*PsiOddIn+b*dPsiOddIndr)/conjdenom;
 			dPsiUp=(-12*\[Omega]^2/(1-2/r0)*PsiOddUp+c*PsiOddUp+b*dPsiOddUpdr)/denom;
 		,
-			PsiIn=RF["In"]["Psi"][r0];
-			dPsiIn=RF["In"]["dPsidr"][r0];
-			PsiUp=RF["Up"]["Psi"][r0];
-			dPsiUp=RF["Up"]["dPsidr"][r0];
+			PsiIn=RF["In"][r0];
+			dPsiIn=RF["In"]'[r0];
+			PsiUp=RF["Up"][r0];
+			dPsiUp=RF["Up"]'[r0];
 		];
+
 		Wronskian = PsiIn*dPsiUp - PsiUp*dPsiIn;
 		deltaPsi=SO["deltaPsi"];
 		deltadPsidr=SO["deltadPsidr"];
 		ZIn = (PsiUp*deltadPsidr - deltaPsi*dPsiUp)/Wronskian;
 		ZUp = (PsiIn*deltadPsidr - deltaPsi*dPsiIn)/Wronskian;
-		<|"ZInf"->ZUp,"ZHor"->ZIn|>
+		<|"ZInf"->ZUp, "ZHor"->ZIn|>
 	];
 
 

NumericalIntegration.m

@@ -27,10 +27,12 @@
 				psiBC=boundaryconditions["Psi"];
 				dpsidxBC=boundaryconditions["dPsidx"];
 				xBC=boundaryconditions["xBC"];
+
 				soln = Integrator[s,l,\[Omega],psiBC,dpsidxBC,xBC,xmax,ReggeWheelerPotential,$MachinePrecision];
+
 				newassoc=<|
-					"Psi"->Function[y,y1[y]/.soln],
-					"dPsidr"->Function[y,y2[y]/.soln],
+					"Psi"->soln[[1,1,2]],
+					"dPsidr"->soln[[1,2,2]],
 					"xmin"->xBC,
 					"xmax"->xmax					
 					|>;
@@ -77,8 +79,8 @@
 					xBC=boundaryconditions["xBC"];
 					soln = Integrator[s,l,\[Omega],psiBC,dpsidxBC,xBC,xmin,ReggeWheelerPotential,$MachinePrecision];
 					newassoc= <|
-						"Psi"->Function[y,y1[y]/.soln],
-						"dPsidr"->Function[y,y2[y]/.soln],
+						"Psi"->soln[[1,1,2]],
+						"dPsidr"->soln[[1,2,2]],
 						"xmin"->xmin,
 						"xmax"->xBC
 					|>
@@ -106,7 +108,7 @@
 
 
 (*should this be in a module for y1 and y2?*)
-Integrator[s_,l_,\[Omega]_,y1BC_,y2BC_,xBC_,xend_,potential_,precision_]:=
+Integrator[s_,l_,\[Omega]_,y1BC_,y2BC_,xBC_,xend_,potential_,precision_]:=Module[{y1,y2,x},
 	NDSolve[
 		{y1'[x]==y2[x],(1-2/x)^2*y2'[x]+2(1-2/x)/x^2*y2[x]+(\[Omega]^2-potential[s,l,x])*y1[x]==0,y1[xBC]==y1BC,y2[xBC]==y2BC},
 		{y1,y2},
@@ -115,15 +117,16 @@
 		MaxSteps->Infinity,
 		WorkingPrecision->precision,
 		InterpolationOrder->All
-	];
+		]
+	]
 
 
 (*boundary conditions for odd-parity Regge-Wheeler eqn.*)
 (*currently only working for s=2*)
 (*SetAttributes[ReggeWheelerInBC, {NumericFunction}];
 SetAttributes[ReggeWheelerUpBC, {NumericFunction}];*)
 
-ReggeWheelerInBC[s_Integer,l_Integer,\[Omega]_,workingprecision_Integer]:=
+ReggeWheelerInBC[s_Integer,l_Integer,\[Omega]_,workingprecision_]:=
 	Module[{rm2M,p,ptrys,expeh,Dexpeh,done=False,delReh,Reh,rstar,drstardr,
 	nmax,Xn,n,bk,denominator,f1=0,f2=0,f3=0,last,Bkm1=1,Bkm2=0,Bkm3=0,next,psi,dpsidr,om,precision=workingprecision+10,count},
 		rm2M=4*^-2*\[Omega]^2/(l*(l+1));
@@ -186,7 +189,7 @@
 	<|"Psi"->N[psi,precision],"dPsidx"->N[dpsidr,precision],"xBC"->N[Reh,precision]|>
 ]
 
-ReggeWheelerUpBC[s_Integer,l_Integer,\[Omega]_,xmax_,workingprecision_Integer]:=
+ReggeWheelerUpBC[s_Integer,l_Integer,\[Omega]_,xmax_,workingprecision_]:=
 	Module[{An=1,Anm1=0,Anm2=0,Anm3=0,Nmax=75,NNmax=1000,n,nn,rstart,rstar,drstardr,increment=0,
 	lastincrement=1*^40,S=0,lastS=0,dS=0,lastdS=0,count=0,r,rn,np,continue=True,precision=workingprecision+10,om,BCinc},
 		(*rstarstart=xmax+2*Log[xmax/2-1]+5*Pi/Abs[om];*)
@@ -242,10 +245,10 @@
 
 
 (*radial potentials for ODE*)
-SetAttributes[ReggeWheelerPotential,{NumericFunction}];
+(*SetAttributes[ReggeWheelerPotential,{NumericFunction}];*)
 SetAttributes[ZerilliPotential,{NumericFunction}];
 
-ReggeWheelerPotential[s_Integer,l_Integer,x_Numeric]:=(1-2/x)*(2(1-s^2)+l*(l+1)*x)/x^3;
+ReggeWheelerPotential[s_Integer,l_Integer,x_]:=(1-2/x)*(2(1-s^2)+l*(l+1)*x)/x^3;
 
 (*only for spin s=2*)
 ZerilliPotential[l_Integer,x_Numeric]:=

ReggeWheelerMode.m

@@ -4,30 +4,36 @@
   {"ReggeWheeler`ReggeWheelerSource`",
    "ReggeWheeler`ReggeWheelerRadial`",
    "ReggeWheeler`ConvolveSource`",
-   "KerrGeodesics`KerrGeoOrbit`"}
+   "KerrGeodesics`KerrGeoOrbit`",
+   "KerrGeodesics`OrbitalFrequencies`"}
 ];
 
 ReggeWheelerModeObject::usage = "ReggeWheelerModeObject[assoc] an object which contains a Regge Wheeler mode.";
 
-ReggeWheelerPointParticleMode::usage = "ReggeWheelerPointParticleMode[s, l, m, n, k, orbit] solves the Regge Wheeler equation with a point particle source.";
+ReggeWheelerPointParticleMode::usage = "ReggeWheelerPointParticleMode[s, l, m, n, orbit] solves the Regge Wheeler equation with a point particle source.";
 
 Begin["`Private`"];
 
 
 (*defined currently for circular orbits*)
-ReggeWheelerPointParticleMode[s_Integer, l_Integer, m_Integer, n_Integer, k_Integer, orbit_KerrGeoOrbitFunction] :=
+ReggeWheelerPointParticleMode[s_Integer, l_Integer, m_Integer, n_Integer, orbit_KerrGeoOrbitFunction] :=
 	Module[{source},
+
+		If[!PossibleZeroQ[orbit["a"]], Print["Regge-Wheeler perturbations are only for Schwarzschild (a=0) black holes"]; Return[];];
+
 		source = ReggeWheelerPointParticleSource[s, l, m, orbit];
-		Return[ReggeWheelerPointParticleMode[s, l, m, n, k, orbit, source]];
+		Return[ReggeWheelerPointParticleMode[s, l, m, n, orbit, source]];
 	];
 
 
-ReggeWheelerPointParticleMode[s_Integer, l_Integer, m_Integer, n_Integer, k_Integer, orbit_KerrGeoOrbitFunction, source_ReggeWheelerSourceObject] :=
+ReggeWheelerPointParticleMode[s_Integer, l_Integer, m_Integer, n_Integer, orbit_KerrGeoOrbitFunction, source_ReggeWheelerSourceObject] :=
  Module[{assoc, R, S, \[Omega], \[CapitalOmega]r, \[CapitalOmega]\[Phi], \[CapitalOmega]\[Theta], Z, Fluxes},
-  {\[CapitalOmega]r, \[CapitalOmega]\[Theta], \[CapitalOmega]\[Phi]} = orbit["Frequencies"];
+  (*{\[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"]]];
   \[Omega] = m \[CapitalOmega]\[Phi] + n \[CapitalOmega]r;
 
   R = ReggeWheelerRadial[s, l, \[Omega]];
+
 
   Z = ReggeWheeler`ConvolveSource`Private`ConvolveSource[R, source];
 
@@ -43,10 +49,10 @@
 	|>
 	];
 
-  assoc = <| "l" -> l,
+  assoc = <| "s" -> s, 
+             "l" -> l,
              "m" -> m,
              "n" -> n,
-             "k" -> k,
              "Type" -> "PointParticleCircular",
              "Homogeneous" -> R,
              "Particular" -> Z,
@@ -57,7 +63,7 @@
 ]
 
 
-Format[ReggeWheelerModeObject[assoc_]] := "ReggeWheelerModeObject[<<>>]";
+Format[ReggeWheelerModeObject[assoc_]] := "ReggeWheelerModeObject["<>ToString[assoc["s"]]<>","<>ToString[assoc["l"]]<>","<>ToString[assoc["m"]]<>","<>ToString[assoc["n"]]<>","<>"<<>>]";
 
 ReggeWheelerModeObject[assoc_][string_] := assoc[string]
 

ReggeWheelerRadial.m

@@ -19,10 +19,10 @@
       method = {"MST"};
       solFuncs = $Failed,
     {"NumericalIntegration", "rmin" -> _, "rmax" -> _},
-      method = Association[OptionValue[Method]];
-      solFuncs = OptionValue["BoundaryConditions"] /.
-        {"In" -> ReggeWheeler`NumericalIntegration`Private`PsiIn[s, l, \[Omega], method["rmin"], method["rmax"]],
-         "Up" -> ReggeWheeler`NumericalIntegration`Private`PsiUp[s, l, \[Omega], method["rmin"], method["rmax"]]};
+      method = Association[OptionValue[Method][[2;;]]];
+      solFuncs = 
+        {ReggeWheeler`NumericalIntegration`Private`PsiIn[s, l, \[Omega], method["rmin"], method["rmax"]],
+         ReggeWheeler`NumericalIntegration`Private`PsiUp[s, l, \[Omega], method["rmin"], method["rmax"]]};
   ];
 
   assoc = Association[
@@ -54,13 +54,13 @@
   ];  
 ];
 
-(*Derivative[n_][ReggeWheelerRadialFunction[s_, l_, \[Omega]_, assoc_]][r_?NumericQ] := Module[{},
+Derivative[n_][ReggeWheelerRadialFunction[s_, l_, \[Omega]_, assoc_]][r_?NumericQ] := Module[{},
   If[
     Head[assoc["BoundaryConditions"]] === List,
-    Return[Association[MapThread[#1 -> Derivative[n][#2][r] &, {assoc["BoundaryConditions"], assoc["SolutionFunctions"]}]]], 
-    Return[Derivative[n][assoc["SolutionFunctions"]][r]]
+    Return[Association[MapThread[#1 -> #2["dPsidr"][r] &, {assoc["BoundaryConditions"], assoc["SolutionFunctions"]}]]], 
+    Return[assoc["SolutionFunctions"]["dPsidr"][r]]
   ];  
-];*)
+];
 
 
 End[]

ReggeWheelerSource.m

@@ -10,7 +10,7 @@
 
 
 ReggeWheelerPointParticleSource[s_,l_,m_, orbit_] :=
- If[s==-2 && orbit["e"] == 0 && Abs[orbit["Inclination"]] == 1,
+ If[orbit["e"] == 0 && Abs[orbit["Inclination"]] == 1,
          Return[ReggeWheelerPointParticleSourceCircular[s,l,m,orbit]],
          Print["No point-particle source yet available for those parameters"];]