62 lines
2.0 KiB
Mathematica
62 lines
2.0 KiB
Mathematica
(*
|
|
Usage: Please set the fowllowing variables from outside as
|
|
command line arguments, like this:
|
|
wolframscript -file calculate_modes.wl psi4_file r rfilename ifilename
|
|
for example:
|
|
wolframscript -file calculate_modes.wl psi4.txt 2.5 rmodes.txt imodes.txt
|
|
*)
|
|
|
|
(* Parsing arguments *)
|
|
argv = Rest@$ScriptCommandLine;
|
|
(* Print[argv] *)
|
|
Check[
|
|
psi4file = argv[[1]];
|
|
(* {l, m, r} = Interpreter["Number"]/@argv[[2;;4]]; *)
|
|
r = Interpreter["Number"]@argv[[2]];
|
|
rfilename = argv[[3]];
|
|
ifilename = argv[[4]];
|
|
,
|
|
Print["Error: Check your command line arguments!"];Abort[]
|
|
];
|
|
|
|
(* Set Assumptions *)
|
|
$Assumptions = Pi > theta > 0 && 2 Pi >= phi >= 0 &&\
|
|
Element[l, Reals] && Element[m, Reals] && -l<=m<=l;
|
|
|
|
Y[s_, l_, m_, theta_, phi_] :=
|
|
(-1)^m*
|
|
(Sqrt[(((l + m)!*(l - m)!)*(2*l + 1))/((((l + s)!*(l - s)!)*4)*Pi)]*
|
|
Sin[theta/2]^(2*l))*
|
|
Sum[
|
|
(((Binomial[l - s, r]*Binomial[l + s, r + s - m])*(-1)^(l - r - s))*
|
|
E^((I*m)*phi))*Cot[theta/2]^(2*r + s - m),
|
|
{r, 0, l - s}
|
|
]
|
|
|
|
(* Y[l_, m_, theta_, phi_] := ((-1)^(l + m) E^(I m phi)
|
|
Binomial[-2 + l, -2 - m] Cot[theta/2]^(-2 -m)
|
|
Sqrt[((1 + 2 l) (l - m)! (l + m)!)/((-2 + l)! (2 + l)!)]
|
|
Hypergeometric2F1[-2 - l, -l - m, -1 - m,
|
|
-Cot[theta/2]^2] Sin[theta/2]^(2 l))/(2 Sqrt[Pi]) *)
|
|
|
|
psi4[theta_, phi_] = ToExpression[Import[psi4file, "Text"] ]/.{t->0};
|
|
|
|
mode[l_, m_] := Module[{f},
|
|
(* f = Simplify[Conjugate[Y[-2, l,m,theta,phi] ]*psi4[r,theta,phi]*Sin[theta] ]; *)
|
|
f = Conjugate[Y[-2, l,m,theta,phi] ]*psi4[theta,phi]*Sin[theta];
|
|
Quiet[NIntegrate[f,{theta, 0, Pi}, {phi, 0, 2 Pi}] ]
|
|
]
|
|
|
|
(* mode[0,0,2.2]//Print *)
|
|
|
|
(* result = ParallelTable[mode[l,m],{l,1,8},{m,-l,l}]//Flatten; *)
|
|
lmlist = Flatten[Table[{l, m}, {l, 1, 8}, {m, -l, l}], 1]
|
|
result = ParallelMap[mode@@#&, lmlist]
|
|
rmodes = Re/@result;
|
|
imodes = Im/@result;
|
|
Export[rfilename, rmodes, "Table"]
|
|
Export[ifilename, imodes, "Table"]
|
|
|
|
(* result = mode[l,m]
|
|
Print[Re@result, " ", Im@result] *)
|