Re: Efficient BesselJ and LegendreP Evaluation
- To: mathgroup at smc.vnet.net
- Subject: [mg74717] Re: [mg74706] Efficient BesselJ and LegendreP Evaluation
- From: "Antonio Neves" <aneves at gmail.com>
- Date: Mon, 2 Apr 2007 06:53:24 -0400 (EDT)
- References: <200704010818.EAA07376@smc.vnet.net>
Sorry for the greek letters in the code below I have replaced them
with their names.
Clear["Global`*"];
Off[General::spell];
Off[General::spell1];
lambda = 0.8;
w = 2500;
NA = 1.25;
f = 1700;
d = 50.;
nGlass = 1.4883 + 1.9023/(10^2*lambda^2) - 4.25016/(10^3*lambda^4);
nWater = 1.3253 + 2.7553/(10^3*lambda^2) + 3.779/(10^5*lambda^4);
nAir = 1.0001;
alphaMax = ArcSin[NA/nWater];
ko = (2*Pi*nAir)/lambda;
Nmax = 20;
GridNM = Table[{n, m}, {n, 1, Nmax}, {m, 0, n}];
alphaCrit = ArcSin[nWater/nGlass];
If[alphaCrit < alphaMax, alphaLim = alphaCrit, alphaLim = alphaMax];
cosA2[(alpha_)?NumericQ] := Sqrt[1 - (nGlass*(Sin[alpha]/nWater))^2];
ts[(alpha_)?NumericQ] := 2*nGlass*(Cos[alpha]/(nGlass*Cos[alpha] +
nWater*cosA2[alpha]));
tp[(alpha_)?NumericQ] := 2*nGlass*(Cos[alpha]/(nWater*Cos[alpha] +
nGlass*cosA2[alpha]));
DBessel[m_, x_] := (1/2)*(BesselJ[-1 + m, x] - BesselJ[1 + m, x]);
DLegendre[n_, m_, x_] := ((-1 - n)*x*LegendreP[n, m, x] + (1 - m +
n)*LegendreP[1 + n, m, x])/(-1 + x^2);
ITM[(n_Integer)?Positive, (m_Integer)?NonNegative, zo_Real, rhoo_Real] :=
Block[{func}, func[(alpha_)?NumericQ] :=
Sqrt[Cos[alpha]]*Exp[-(f*(Sin[alpha]/w))^2]*Exp[(-I)*ko*nWater*zo*cosA2[alpha]]*
Exp[(-I)*ko*d*(nGlass*Cos[alpha] -
nWater*cosA2[alpha])]*(m^2*ts[alpha]*LegendreP[n, m,
cosA2[alpha]]*(BesselJ[m,
ko*nGlass*rhoo*Sin[alpha]]/(ko*nGlass*rhoo*Sin[alpha])) -
tp[alpha]*(nGlass/nWater)^2*Sin[alpha]^2*DBessel[m,
ko*nGlass*rhoo*Sin[alpha]]*DLegendre[n, m, cosA2[alpha]] +
I*m*(ts[alpha]*DBessel[m,
ko*nGlass*rhoo*Sin[alpha]]*LegendreP[n, m, cosA2[alpha]] -
tp[alpha]*(nGlass/nWater)^2*Sin[alpha]^2*DLegendre[n, m,
cosA2[alpha]]*
(BesselJ[m,
ko*nGlass*rhoo*Sin[alpha]]/(ko*nGlass*rhoo*Sin[alpha]))));
NIntegrate[func[alpha], {alpha, 0, alphaLim}, Compiled -> True, Method
-> DoubleExponential, MaxRecursion -> 1000, SingularityDepth ->
1000]];
zo = Random[Real, {-10., 10.}];
rhoo = Random[Real, {0., 10.}];
Timing[MI1 = Apply[ITM[#1, #2, zo, rhoo] & , GridNM, {2}];]
On 4/1/07, Antonio <aneves at gmail.com> wrote:
> Dear group members,
>
> I have a function that I want to integrate for various variable
> parameter in the most precise and fastest way using mathematica
> instead of external compile C program. Below is a code to evaluate a
> numerical integral, which has BesselJ, LegendreP and their derivative.
> These needs to be evaluated for a set of {n,m}. Below in this example
> I have {n,1,20} and {m,0,n}.
>
> There is also a singularity issue to be solved:"NIntegrate::singd :
> NIntegrate's singularity handling has failed at point
> {=CE=B1}={6.007268658180993`*^-9} for the specified precision goal. Try
> using larger values for any of $MaxExtraPrecision or the options
> WorkingPrecision, or SingularityDepth and MaxRecursion."
>
> Any comments are welcomed.
>
> Antonio
>
>
>
> Clear["Global`*"];
> Off[General::spell];
> Off[General::spell1];
> =CE=BB = 0.8;
> w = 2500;
> NA = 1.25;
> f = 1700;
> d = 50.;
> nGlass = 1.4883 + 1.9023/(10^2*=CE=BB^2) - 4.25016/(10^3*=CE=BB^4);
> nWater = 1.3253 + 2.7553/(10^3*=CE=BB^2) + 3.779/(10^5*=CE=BB^4);
> nAir = 1.0001;
> =CE=B1Max = ArcSin[NA/nWater];
> ko = (2*Pi*nAir)/=CE=BB;
> Nmax = 20;
> GridNM = Table[{n, m}, {n, 1, Nmax}, {m, 0, n}];
> =CE=B1Crit = ArcSin[nWater/nGlass];
> If[=CE=B1Crit < =CE=B1Max, =CE=B1Lim = =CE=B1Crit, =CE=B1Lim = =CE=B1Ma=
> x];
> cosA2[(=CE=B1_)?NumericQ] := Sqrt[1 - (nGlass*(Sin[=CE=B1]/nWater))^2];
> ts[(=CE=B1_)?NumericQ] := 2*nGlass*(Cos[=CE=B1]/(nGlass*Cos[=CE=B1] +
> nWater*cosA2[=CE=B1]));
> tp[(=CE=B1_)?NumericQ] := 2*nGlass*(Cos[=CE=B1]/(nWater*Cos[=CE=B1] +
> nGlass*cosA2[=CE=B1]));
> DBessel[m_, x_] := (1/2)*(BesselJ[-1 + m, x] - BesselJ[1 + m, x]);
> DLegendre[n_, m_, x_] := ((-1 - n)*x*LegendreP[n, m, x] + (1 - m +
> n)*LegendreP[1 + n, m, x])/(-1 + x^2);
> ITM[(n_Integer)?Positive, (m_Integer)?NonNegative, zo_Real,
> =CF=81o_Real] :=
> Block[{func}, func[(=CE=B1_)?NumericQ] := Sqrt[Cos[=CE=B1]]*Exp[-(f*(S=
> in[=CE=B1]/
> w))^2]*Exp[(-I)*ko*nWater*zo*cosA2[=CE=B1]]*
> Exp[(-I)*ko*d*(nGlass*Cos[=CE=B1] -
> nWater*cosA2[=CE=B1])]*(m^2*ts[=CE=B1]*LegendreP[n, m, cosA2[=CE=B1]]*(Bess=
> elJ[m,
> ko*nGlass*=CF=81o*Sin[=CE=B1]]/(ko*nGlass*=CF=81o*Sin[=CE=B1])) -
> tp[=CE=B1]*(nGlass/nWater)^2*Sin[=CE=B1]^2*DBessel[m,
> ko*nGlass*=CF=81o*Sin[=CE=B1]]*DLegendre[n, m, cosA2[=CE=B1]] +
> I*m*(ts[=CE=B1]*DBessel[m, ko*nGlass*=CF=81o*Sin[=CE=B1]]*LegendreP=
> [n, m,
> cosA2[=CE=B1]] - tp[=CE=B1]*(nGlass/nWater)^2*Sin[=CE=B1]^2*DLegendre[n, m,
> cosA2[=CE=B1]]*
> (BesselJ[m, ko*nGlass*=CF=81o*Sin[=CE=B1]]/(ko*nGlass*=CF=81o*Si=
> n[=CE=B1]))));
> NIntegrate[func[=CE=B1], {=CE=B1, 0, =CE=B1Lim}, Compiled -> True, Method ->
> DoubleExponential,
> MaxRecursion -> 1000, SingularityDepth -> 1000]];
> zo = Random[Real, {-10., 10.}];
> =CF=81o = Random[Real, {0., 10.}];
> Timing[MI1 = Apply[ITM[#1, #2, zo, =CF=81o] & , GridNM, {2}];]
>
>
>
- References:
- Efficient BesselJ and LegendreP Evaluation
- From: "Antonio" <aneves@gmail.com>
- Efficient BesselJ and LegendreP Evaluation