I am trying to numerically integrate a function using nested NIntegrate:

$$F(N,x,s)=int_{-infty}^s int_{-infty}^{+infty} K(N,z’,x,x’) g_{x’,s’} dx’ds’ $$

where the kernel of the integration, $K(N,z’x,x’)$, is a messy expression defined in the mathematica code below, and $g_{x’,s’}$ is a bi-variate gaussian defined by:

$$ g_{x,s’}=frac{n}{2pisigma_{x’}sigma_{s’}}expleft({ -frac{x’^2}{2sigma_{x’}^2} }right)expleft({ -frac{s’^2}{2sigma_{s’}^2} }right).$$

The tricky part(s) is that:

- $z’$ in the $K(N,z’,x,x’)$ needs to be solved for numerically using FindRoot and will have a $s’$ dependence.
- The integration upper limit over $ds’$ is a variable $s$.
- I suspect the kernel is oscillatory with $N$ (denoted “Kernel” in the code below) so maybe an averaging of the kernel over $N$ can be done to simplify the kernel and eliminate $N$ if the integrations prove to be too time consuming.

At the end, I would like a function, F(N,x,s), that would be able to plot across $s$ for a given $(N,x)$ values i.e. Plot(F(a,b,s,{s,-1e-5,1e-5}).

```
(*Constants*)
e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;
re = 2.81794092*10^-15;
c = 2.99792458*10^8;
n = -10^-10/e;
KK = 1;
lw = 0.026;
kw = (2 Pi)/lw;
gamma = 4000/0.511;
beta = Sqrt(1 - 1/gamma^2);
sigmaS = 10^-5;
sigmaX = 30*10^-6;
coeff = n/(2 Pi*sigmaS*sigmaX) Exp(-(xprime^2/(2 sigmaX^2)))*
Exp(-(sprime^2/(2 sigmaS^2)));
(*Preliminary Equations*)
rs2 = {zprime, xprime + KK/(gamma*kw) Sin(kw*zprime), 0};
ro2 = {(NN + 10000)*lw, x + KK/(gamma*kw) Sin(kw*(NN + 10000)*lw), 0};
betas = {beta - KK^2/(2 gamma^2) Cos(kw*zprime)^2,KK/gamma Sin(kw*zprime), 0};
betao = {beta - KK^2/(2 gamma^2) Cos(kw*(NN + 10000)*lw)^2,KK/gamma Sin(kw*(NN + 10000)*lw), 0};
betaDot = {(c*KK^2*kw)/(2 gamma^2)Sin(2 kw*zprime), -((KK*c*kw)/gamma) Sin(kw*zprime), 0};
deltar2 = ro2 - rs2;
Rgam2 = Sqrt(deltar2((1))^2 + deltar2((2))^2);
Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross(deltar2/Rgam2,Cross(deltar2/Rgam2 - betas, betaDot))/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);
sumElong = (Ec2((1)) + Erad2((1)));
sumEtran = (Ec2((2)) + Erad2((2)));
(*Numerical Functions*)
ZPRIME(NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ) := zprime /.FindRoot(s - sprime == (Sqrt(gamma^2 + KK^2) (EllipticE(kw*(NN + 10000)*lw,KK^2/(gamma^2 + KK^2)) - EllipticE(kw zprime, KK^2/(gamma^2 + KK^2))))/(gamma kw) -beta Sqrt(((NN + 10000)*lw - zprime)^2 + (x - xprime + (KK Sin(kw *(NN + 10000)*lw))/(gamma kw) - (KK Sin(kw zprime))/(gamma kw))^2), {zprime, 0})
Kernel = coeff re/gamma (sumElong*betao((1)) + sumEtran*betao((2)))/.{zprime -> ZPRIME(NN, x, xprime, s, sprime)};
FNxprimesprime(NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ):= Kernel
FNsprime(NN_?NumericQ, x_?NumericQ, s_?NumericQ, sprime_?NumericQ) :=NIntegrate(FNxprimesprime(NN, x, xprime, s, sprime), {xprime, -300/10^6, 300/10^6})
FN(NN_?NumericQ,x_?NumericQ, s_?NumericQ) := NIntegrate(FNsprime(NN,x, s, sprime), {sprime,-10^-4, s})
lst1 = Table({ss, FN(0,0, ss), PrecisionGoal -> 5) // Quiet}, {ss, -10^-5, 10^-5, 10^-6})
ListPlot(lst1)
```