<?xml version="1.0" encoding="UTF-8"?>
<simulation xmds-version="2">
  <name>wigner</name>
  <author>Graham Dennis</author>
  <description>
    Simulation of the Wigner function for an anharmonic oscillator with the initial state
    being a coherent state.
    
    WARNING: This simulation will take a couple of hours.
  </description>
  
  <features>
    <benchmark />
    <!-- <error_check /> -->
    <bing />
    <fftw plan="patient" />

    <globals>
      <![CDATA[
      /* physical constants */
        const real hbar = 1.05457148e-34;
        
      /* system constants */
        const real omega = 0.0;
        const real Uint = 1.0*hbar;
        
      /* initial state constants */
        const real alpha_0 = 3.0;
        
      /* used for calculating constants derived from arguments */
        bool firstTime = true;
        
      /* absorbing boundary constants */
       const real absorb = 8.0;
       const real width = 0.3;
       
      /* derived constants */
        const complex miOmega = -i*omega;
        const complex miUint_hbar = -i*Uint/hbar;
        const real Uint_hbar = Uint/hbar;
      ]]>
    </globals>
  </features>
  
  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="x" lattice="128"  domain="(-6, 6)" />
      <dimension name="y" lattice="128"  domain="(-6, 6)" />
    </transverse_dimensions>
  </geometry>
  
  <vector name="main" initial_basis="x y" type="complex">
    <components> W </components>
    <initialisation>
      <![CDATA[
        W = 2.0/M_PI * exp(-2.0*(y*y + (x-alpha_0)*(x-alpha_0)));
      ]]>
    </initialisation>
  </vector>
  
  <vector name="dampConstants" initial_basis="x y" type="real">
    <components>damping</components>
    <initialisation>
      <![CDATA[
      real r = sqrt(x*x + y*y);
      
      if (r > _max_x-width)
        damping = 0.0;
      else
        damping = 1.0;
      ]]>
    </initialisation>
  </vector>
  
  <sequence>
    <integrate algorithm="ARK89" tolerance="1e-7" interval="7.0" steps="100000">
      <samples>200</samples>
      <operators>
        <operator kind="ex">
          <operator_names>Lx Ly Lxxx Lxxy Lxyy Lyyy</operator_names>
          <![CDATA[
            Lx = i*kx;
            Ly = i*ky;
            Lxxx = -i*kx*kx*kx;
            Lxxy = -i*kx*kx*ky;
            Lxyy = -i*kx*ky*ky;
            Lyyy = -i*ky*ky*ky;
          ]]>
        </operator>
        <integration_vectors>main</integration_vectors>
        <dependencies>dampConstants</dependencies>
        <![CDATA[
          dW_dt = damping * (
                  (omega + Uint_hbar*(-1.0 + x*x + y*y))*x*Ly[W]
                - (omega + Uint_hbar*(-1.0 + x*x + y*y))*y*Lx[W]
                - Uint_hbar/16.0*(x*(Lxxy[W] + Lyyy[W])
                                 -y*(Lxyy[W] + Lxxx[W]))
                );
        ]]>
      </operators>
    </integrate>
  </sequence>
  
  <output>
      <sampling_group basis="x y" initial_sample="yes">
        <moments>WR WI</moments>
        <dependencies>main</dependencies>
        <![CDATA[
          WR = W.Re();
          WI = W.Im();
        ]]>
      </sampling_group>
  </output>
</simulation>