# Plot3D including a double integral

I was trying to run the following code that should generate a 3D plot. The code seems to be correct, but for some reason it doesn't give a plot! I have already the plot with MATLAB (numbers on z axis might be incorrect) and wanted to redraw it with Mathematica, but I couldn't. First, I was trying to calculate the integral but it includes two integration variables $$t_2$$ and $$t_3$$, so numerical integration is not possible, unless I put the integral inside Plot3D. I have attached the MATLAB plot below.

Here is the code:

c = 3 10^8;
L = 0.07;
\[Epsilon]0 = 8.85 10^-12;
hbar = 10^-34;
NN = 1.2*10^17;

\[Mu]13 = 2.06 10^-29;
\[Mu]23 = 2.06 10^-29;
\[Mu]24 = 2 10^-29;
\[Mu]14 = 2 10^-29;

\[CapitalGamma]41 = 18840000;
\[CapitalGamma]31 = 0.5 \[CapitalGamma]41;
\[CapitalGamma]11 = 0.02 \[CapitalGamma]41;
\[CapitalGamma]21 = 0.04 \[CapitalGamma]41;
\[CapitalGamma]22 = 0.02 \[CapitalGamma]41;
\[CapitalGamma]24 = \[CapitalGamma]41;
\[CapitalGamma]42 = \[CapitalGamma]24;

\[CapitalDelta]1 = 130 \[CapitalGamma]41;
\[CapitalDelta]3 = -3 \[CapitalGamma]41;
\[CapitalDelta]2 = 3 \[CapitalGamma]41;

\[CapitalOmega]2 = 15 \[CapitalGamma]41;
\[CapitalOmega]3 = 15 \[CapitalGamma]41;

\[Lambda]13 = 795 10^-9;
\[Lambda]23 = 795 10^-9;
\[Lambda]14 = 780 10^-9;
\[Lambda]24 = 780 10^-9;

\[Omega]24 = 2 \[Pi]/\[Lambda]24;
\[Omega]14 = 2 \[Pi]/\[Lambda]14;
\[Omega]13 = 2 \[Pi]/\[Lambda]13;
\[Omega]23 = 2 \[Pi]/\[Lambda]23;

\[Omega]1 = \[Omega]13 - \[CapitalDelta]1;
\[Omega]2 = \[Omega]24 - \[CapitalDelta]2;
\[Omega]3 = \[Omega]24 - \[CapitalDelta]3;
\[Omega]s1 = \[Omega]23 - \[CapitalDelta]1;
\[Omega]s2 = \[Omega]24 - \[CapitalDelta]2;
\[Omega]s3 = \[Omega]14 - \[CapitalDelta]3;

k1 = \[Omega]1/c;
k2 = \[Omega]2/c;
k3 = \[Omega]3/c;

Den1 = \[Epsilon]0 hbar^5 ((\[CapitalGamma]31 +
I \[CapitalDelta]1)) ((\[CapitalGamma]21 + I  \[Delta]2 +
I \[Delta]3) (\[CapitalGamma]41 + I  \[Delta]2 +
I  \[Delta]3 + I \[CapitalDelta]2) + \[CapitalOmega]2*
Conjugate[\[CapitalOmega]2]) ((\[CapitalGamma]11 +
I  \[Delta]3) (\[CapitalGamma]41 + I  \[Delta]3 +
I \[CapitalDelta]3) + \[CapitalOmega]3*
Conjugate[\[CapitalOmega]3]);

Den2 = \[Epsilon]0 hbar (4 (\[Delta]2 - \[CapitalDelta]2 +
I \[CapitalGamma]42) (\[Delta]2 +
I \[CapitalGamma]22) + \[CapitalOmega]2*
Conjugate[\[CapitalOmega]2]);

Den3 = \[Epsilon]0 hbar (4 (\[Delta]3 - \[CapitalDelta]3 +
I \[CapitalGamma]41) (\[Delta]3 +
I \[CapitalGamma]11) + \[CapitalOmega]3*
Conjugate[\[CapitalOmega]3]);

\[Kappa] = (2 NN \[Mu]13 \[Mu]24 \[Mu]23 (\[Mu]14^3) )/Den1;

\[Chi]S2 = (-4 I NN  (\[Mu]24^2) (\[Delta]2 + I \[CapitalGamma]22) )/
Den2;

\[Chi]S3 = (-4 I NN  (\[Mu]14^2) (\[Delta]3 + I \[CapitalGamma]11) )/
Den3;

ks1 = (\[Omega]s1 - (\[Delta]2 + \[Delta]3))/c;
ks2 = (\[Omega]s2 Sqrt[1 + \[Chi]S2])/c;
ks3 = (\[Omega]s3 Sqrt[1 + \[Chi]S3])/c;

\[CapitalDelta]k = ks1 - ks2 + ks3 - k1 + k2 - k3;

\[CapitalPhi] =
Sinc[\[CapitalDelta]k L/2] Exp[-I (ks1 + ks2 + ks3) L/2];

Plot3D[Abs[(L/
2 \[Pi]) NIntegrate[\[Kappa] \[CapitalPhi] Exp[-I (\[Delta]2 \
t2 + \[Delta]3 t3)], {\[Delta]2, -1 10^9, 1 10^9}, {\[Delta]3,
0, -1 10^9, 1 10^9}, WorkingPrecision -> 15,
PrecisionGoal -> 15]]^2, {t2, 0, 30 10^-9}, {t3, 0, 30 10^-9},
PlotRange -> All]

• Note that the complex exponential in the integrand has $\delta_2$ twice. The variables Den2 and Den3 are defined but not used. The limits of the $\delta_3$ integration are suspicious. Also, the first line Clear; does not do anything. Did you mean Clear["Global*"]? Commented Dec 10, 2021 at 3:18
• @LouisB I edited the code and corrected the mistakes. But it still doesn't generate a plot. The limits of integral should be infinity, but the function kappa is zero beyond 10^9. I have no idea why it doesn't give a plot. Commented Dec 10, 2021 at 14:19

The integrand takes small values, e.g.

\[Kappa] \[CapitalPhi] Exp[-I (\[Delta]2 t2 + \[Delta]3 t3)] /.
{\[Delta]2 -> 100, \[Delta]3 -> 100, t2 -> 100, t3 -> -100}


3.3303*10^-18 - 1.13036*10^-17 I

and ranges are big, these circumstances cause problems. Something can be done by rescaling t2 and t3 and other changes:

Plot3D[Abs[(L/2 \[Pi]) NIntegrate[10^21*\[Kappa] \[CapitalPhi] Exp[-I (\[Delta]2* t2 + \[Delta]3 *t3)] /.
{t2 -> t2*10^7, t3 -> t3*10^7}, {\[Delta]2, -1 10^8,1 10^8},
{\[Delta]3, 0, -1 10^8, 1 10^8}]]^2, {t2, 0, 1 }, {t3, 0, 1}]
`

It should be noticed that takes some time. The whole executed code on demand through Dropbox.