I have defined the following functions:
γ[r_, v_, rDet_] :=
Which[
Abs[r - v tDet] >= rDet, 0,
r + v tDet <= rDet, π,
True, ArcCos[((v tDet)^2 + r^2 - rDet^2)/(2 v tDet r)]];
ρ[r_, v_, v0_] := Exp[-(1/2) (r/rCloud)^2] r Exp[-(1/2) (v/v0)^2] v;
ΔΦ[v_] := (2 π )/λ c (4 (v*t)^2)/(rBeam)^2;
I want to weight ΔΦ
by ρ
and γ
. I do so by defining a function Awfc
that numerically integrates the product of the three functions over v
and r
and divides by the numerical integral over the two weighting functions:
Awfc[rDet_?NumericQ, v0_?NumericQ] :=
NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0]* Exp[I ΔΦ[v]], {r,0,∞}, {v,0,∞ }}]/
NIntegrate[ γ[r,v,rDet]*ρ[r,v,v0], {r,0,∞}, {v,0,∞}}];
I calculate Awfc
for some parameters rDet
and v0
AwfcTable =
ParallelTable[
{rDet, Awfc[rDet, 0.5 v0], Awfc[rDet, v0], Awfc[rDet,2 v0], Awfc[rDet,3.7 v0]},
{rDet,0.0005,0.010,0.0005}];
using these values for the other parameters:
tDet = 0.7;
t = 0.230;
rCloud = 0.0025;
λ = 780 10 ^-9;
c = λ/20;
v0 = 0.00588;
rDet = 0.008/2;
rBeam = 0.015;
All kernels throw warnings:
NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
I think that ΔΦ
is a hightly osciallating function so I don't see a way to fix this problem.
I use TableForm[Abs[AwfcTable]]
to get the amplitude of the complex numbers in AwfcTable
. The problem with the output is, that this amplitude is significantly larger than 1, especially for larger v0
and small rDet
. The function ΔΦ
I am calculating the weighted mean of, has a maximum amplitude of 1 so I think the weighted mean should also have a maximum of 1. (The complex phase is consistent with what I expect from a different approach to the problem I did using MATLAB)
Is there something wrong with my reasoning (meaning that the result Mathematica gives is actually correct), or is does the fault lie with my implementation of the problem?