(* 

This notebook contains the modal solution to Lilley's equation for a point source (1). Also included here are the symmetric ring source (2)and weighted ring source (3).

This work is used in the paper "Stone et al., Aeroacoustic catastrophes: upstream beaming of jet noise."

In all cases the single simple turning-point solution is used. Extensions can be made using the methods described in the above paper. 

  *)

(* ------------------------- 1. Point source calculation -------------------------- *)

(* The modal solution is calculated quickly by noting the separable nature of the sum. Compute the azimuthal components first using. *)

ExpFunc[Δϕ_,N_]:=Table[Exp[I*n*Δϕ],{n,-N,N}]

(* Compute the radial solutions *)

myArgm[x_]:=Mod[Arg[x],2 π,-(π/2)]

q[θ_?NumericQ,r_?NumericQ,MJ_?NumericQ]:=(Sqrt[(((1-MJ ((Sech[2 r])^2) Cos[θ]))^2)-((Cos[θ])^2)])

Q[ν_?NumericQ,θ_?NumericQ,r_?NumericQ,MJ_?NumericQ]:=(((Sqrt[Abs[(r^2) ((q[θ,r,MJ])^2)-(ν^2)]]) Exp[(1/2) I (myArgm[r q[θ,r,MJ]-ν]+myArgm[r q[θ,r,MJ]+ν])])/r)

rδ[ν_?NumericQ,θ_?NumericQ,M_?NumericQ]:=p/. FindRoot[Abs[(p^2) ((q[θ,p,M])^2)-(ν^2)] Exp[I (myArgm[p q[θ,p,M]-ν]+myArgm[p q[θ,p,M]+ν])],{p,1}]

rδ3[ν_,θ_,M_]:=Abs[rδ[ν,θ,M]]

DiffP[θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ]:=(∂_x)((x^2) ((((1-MJ ((Sech[2 x])^2) Cos[θ]))^2)-((Cos[θ])^2)))/. x->rs

P[θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ]:=(1/(rs DiffP[θ,rs,MJ]))

ZetaN[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ]:=NIntegrate[Abs[Q[n,θ,r,MJ]],{r,rδ3[n,θ,MJ],rs}]

EtaN[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=With[{sign=Sign[rs-rδ[n,θ,MJ]]},-sign ((((3/2) sign k ZetaN[n,θ,rs,MJ]))^(2/3))]

AmplitudeFunc[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=With[{sign=Sign[rs-rδ[n,θ,MJ]]},If[sign==0,(Sqrt[(2/k)]) (((k P[θ,rs,MJ]))^(1/3)),(Sqrt[(2/k)]) (Sqrt[((Sqrt[-sign EtaN[n,θ,rs,MJ,k]])/(Abs[rs Q[n,θ,rs,MJ]]))])]]

AiMode[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=AiryAi[EtaN[n,θ,rs,MJ,k]]

ReverseFunc[array1_, array0_]:=(* array1 is a sequence of terms a[1] to a[N]*)Flatten[{Reverse[array1],array0,array1}]

NTermQuick[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=AmplitudeFunc[(n/k),θ,rs,MJ,k] AiMode[(n/k),θ,rs,MJ,k] Exp[I k (ZetaN[(n/k),θ,100 Sin[θ],(9/10)]-100 ((Sin[θ])^2))]

NTermQuickArray[N_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=Table[NTermQuick[n,θ,rs,MJ,k],{n,1,N}]

NTermsForAPoint[N_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=ReverseFunc[NTermQuickArray[N,θ,rs,MJ,k], NTermQuick[0,θ,rs,MJ,k]]

(* ------------ Example computations of point source ------------  *)

(* N.B. Use radians, not degrees. *)

k0=28.27; rs = 0.75; Ns = 30; MJ = 0.9; Δϕ = 0;cInf = 340;

radialSol = Table[( (I/(cInf*k0*(((1-MJ*((Sech[2*rs])^2)*Cos[i*π/180]))^2))))*NTermsForAPoint[Ns,i*π/180,rs,MJ,k0],{i, 2,100, 1}];

aziAngles = ExpFunc[Δϕ,30];

totSolΔϕ = radialSol.aziAngles ;

ListLinePlot[Transpose[{Table[i*π/180,{i,1,100,2}],Abs[totSolΔϕ]}], AxesLabel->{"θ","Abs[(Subscript[G, ω])/(Subscript[g, ω])]"}]

(* ------------------------- 2. Ring source calculation -------------------------- *)

(* The ring source is calculated using the point source solution of the previous section. *)

AmplitudeFuncRS[n_?NumericQ,θ_?NumericQ,rs_?NumericQ,MJ_?NumericQ,k_?NumericQ]:=With[{sign=Sign[rs-rδ[n,θ,MJ]]},((Sqrt[-sign EtaN[n,θ,rs,MJ,k]])/(Abs[rs Q[n,θ,rs,MJ]]))]

Dω[θ_, rs_, MJ_, k0_]:=((4*π)/((cInf^2)*(k0^3)*(((1 - MJ*((Sech[2*rs])^2)*Cos[θ]))^4)))*Sum[AmplitudeFuncRS[n/k0,θ,rs,MJ,k0]*((AiMode[n/k0, θ,rs,MJ,k0])^2),{n,-30, 30}]

ringSource = Dω[#, rs, MJ, k0]&/@Table[i*π/180,{i,1, 158,0.5}];

(* ------------------- 3. Weighted ring source calculation ---------------------- *)

(* We need the radial solutions from the point source solution. I.e., needs 'radialSol' from section 1. This code generates one cross-section at ϕ=const. of the weighted ring source calculation.  *)

(* Example parameters *)

αT = 1; ϕ = π/2+28*π/180;

YNewFunc[ αT_,ϕ_ ]:=(* ϕP is at π/2 *)Table[((Sqrt[π])/(2 (Sqrt[αT])))I (E^(-(((m-n) (m-n-4 I *αT*(ϕ - π/2)))/(4 *αT)))) (-Erfi[((-m+n+2 I *αT* (π-π/2))/(2 (Sqrt[αT])))]+Erfi[((-m+n-2 I *αT*(π+π/2))/(2 (Sqrt[αT])))]),{n,-40,40},{m,-40,40}]//N;

siftFunc[a_,n_]:=a[[1]][[n]]

D1 = {Outer[Times, radialSol[[#]],Conjugate[radialSol[[#]]]]&/@Table[i,{i, 1, Length[radialSol]}]};

D1 = Flatten[D1,1];

Yarrays = YNewFunc[αT,ϕ];

dϕExample ={D1[[#]]*Yarrays&/@Table[i,{i,1,Length[radialSol]}]};;

dϕExample2=Flatten[{Total[Flatten[siftFunc[dϕExample,#]]]}&/@Table[i,{i,1,Length[radialSol]}]];

ListLinePlot[Transpose[{Table[i*π/180,{i,1,100,2}],Abs[dϕExample2]}] , AxesLabel->{"θ","Abs[(Subscript[D, ω])]"}]