# Optimize regarding to one variable in a contour plot involving a complex integral

I want to do a contour plot of a function involving numerical integration of a fast oscillating complex exponential and optimize (maximize) the plot a every point regarding to the phase, which in the special case here is 3.25. Would be very thankful for any help here!

U = 1/2 +
I/4 (1/((x* y*(1 + x/y))/(y + x/y *x)) + 1/((x + x/y *y)/(
1 + x/y)));
P = -(1/16) (1/((x* y*(1 + x/y))/(y + x/y*x)) - 1/((x + x/y *y)/(
1 + x/y)))^2 ;

h = ContourPlot[
1/(4*((x + x/y *y)/(1 + x/y)) )*
NIntegrate[
Exp[I *3.25*(z1 - z2)]/((z1 - U) (z2 - Conjugate[U]) + P), {z1, 0,
1}, {z2, 0, 1}, Method -> "LocalAdaptive"] , {x, 1, 4}, {y, 1,
4}, PlotLegends -> Automatic, Contours -> 80]


As background: the plotted h is proportional to the efficiency of a photon conversion process in a nonlinear crystal and x and y are the crystal length L devided by the rayleigh ranges of the participating gaussian beams (strong laser beam + beam with a single photon in it). The phase is the relative phase of the beams which can practially always be freely choosen. The phase is different from zero to compensate for the gouy-phase shift near the beamwaist.

Guha, Shekhar, and Joel Falk. "The effects of focusing in the three‐frequency parametric upconverter." Journal of applied physics 51.1 (1980): 50-60.

Boyd, G. D., and D. A. Kleinman. "Parametric interaction of focused Gaussian light beams." Journal of Applied Physics 39.8 (1968): 3597-3639.

Lastzka, Nico, and Roman Schnabel. "The Gouy phase shift in nonlinear interactions of waves." Optics Express 15.12 (2007): 7211-7217.

• You show the contoutplot! What's your question? Jul 26 at 15:02

Do some simplifications, define nint[x,y,aa], Chop imaginary part, which is lower 10^-8.

U0 = 1/2 +
I/4 (1/((x*y*(1 + x/y))/(y + x/y*x)) +
1/((x + x/y*y)/(1 + x/y))) // Together;
P0 = -(1/16) (1/((x*y*(1 + x/y))/(y + x/y*x)) -
1/((x + x/y*y)/(1 + x/y)))^2 // Together;

cond = 1 < x < 4 && 1 < y < 4 && 0 < z1 < 1 && 0 < z2 < 1 && aa > 0;

U = U0 // ComplexExpand[#, TargetFunctions -> {Re, Im}] & //
FullSimplify[#, cond] & // Expand

P = P0 // ComplexExpand[#, TargetFunctions -> {Re, Im}] & //
FullSimplify[#, cond] &

igd[x_, y_, z1_, z2_, aa_] =
Exp[I*aa*(z1 - z2)]/((z1 - U) (z2 - Conjugate[U]) + P) //
FullSimplify[#, cond] &

(*   (64 E^(I aa (z1 - z2)))/(-((x - y)^4/(
x^2 y^2 (x + y)^2)) + (-4 - (3 I)/x - (3 I)/y + (4 I)/(x + y) +
8 z1) (-4 + (3 I)/x + (3 I)/y - (4 I)/(x + y) + 8 z2))   *)

nint[x_, y_, aa_] :=
NIntegrate[igd[x, y, z1, z2, aa], {z1, 0, 1}, {z2, 0, 1}] //
Chop[#, 10^-8] &

Manipulate[
ContourPlot[
1/(4*((x + x/y*y)/(1 + x/y)))*nint[x, y, aa], {x, 1, 4}, {y, 1, 4},
Contours -> 80, PerformanceGoal -> "Speed"], {{aa, 3.25}, 1, 5,
Appearance -> "Labeled"}, ContinuousAction -> False]


For high quality plot set: PerformanceGoal -> "Quality" , which is much slower.

Do you realy need 80 Contours? How about Contours -> Range[.4, 1.2, .025] ?

Edit

Answer to the comment of 'Felix Mann': "Originally I wanted to do one contour plot where at every point (x,y) the phase aa takes the value which maximizes the plotted value h. "

Let me give a quick solution with FindArgMax(think, there could be a more elegant way with some thinking about)

Define function to maximize h and find the corresponding aa with FindArgMax. (Using Re@h , to avoid problems due to little intermediate imaginary search points.)

ListContourPlot with not much points is faster.

h[x_, y_, aa_?NumericQ] :=
1/(4*((x + x/y*y)/(1 + x/y)))*nint[x, y, aa]

fam[x_, y_] := First@FindArgMax[{Re@h[x, y, aa], 1 < aa < 5}, aa]

Plot[Re@h[1.9, 2.8, aa], {aa, 1, 5},
Epilog -> {Red, PointSize[.03],
Point[{fam[1.9, 2.8], Re@h[1.9, 2.8, fam[1.9, 2.8]]}]},
GridLines -> {{fam[1.9, 2.8]}, Automatic}]

(tab = Flatten[Table[{x, y, fam[x, y]}, {x, 1, 4, .5}, {y, 1, 4, .5}],
1];) // Timing

ListPlot3D[tab]

ListContourPlot[tab, Contours -> Range[1, 4, .1]]


• Thank you so much! That is already very helpful for me. Originally I wanted to do one contour plot where at every point (x,y) the phase aa takes the value which maximizes the plotted value h. How could this be done? So instead of a global phase for all points (x,y) I want the phase to always take the value which maximizes h at that point Jul 27 at 13:13
• @Felix, these are details that should have been included in your question; please edit your question to include this. Jul 27 at 17:50
• Thank you so much for all your help! Somehow h takes values of about 3.7 in that new polts, but from theory it is bounded to h_max=1.068. It takes that value at x=y=2.84 for a phase of 3.25. That is due to plain waves having h=1 and gaussian beams being just slighlty more efficient. Jul 28 at 10:56
• @FelixMann you say, "Somehow h takes values of about 3.7 in that new plot, but ....." But what you see in the second plot, is not h, but the Max of aa at the given x,y point. Look at the defintion of fam. That is, what you wanted in you first comment: " Originally I wanted to do one contour plot where at every point (x,y) the phase aa takes the value which maximizes the plotted value h. " Please confirm step by step, waht every command is doing. Jul 28 at 18:41