<?xml version="1.0"?>
<!-- Example grand-canonical, imaginary time calculation for an -->
<!-- anharmonic oscilator using the thermal representation with gauges -->

<simulation>

  <name>canonical_bose</name>

  <author>Joel Corney</author>
  <description>
    Example grand-canonical, imaginary time calculation for an
    anharmonic oscilator using the thermal representation with gauges
  </description>

  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>1000</paths>
  <use_mpi>no</use_mpi>
  <seed>1 2</seed>
  <noises>1</noises>

  /***********/
  <globals>
  /***********/

  <![CDATA[
    const double Pi = 2.0*asin(1.0);
    const double N = 10.0;        /*Number of initial particles*/
    const double kappa = 0.1; /*Interaction strength*/
    const double omega = 1;  /*Oscillator strength*/
    ]]>
  </globals>

  /*************/
  <field>
  /*************/
    <name>main</name>
    <samples>1</samples>
 
    /***MAIN***/
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>n Omega</components>
      <fourier_space>no</fourier_space>    
      <![CDATA[
        /* NB: Initial conditions */
        n = rcomplex(N,0.0);
        Omega = 1.0;
        ]]>
    </vector>        
  </field>

  /******************/
  <sequence>
  /******************/
 
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>5</interval>
      <lattice>1000</lattice>
      <samples>200</samples>
      <iterations>4</iterations>
      <vectors>main</vectors>
      <![CDATA[
        /*correlations */
        double x = real(n);
        double y = imag(n);
        complex bose = n*(1+n);  /*Bose enhancement factor*/
        double X = real(bose);
        double Y = imag(bose);
        double Theta = X<0;
        complex K1 = c_sqrt(-rcomplex(2.0,0.0)*kappa);
        complex g = -4*K1*Theta*X/(1+n);
        complex strat = 4*Theta *Omega*kappa*(-i*(1+2*x)*Y + (~n + 4*X/(1+n))*X)/(1+n);

        complex s_1 = K1*bose;  /*Noise Correlations*/
        complex s_0 = Omega*g;  /*Noise Correlations*/

        dOmega_dt = -Omega*(omega*n +kappa*2*n*n) + strat + s_0*n_1;
        dn_dt = -(omega*n + kappa*n*(4*(1-2*Theta)*X + 4*i*Y -(1+n)*(2*n+1))) + s_1*n_1;
      ]]>
    </integrate>
  </sequence>

  /**************/
  <output>
    <filename>canonical_bose.xsil</filename>
    <group>
      <sampling>
       <vectors>main</vectors>
        <moments>x y thet N NN Weight</moments>
        <![CDATA[
          double xx = real(n);
          double yy = imag(n);
          x = xx;                   /*real part of unnormalised mean particle number*/
          y = yy;                   /*imag part of unnormalised mean particle number*/

          thet = xx*xx<yy*yy; /*=one when gauge is active, =zero otherwise*/
          N = n*Omega;      /*Weighted mean number of particles*/
          NN = n*n*Omega;           /*Weighted mean number of particles*/
          Weight = Omega;           /*Mean Weight*/
        ]]>
      </sampling>
    </group>
   </output>
</simulation>
