Nested NIntegrate with FindRoot

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}{2\pi\sigma_{x’}\sigma_{s’}}\exp\left({ -\frac{x’^2}{2\sigma_{x’}^2} }\right)\exp\left({ -\frac{s’^2}{2\sigma_{s’}^2} }\right).$ $

The tricky part(s) is that:

  1. $ z’$ in the $ K(N,z’,x,x’)$ needs to be solved for numerically using FindRoot and will have a $ s’$ dependence.
  2. The integration upper limit over $ ds’$ is a variable $ s$ .
  3. 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]