(* Exact quantum and primitive semiclassical wave functions *) (* Before reading of this package, there should be defined (see examples below): en - energy m - mass v11, v22, v12 - diabatic potentials x0, x1 - initial and final points of integration, if x1=x1=0, then will be estimated atomatically xp0, xp1 - lower and upper bounds of the interval where results will be displayed *) (* After execution of this package, the results are: v11, v22, v12 - diabatic potentials wf1, wf2 - diabatic wave functions w1, w2, theta, dtheta - adiabatic potentials, rotation angle, and its derivative psi1abs, psi2abs - abs. values of primitive semiclassical w. f. psi1, psi2 - primitive semiclassical w. f. (complex-valued) psi1d, psi2d - primitive semiclassical w. f., diabatic representation Plots: pltv - plot of diabatic potentials pltd1, pltd2 - exact diabatic w. f. pltw - adiabatic potentials plta1, plta2 - exact adiabatic w. f. pltap1, pltap2 - abs. val. of primitive adiabatic w. f. pltac1, pltac2 - combined (exact/primitive) plots of adiabatic w. f. pltdp1, pltdp2 - abs. val. of primitive w. f. in diabatic representation pltdc1, pltdc2 - combined (exact/primitive) plots of diabatic w. f. *) Off[General::"spell", General::"spell1", FindMinimum::"fmgz", CompiledFunction::"cfse", CompiledFunction::"cfex"]; \ Get["C:\\sergeev\\tulane\\math\\packages\\argcolorplot.m"]; saturat = 0.2; (*0.1*) << Graphics`FilledPlot`; (* Automatic finding x0, x1 *) eps = 10.^-14; xstart = 0; dxstart = 0.1; fx = 1.1; dxmax = 10000; mdx = Log[dxmax/dxstart]/Log[fx]; xt = 0; dxtab = Table[xt = xt + dxstart fx^ndx, {ndx, 0, mdx}]; dxtab = Select[dxtab, # < dxmax &]; xtab = Sort[Join[xstart - dxtab, xstart + dxtab]]; mx = Length[xtab]; findx0x1 := ( Off[NMinimize::"cvmit"]; ftab = -Table[ xt0 = xtab[[nx]]; xt1 = xtab[[nx + 1]]; NMinimize[{-Abs[fun[x]], xt0 <= x <= xt1}, { x}, AccuracyGoal -> 4, PrecisionGoal -> 4] // First, {nx, mx - 1}]; On[NMinimize::"cvmit"]; fmax = Max[ftab]; Do[ If[ftab[[nx]] > eps fmax, Break[]]; x0 = xtab[[nx + 1]], {nx, mx - 1}]; Do[ If[ftab[[nx]] > eps fmax, Break[]]; x1 = xtab[[nx]], {nx, mx - 1, 1, -1}]; (*Print["x0 = ", x0, " x1 = ", x1];*) {x0, x1} ); If[x0==0 && x1==0, fun = v11'; {x0a, x1a} = findx0x1; fun = v22'; {x0b, x1b} = findx0x1; fun = v12; {x0c, x1c} = findx0x1; x0 = Min[x0a, x0b, x0c]; x1 = Max[x1a, x1b, x1c]; Print["eps = ", eps, " ", { {"", "x0", "x1"}, {"V11'", x0a, x1a}, {"V22'", x0b, x1b}, {"V12", x0c, x1c}, {"Total", x0, x1} } // TableForm]; ]; mpts = 22; {en, m, x0, x1} = N[{en, m, x0, x1}]; pts[mpts_] := Table[ N[x0 + npt(x1 - x0)/(mpts + 1)], {npt, mpts}]; (* Analyzing potentials *) (* Finding interval of function values *) intfun[fun_] := Interval[Table[ vtab = Table[ xt0 = x0 + (n - 1)(x1 - x0)/(mpts + 1); xt1 = x0 + (n + 1)(x1 - x0)/(mpts + 1); Minimize[{imm fun[x], xt0 <= x <= xt1}, {x}] // First, {n, mpts}]; Chop[imm Min[vtab]], {imm, 1, -1, -2}]]; int1 = intfun[v11]; int2 = intfun[v22]; int3 = intfun[v12]; Print[{{"Potential", "Range of change"}, {"v11", int1}, {"v22", int2}, {"v12", int3}} // TableForm]; (* Plot range in y - axes *) intu = IntervalUnion[int1, int2]; ymm = {Min[intu], Max[intu]}; yav = (Plus @@ ymm)/2; {ymin, ymax} = ymm = yav + 1.05(ymm - yav); If[ymax > 0 && 0 < ymin < 0.3 ymax, ymin = 0]; If[ymin < 0 && 0.3 ymin < ymax < 0, ymax = 0]; ymm = {ymin, ymax}; int = Interval[{ymin, ymax}]; yav = (Plus @@ ymm)/2; int0 = Interval[yav + 0.4(ymm - yav)]; (* Scaling of v12 *) {s, g} = {0, 0}; time = Timing[ Do[ s1 = isign Switch[ir1, 1, 1, 2, 2, 3, 5]10.^imant; int3s = s1 int3; inti = IntervalIntersection[int0, int3s]; g1 = Max[inti] - Min[inti]; intu = IntervalUnion[int0, int3s]; {minu, maxu} = {Min[intu], Max[intu]}; gu = maxu - minu; g1 = g1/gu; If[isign == -1, g1 = 0.8 g1]; If[g1 > g, s = s1; g = g1]; , {isign, -1, 1, 2}, {imant, -6, 6}, {ir1, 3}]; ] // First; intu = IntervalUnion[int, s int3]; ymm = {Min[intu], Max[intu]}; int = Interval[ymm]; yav = (Plus @@ ymm)/2; int0 = Interval[yav + 1.6(ymm - yav)]; If[IntervalMemberQ[int0, en], (* Incorporate E *) int = IntervalUnion[int, Interval[{en, en}]]; ymm = {Min[int], Max[int]}; yav = (Plus @@ ymm)/2; int = Interval[yav + 1.08(ymm - yav)]; ]; {yp0, yp1} = {Min[int], Max[int]}; ifen = IntervalMemberQ[int, en]; sstr = NumberForm[s, ExponentFunction -> (Null &)] // ToString; If[sstr === "1", sstr = ""]; pltv = If[ifen, Plot[{v11[x], v22[x], s v12[x], en}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.02, 0.02}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, (*FrameLabel -> {"x", "Diabat. pot. v11, v22, " <> sstr <> " v12 and E"},*) PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity], Plot[{v11[x], v22[x], s v12[x]}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.02, 0.02}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, (*FrameLabel -> {"x", "Diabat. pot. v11, v22, and " <> sstr <> " v12"},*) PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity] ]; pltv1 = If[ifen, Plot[{v11[x], v22[x], s v12[x], en}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.05, 0.05}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, (*FrameLabel -> {"x", "Diabat. pot. v11, v22, " <> sstr <> " v12 and E"},*) PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, FrameTicks -> None, DisplayFunction -> Identity], Plot[{v11[x], v22[x], s v12[x]}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.05, 0.05}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, (*FrameLabel -> {"x", "Diabat. pot. v11, v22, and " <> sstr <> " v12"},*) PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, FrameTicks -> None, DisplayFunction -> Identity] ]; slv := ( Clear[psi1, psi2]; s = NDSolve[{ -1/2/m psi1''[x] + (v11[x] - en)psi1[x] == -v12[x]psi2[x], -1/2/m psi2''[x] + (v22[x] - en)psi2[x] == -v12[x] psi1[x], psi1[x0] == psi10, psi1'[x0] == dpsi10, psi2[x0] == psi20, psi2'[x0] == dpsi20}, {psi1, psi2}, {x, x0, x1}, MaxSteps -> 1000000] // First; wf = {psi1, psi2} /. s; {psi11, dpsi11, psi21, dpsi21} = {psi1[x1], psi1'[x1], psi2[x1], psi2'[x1]} /. s; r = {psi11 + dpsi11/I/pd11, psi11 - dpsi11/I/pd11, psi21 + dpsi21/I/pd21, psi21 - dpsi21/I/pd21}/2; {wf, r} ); pd1[x_] = Sqrt[2m(en - v11[x])]; pd2[x_] = Sqrt[2m(en - v22[x])]; pd10 = pd1[x0]; pd20 = pd2[x0]; pd11 = pd1[x1]; pd21 = pd2[x1]; (* (+, 0) -> ... *) psi10 = Sqrt[m/pd10]; dpsi10 = I pd10 psi10; psi20 = 0; dpsi20 = 0; {{wf1i1, wf1i2}, r1i} = slv; Print["r1i = ", r1i // N]; (* (-, 0) -> ... *) psi10 = Sqrt[m/pd10]; dpsi10 = -I pd10 psi10; psi20 = 0; dpsi20 = 0; {{wf1r1, wf1r2}, r1r} = slv; Print["r1r = ", r1r // N]; (* (0, -) -> ... *) psi20 = Sqrt[m/pd20]; dpsi20 = -I pd20 psi20; psi10 = 0; dpsi10 = 0; {{wf2r1, wf2r2}, r2r} = slv; Print["r2r = ", r2r // N]; Clear[c1, c2]; r = r1i + c1 r1r + c2 r2r; swf = Solve[r[[{2, 4}]] == {0, 0}, {c1, c2}][[1]]; Print[swf]; {r1t, r2t} = r[[{1, 3}]] /. swf; {r1t, r2t} = {r1t, r2t}Sqrt[{pd11, pd21}/m]; {pp11, pp12} = Abs[{r1t, r2t}]^2; Print["P11 = ", pp11, " P12 = ", pp12, " P11+P12 = ", pp11 + pp12]; {rr11, rr12} = Abs[{c1, c2} /. swf]^2; Print["R11 = ", rr11, " R12 = ", rr12, " R11+R12 = ", rr11 + rr12]; Print["1-P11-P12-R11-R12 = ", 1 - pp11 - pp12 - rr11 - rr12]; (* Wave functions *) wf1[x_] := wf1i1[x] + c1 wf1r1[x] + c2 wf2r1[x] /. swf; wf2[x_] := wf1i2[x] + c1 wf1r2[x] + c2 wf2r2[x] /. swf; pltd1 = ArgColorPlot[wf1[x], {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}}, (*FrameLabel -> {"x", "Diabatic wave funcion-1"},*) PlotRange -> All, PlotPoints -> 99, Saturation -> saturat, Brightness -> 1, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity]; pltd2 = ArgColorPlot[wf2[x], {x, xp0, xp1}, PlotStyle -> {{RGBColor[0, 0, 1]}}, (*FrameLabel -> {"x", "Diabatic wave funcion-2"},*) PlotRange -> All, PlotPoints -> 99, Saturation -> saturat, Brightness -> 1, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity]; (* Adiabatic wavefunctions *) thet[x_] = (* Defined up to n/2 Pi *) 1/2 ArcTan[2 v12[x]/(v11[x] - v22[x])]; theta0 = thet[x0]; nthet0 = Round[2/Pi theta0]; theta0 = theta0 - Pi/2 nthet0; dtheta[x_] := ((v11[x] - v22[x])*Derivative[1][v12][x] + v12[x]* (-Derivative[1][v11][x] + Derivative[1][v22][x]))/(4* v12[x]^2 + (v11[x] - v22[x])^2); stheta = NDSolve[{thetax'[x] == dtheta[x], thetax[x0] == theta0}, {thetax}, {x, x0, x1}, MaxStepFraction -> 1/2000, MaxSteps -> 1000000] // First; theta[x_] = thetax[x] /. stheta; wf1a[x_] := Cos[theta[x]]wf1[x] + Sin[theta[x]]wf2[x]; wf2a[x_] := -Sin[theta[x]]wf1[x] + Cos[theta[x]]wf2[x]; w1[x_] := Cos[theta[x]]^2 v11[x] + Sin[theta[x]]^2 v22[x] + 2Cos[theta[x]]Sin[theta[x]] v12[x]; w2[x_] := Cos[theta[x]]^2 v22[x] + Sin[theta[x]]^2 v11[x] - 2Cos[theta[x]]Sin[theta[x]] v12[x]; pp12exact = If[Abs[Cos[theta[x1]]] < 0.5, pp11, pp12]; int1 = intfun[w1]; int2 = intfun[w2]; int3 = intfun[dtheta]; Print[{{"Potential", "Range of change"}, {"w1", int1}, {"w2", int2}, {"theta'", int3}} // TableForm]; enmin = Max[int1, int2]; (* Plot range in y - axes *) intu = IntervalUnion[int1, int2]; ymm = {Min[intu], Max[intu]}; yav = (Plus @@ ymm)/2; {ymin, ymax} = ymm = yav + 1.05(ymm - yav); If[ymax > 0 && 0 < ymin < 0.3 ymax, ymin = 0]; If[ymin < 0 && 0.3 ymin < ymax < 0, ymax = 0]; ymm = {ymin, ymax}; int = Interval[{ymin, ymax}]; yav = (Plus @@ ymm)/2; int0 = Interval[yav + 0.4(ymm - yav)]; (* Scaling of v12 *) {s, g} = {0, 0}; time = Timing[ Do[ s1 = isign Switch[ir1, 1, 1, 2, 2, 3, 5]10.^imant; int3s = s1 int3; inti = IntervalIntersection[int0, int3s]; g1 = Max[inti] - Min[inti]; intu = IntervalUnion[int0, int3s]; {minu, maxu} = {Min[intu], Max[intu]}; gu = maxu - minu; g1 = g1/gu; If[isign == -1, g1 = 0.8 g1]; If[g1 > g, s = s1; g = g1]; , {isign, -1, 1, 2}, {imant, -6, 6}, {ir1, 3}]; ] // First; intu = IntervalUnion[int, s int3]; ymm = {Min[intu], Max[intu]}; int = Interval[ymm]; yav = (Plus @@ ymm)/2; int0 = Interval[yav + 1.6(ymm - yav)]; If[IntervalMemberQ[int0, en], (* Incorporate E *) int = IntervalUnion[int, Interval[{en, en}]]; ymm = {Min[int], Max[int]}; yav = (Plus @@ ymm)/2; int = Interval[yav + 1.08(ymm - yav)]; ]; {yp0, yp1} = {Min[int], Max[int]}; ifen = IntervalMemberQ[int, en]; sstr = NumberForm[s, ExponentFunction -> (Null &)] // ToString; If[sstr === "1", sstr = ""]; sv12 = s; pltw = If[ifen, Plot[{w1[x], w2[x], s dtheta[x], en}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.02, 0.02}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity], Plot[{w1[x], w2[x], s dtheta[x]}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.02, 0.02}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity] ]; pltw1 = If[ifen, Plot[{w1[x], w2[x], s dtheta[x], en}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.05, 0.05}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, FrameTicks -> None, DisplayFunction -> Identity], Plot[{w1[x], w2[x], s dtheta[x]}, {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}, {RGBColor[0.5, 0.5, 0.5]}, {Dashing[{0.05, 0.05}]}}, PlotRange -> {{xp0, xp1}, {yp0, yp1}}, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, FrameTicks -> None, DisplayFunction -> Identity] ]; gr1 = GraphicsArray[{{pltv1, pltw1}}]; Show[gr1, ImageSize -> 200, DisplayFunction -> $DisplayFunction]; plta1 = ArgColorPlot[wf1a[x], {x, xp0, xp1}, PlotStyle -> {{RGBColor[1, 0, 0]}}, (*FrameLabel -> {"x", "Diabatic wave funcion-1"},*) PlotRange -> All, PlotPoints -> 99, Saturation -> saturat, Brightness -> 1, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity]; plta2 = ArgColorPlot[wf2a[x], {x, xp0, xp1}, PlotStyle -> {{RGBColor[0, 0, 1]}}, (*FrameLabel -> {"x", "Diabatic wave funcion-2"},*) PlotRange -> All, PlotPoints -> 99, Saturation -> saturat, Brightness -> 1, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity]; (* Primitive semiclassical *) pa1[x_] = Sqrt[2m(en - w1[x])]; pa2[x_] = Sqrt[2m(en - w2[x])]; pa10 = pa1[x0]; pa20 = pa2[x0]; pa11 = pa1[x1]; pa21 = pa2[x1]; pfact[x_] = ( p1x = pa1[x]; p2x = pa2[x]; (p1x + p2x)/2/Sqrt[p1x p2x] ); tau[x_] = pfact[x] dtheta[x]; p12[x_] = (pa1[x] - pa2[x])/2; eq1 = phi'[x] - 2(p12[x] - tau[x]Sin[phi[x]]Cot[2 bet[x]]); eq2 = bet'[x] - tau[x]Cos[phi[x]]; time = Timing[ s = NDSolve[{eq1 == 0, eq2 == 0, phi[x0] == 0, bet[x0] == 10.^(-10)}, {phi, bet}, {x, x0, x1}, MaxSteps -> 100000, AccuracyGoal -> 12, PrecisionGoal -> 12][[1]]; ][[1]]; Print["Time of calculation of beta and phim: ", time]; phim[x_] = phi[x] /. s; beta[x_] = bet[x] /. s; pp12 = Sin[beta[x1]]^2; psi1abs[x_] = Sqrt[m/pa1[x]]Abs[Cos[beta[x]]]; psi2abs[x_] = Sqrt[m/pa2[x]]Abs[Sin[beta[x]]]; (* Plotting *) Print["Plotting primitive WKB wave functions ..."]; pltap1 = Plot[psi1abs[x], {x, xp0, xp1}, PlotStyle -> {RGBColor[0.3, 0, 0], Dashing[{0.02, 0.02}]}, PlotPoints -> 99, PlotRange -> All, DisplayFunction -> Identity]; pltap2 = Plot[psi2abs[x], {x, xp0, xp1}, PlotStyle -> {RGBColor[0, 0, 0.3], Dashing[{0.02, 0.02}]}, PlotPoints -> 99, PlotRange -> All, DisplayFunction -> Identity]; Print["P12 = ", pp12]; (* finding phip *) time = Timing[ s1 = NDSolve[{php'[x] == 2 tau[x] Sin[phim[x]]/Sin[2 beta[x]], php[x0] == 0}, php, {x, x0, x1}, MaxSteps -> 100000, AccuracyGoal -> 12, PrecisionGoal -> 12][[1]]; ][[1]]; Print["Time of calculation of phip: ", time]; phip[x_] = php[x] /. s1; phia[x_] = (phip[x] + phim[x])/2; phib[x_] = (phip[x] - phim[x])/2; (* finding phi0 *) time = Timing[ s2 = NDSolve[{ph0'[x] == (pa1[x] + pa2[x])/2 , ph0[x0] == 0}, ph0, {x, x0, x1}, MaxSteps -> 100000, AccuracyGoal -> 12, PrecisionGoal -> 12][[1]]; ][[1]]; Print["Time of calculation of phi0: ", time]; phi0[x_] = ph0[x] /. s2; psi1[x_] = Exp[I phi0[x]]Sqrt[m/pa1[x]]Exp[I phia[x]]Cos[beta[x]]; psi2[x_] = -Exp[I phi0[x]]Sqrt[m/pa2[x]]Exp[I phib[x]]Sin[beta[x]]; Print["Plotting primitive WKB wave functions in diabatic representation"]; psi1d[x_] := Cos[theta[x]]psi1[x] - Sin[theta[x]]psi2[x]; psi2d[x_] := Sin[theta[x]]psi1[x] + Cos[theta[x]]psi2[x]; pltdp1 = Plot[Abs[psi1d[x]], {x, xp0, xp1}, PlotStyle -> {RGBColor[0.3, 0, 0], Dashing[{0.02, 0.02}]}, PlotPoints -> 99, PlotRange -> All, DisplayFunction -> Identity]; pltdp2 = Plot[Abs[psi2d[x]], {x, xp0, xp1}, PlotStyle -> {RGBColor[0, 0, 0.3], Dashing[{0.02, 0.02}]}, PlotPoints -> 99, PlotRange -> All, DisplayFunction -> Identity]; dp1[x_] = pa1'[x]; dp2[x_] = pa2'[x]; (*** ** ** ** ** ** ********) (* finding phases *) (* Combining plots of exact and primitive w. functions *) pltac1 = Show[plta1, pltap1, DisplayFunction -> Identity]; pltac2 = Show[plta2, pltap2, DisplayFunction -> Identity]; pltdc1 = Show[pltd1, pltdp1, DisplayFunction -> Identity]; pltdc2 = Show[pltd2, pltdp2, DisplayFunction -> Identity]; Print["------------\n"]; Print[name]; Print["Energy = ", en // N, " Mass = ", m // N, " P12 = ", pp12exact, " (", pp12, ")"]; Print["V(x) = ", {{v11[x], v12[x]}, {v12[x], v22[x]}} // MatrixForm]; psi1g[x3_Real, gamma_Real] := ( eps = 10.^(-32); dx = Sqrt[-Log[eps]/gamma]; If[iprint, Print["psi1g: ", x3]]; xi0 = Max[x0, x3 - dx]; xi1 = Min[x1, x3 + dx]; s1a = NDSolve[{ps1a'[x] == Sqrt[gamma/Pi(1 + I dp1[x]/2/gamma)]*psi1[x]* Exp[-gamma(x3 - x)^2 + I pa1[x](x3 - x)], ps1a[xi0] == 0}, ps1a, {x, xi0, xi1}, MaxSteps -> 100000][[1]]; ps1a[xi1] /. s1a); psi2g[x3_Real, gamma_Real] := ( eps = 10.^(-32); dx = Sqrt[-Log[eps]/gamma]; If[iprint, Print["psi2g: ", x3]]; xi0 = Max[x0, x3 - dx]; xi1 = Min[x1, x3 + dx]; s2a = NDSolve[{ps2a'[x] == Sqrt[gamma/Pi(1 + I dp2[x]/2/ gamma)]*psi2[x]Exp[-gamma(x3 - x)^2 + I * pa2[x](x3 - x)], ps2a[xi0] == 0}, ps2a, {x, xi0, xi1}, MaxSteps -> 100000][[1]]; ps2a[xi1] /. s2a); plotgamma[gamma_] := ( label0 = "\[Gamma] = " <> ToString[NumberForm[gamma, ExponentFunction -> (Null &)]]; box = "\!\(\* StyleBox[ FrameBox[StyleBox[\"" <> label0 <> "\", FontSize->12, FontWeight->\"Bold\", FontColor->RGBColor[1,0,0]], BoxMargins->{{0.2, 0.2}, {0.6, 0.6}}], FontSize->12, FontWeight->\"Plain\", FontColor->GrayLevel[0.700008], Background->GrayLevel[1], FontVariations->{\"CompatibilityType\"->0}]\)"; label1 = Graphics[{Text[box, {0.8, 0.9} // Scaled]}]; box = "\!\(\* StyleBox[ FrameBox[StyleBox[\"" <> label0 <> "\", FontSize->12, FontWeight->\"Bold\", FontColor->RGBColor[0,0,1]], BoxMargins->{{0.2, 0.2}, {0.6, 0.6}}], FontSize->12, FontWeight->\"Plain\", FontColor->GrayLevel[0.700008], Background->GrayLevel[1], FontVariations->{\"CompatibilityType\"->0}]\)"; label2 = Graphics[{Text[box, {0.8, 0.9} // Scaled]}]; pltg1 = FilledPlot[{psi1[x3] // Abs, psi1g[x3 // N, gamma // N] // Abs}, {x3, xp0, xp1}, Fills -> {RGBColor[0.8, 1, 0.8]}, PlotRange -> All, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity, PlotStyle -> {{Dashing[{0.02, 0.02}]}, {RGBColor[1, 0, 0]}}]; pltg2 = FilledPlot[{psi2[x3] // Abs, psi2g[x3 // N, gamma // N] // Abs}, {x3, xp0, xp1}, Fills -> {RGBColor[0.8, 1, 0.8]}, PlotRange -> All, PlotPoints -> 99, Frame -> True, Axes -> True, Ticks -> None, DisplayFunction -> Identity, PlotStyle -> {{Dashing[{0.02, 0.02}]}, {RGBColor[0, 0, 1]}}]; pltg1 = Show[pltg1, label1, DisplayFunction -> Identity]; pltg2 = Show[pltg2, label2, DisplayFunction -> Identity]; {pltg1, pltg2}); (* Exact and primitive semiclassical probabilities *) p12exact[en_Real] := ( pd1[x_] = Sqrt[2m(en - v11[x])]; pd2[x_] = Sqrt[2m(en - v22[x])]; pd10 = pd1[x0]; pd20 = pd2[x0]; pd11 = pd1[x1]; pd21 = pd2[x1]; slv := ( Clear[psi1, psi2]; s = NDSolve[{ -1/2/m psi1''[x] + (v11[x] - en)psi1[x] == -v12[x]psi2[x], -1/2/m psi2''[x] + (v22[x] - en)psi2[x] == -v12[x] psi1[x], psi1[x0] == psi10, psi1'[x0] == dpsi10, psi2[x0] == psi20, psi2'[x0] == dpsi20}, {psi1, psi2}, {x, x0, x1}, MaxSteps -> 1000000] // First; wf = {psi1, psi2} /. s; {psi11, dpsi11, psi21, dpsi21} = {psi1[x1], psi1'[x1], psi2[x1], psi2'[x1]} /. s; r = {psi11 + dpsi11/I/pd11, psi11 - dpsi11/I/pd11, psi21 + dpsi21/I/pd21, psi21 - dpsi21/I/pd21}/2; {wf, r} ); (*(+, 0) -> ...*) psi10 = Sqrt[m/pd10]; dpsi10 = I pd10 psi10; psi20 = 0; dpsi20 = 0; {{wf1i1, wf1i2}, r1i} = slv; (*(-, 0) -> ...*) psi10 = Sqrt[m/pd10]; dpsi10 = -I pd10 psi10; psi20 = 0; dpsi20 = 0; {{wf1r1, wf1r2}, r1r} = slv; (*(0, -) -> ...*) psi20 = Sqrt[m/pd20]; dpsi20 = -I pd20 psi20; psi10 = 0; dpsi10 = 0; {{wf2r1, wf2r2}, r2r} = slv; Clear[c1, c2]; r = r1i + c1 r1r + c2 r2r; swf = Solve[r[[{2, 4}]] == {0, 0}, {c1, c2}][[1]]; {r1t, r2t} = r[[{1, 3}]] /. swf; {r1t, r2t} = {r1t, r2t}Sqrt[{pd11, pd21}/m]; {pp11, pp12} = Abs[{r1t, r2t}]^2; {rr11, rr12} = Abs[{c1, c2} /. swf]^2; If[Abs[Cos[theta[x1]]] < 0.5, pp11, pp12]); p12sc[en_Real] := ( pa1[x_] = Sqrt[2m(en - w1[x])]; pa2[x_] = Sqrt[2m(en - w2[x])]; pa10 = pa1[x0]; pa20 = pa2[x0]; pa11 = pa1[x1]; pa21 = pa2[x1]; pfact[x_] = (p1x = pa1[x]; p2x = pa2[x]; (p1x + p2x)/2/Sqrt[p1x p2x]); tau[x_] = pfact[x] dtheta[x]; p12[x_] = (pa1[x] - pa2[x])/2; eq1 = phi'[x] - 2(p12[x] - tau[x]Sin[phi[x]]Cot[2 bet[x]]); eq2 = bet'[x] - tau[x]Cos[phi[x]]; s = NDSolve[{eq1 == 0, eq2 == 0, phi[x0] == 0, bet[x0] == 10.^(-10)}, {phi, bet}, {x, x0, x1}, MaxSteps -> 100000, AccuracyGoal -> 12, PrecisionGoal -> 12][[1]]; phim[x_] = phi[x] /. s; beta[x_] = bet[x] /. s; pp12 = Sin[beta[x1]]^2; pp12 );