(* Wigner function for a superposition of two Gaussians *)
Clear["@"];
xend = 8;
pend = 1.5;
scpsi = 2;
scrho = 4;
npoints = 100;
{amin, amax, astep} = {0, 10, 1};
SetDirectory["D:\\temp"];
fnames = " ";
batchfile = "conv.bat";
delay = 0.5;
delay100 = Round[delay 100];
animfile = "twogauss.gif";
(* General formula, for any a *)
psi[x_] = E^(-(x - a/2)^2/2) + E^(-(x + a/2)^2/2);
Print["Psi = ", psi[x]];
npsi = Simplify[Integrate[psi[x]^2, {x, -Infinity, Infinity}]];
Print["N = ", npsi];
v[x_] = Simplify[psi''[x]/2/psi[x]];
Print["V = ", v[x]];
rho[x1_, p1_] =
    Simplify[1/Pi Integrate[psi[x1 + x2] psi[x1 - x2]/npsi Exp[-2 I x2 p1],
          {x2, -Infinity, Infinity}]];
Print["Rho = ", rho[x, p]];
(***)
nt = 0;
Do[nt++;
    pltpsi = Plot[scpsi psi[x]/Sqrt[npsi], {x, -xend, xend},
        PlotStyle -> {Thickness[0.007], RGBColor[1, 0, 0]}, 
        DisplayFunction -> Identity];
    pltv = Plot[v[x], {x, -xend, xend},
        PlotStyle -> {Thickness[0.01]}, DisplayFunction -> Identity];
    pltrho = Plot[scrho rho[x1, 0], {x1, -xend, xend},
        PlotStyle -> {Thickness[0.01], RGBColor[0, 1, 0]}, 
        DisplayFunction -> Identity];
    $TextStyle = {FontFamily -> "Times", FontSize -> 18};
    pltup = Show[pltv, pltpsi, pltrho,
        AxesLabel -> {"\!\(\* StyleBox[\"q\", FontSlant->\"Italic\"]\)",
            "\!\(\* StyleBox[\"V\", FontSlant->\"Italic\"]\), " <>
              ToString[scpsi] <> 
              "\[Times]\!\(\* StyleBox[\"\[Psi]\",FontColor->RGBColor[1, 0, 0]]\), " <>
              ToString[scrho] <> 
              "\[Times]\!\(\* StyleBox[\"\[Rho]\",FontColor->RGBColor[0, 1, 0]]\)"},
        DisplayFunction -> Identity,
        PlotRange -> {{-xend, xend}, {-0.6, 1.5}},
        Frame -> True,
        AspectRatio -> 1/4, ImageSize -> 600];
    If[a == 0, zmin0 = 0,
      {zmin0, fmin} = FindMinimum[rho[0, p], {p, 0.001}]];
    zmax0 = rho[0, 0];
    zscale = Max[Abs[zmin0], Abs[zmax0]];
    (* Density plot, takes some time depending on value of npoints *)
    {xmin, xmax} = 0.635{-xend, xend};
    {ymin, ymax} = {-pend, pend};
    z1max = 1.3;
    z1white = 0.7;
    cf = (z = zmin0 + #(zmax0 - zmin0);
          z1 = Abs[z/zscale] z1max;
          
          Hue[If[z > 0, 0, 0.65], 
            If[z1 <= 1, 0.25, 0.25(z1max - z1white z1)/(z1max - z1white)], 
            If[z1 <= 1, z1, 1]]) &;
    pltdown = ContourPlot[rho[x, y]/zscale, {x, xmin, xmax}, {y, ymin, ymax},
        PlotPoints -> npoints + 1,
        Contours -> 200,
        PlotRange -> {zmin0, zmax0}/zscale,
        FrameLabel -> {"\!\(\* StyleBox[\"q\", FontSlant->\"Italic\"]\)", 
            "\!\(\* StyleBox[\"p\", FontSlant->\"Italic\"]\)"},
        ContourLines -> False,
        DisplayFunction -> Identity,
        AspectRatio -> 1/2.5,
        ImageSize -> 600,
        PlotRegion -> {{-0.056, 1}, {0, 1}},
        ColorFunction -> cf];
    comb = GraphicsArray[{{pltup}, {pltdown}}, GraphicsSpacing -> -0.03];
    plt = Show[Graphics[comb],
        ImageSize -> 900,
        DisplayFunction -> $DisplayFunction];
    pictfile = ToString[nt] <> ".bmp";
    fnames = fnames <> pictfile <> " ";
    Display[pictfile, plt, "BMP"], {a, amin, amax, astep}];
os = OpenWrite[batchfile, PageWidth -> Infinity, FormatType -> OutputForm];
Write[os, 
    "D:\\programs\\ImageMagick-5.1.1\\convert -adjoin -loop 4096 -delay " <> 
      ToString[delay100] <> fnames <> animfile];
Close[os];
Run[batchfile];