Twin_beams.xmds

<?xml version="1.0"?>
<!--Twin_beams 1D simulation in xmds version 1.2-->

<simulation>

  <name>Twin_beams</name>

  <author>Karen Kheruntsyan</author>
  <description>
    Twin beams 1D simulation
  </description>

  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>40000</paths>
  <seed>1 2</seed>
  <noises>4</noises>
  <use_mpi>yes</use_mpi>

  <globals>
  <![CDATA[
    const double mass = 1.433e-25;
    const double hbar = 1.0545887e-34;
    const double mol_peak_dens = 4.6296e7;
    const double x0 = 30*1e-6;
    const double d0 = x0;
    const double t0 = 2*mass*x0*x0/hbar;
    const double xi0 = 1;

    const double mol_peak_dens_dimensionless = mol_peak_dens*d0;
    
    const double t1 = 8e-4;
    
    const double Delta_s = -2e4;
    const double chi_1D = 0.1881;     
    const double U_aa_1D = 4.4156e-5;
    const double U_am_1D = -5.6703e-5;
    const double U_mm_1D = 2.2078e-5;

    const double loss1 = 10;
    const double loss2 = 0;
    
    const double delta = t0*Delta_s;
    const double chi_dimensionless = t0*chi_1D/sqrt(d0);
    const double u11 = 0;
    const double u12 = 0;
    const double u22 = t0*U_mm_1D/d0;
    
    const double chi = chi_dimensionless;
  ]]>
  </globals>

  <field>
    <name>main</name>
    <dimensions>   x    </dimensions>
    <lattice>     2000   </lattice>
    <domains>    (-10,10) </domains>
    <samples>1 1</samples>
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>psi1 psi1dag psi2 psi2dag</components>
      <fourier_space>no</fourier_space>
      <![CDATA[
        double theta=1.0-x*x/(xi0*xi0) ;
        double shape=theta;

        if (theta<0) shape=0;

        double psi2_initial=sqrt(mol_peak_dens_dimensionless*shape);
        double psi2dag_initial=sqrt(mol_peak_dens_dimensionless*shape);

        psi1 = pcomplex(0.,0.);
        psi1dag = pcomplex(0.,0.);
        psi2 = pcomplex(psi2_initial,0.);
        psi2dag = pcomplex(psi2dag_initial,0.);
      ]]>
    </vector>
 </field>
 
 <sequence>
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>8e-4</interval>
      <lattice>400</lattice>
      <samples>10 10</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names>Lpsi1 Lpsi1dag Lpsi2 Lpsi2dag</operator_names>
        <![CDATA[
          Lpsi1 = rcomplex(-loss1/2,-delta-kx*kx);
          Lpsi1dag = rcomplex(-loss1/2,delta+kx*kx);
          Lpsi2 = rcomplex(-loss2/2,-kx*kx/2);
          Lpsi2dag = rcomplex(-loss2/2,kx*kx/2);
        ]]>
      </k_operators>
      <vectors>main</vectors>
      <![CDATA[
        dpsi1_dt = Lpsi1[psi1] -i*(u11*psi1dag*psi1)*psi1 +
          chi*psi2*psi1dag + c_sqrt(chi*psi2-i*u11*psi1*psi1)*(n_1);
        dpsi2_dt = Lpsi1dag[psi1dag] -i*(u22*psi2dag*psi2 
          +u22*mol_peak_dens_dimensionless*x*x/(xi0*xi0))*psi2 -
          (chi/2)*psi1*psi1 + c_sqrt(-i*u22*psi2*psi2)*(n_2);

        dpsi1dag_dt = Lpsi2[psi2] +i*(u11*psi1*psi1dag)*psi1dag +
          chi*psi2dag*psi1 +
          c_sqrt(chi*psi2dag+i*u11*psi1dag*psi1dag)*(n_3);
        dpsi2dag_dt = Lpsi2dag[psi2dag] +i*(u22*psi2*psi2dag
          +u22*mol_peak_dens_dimensionless*x*x/(xi0*xi0))*psi2dag -
          (chi/2)*psi1dag*psi1dag + c_sqrt(i*u22*psi2dag*psi2dag)*(n_4);
      ]]>
    </integrate>

    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>0.0102</interval>
      <lattice>100</lattice>
      <samples>100 100</samples>
      <k_operators>
        <operator_names>Lpsi1 Lpsi1dag Lpsi2 Lpsi2dag</operator_names>
        <![CDATA[
          Lpsi1 = rcomplex(-loss1/2,-delta-kx*kx);
          Lpsi1dag = rcomplex(-loss1/2,delta+kx*kx);
          Lpsi2 = rcomplex(-loss2/2,-kx*kx/2);
          Lpsi2dag = rcomplex(-loss2/2,kx*kx/2);
        ]]>
      </k_operators>
      <vectors>main</vectors>
      <![CDATA[
        dpsi1_dt = Lpsi1[psi1];
        dpsi2_dt = Lpsi1dag[psi1dag];
        
        dpsi1dag_dt = Lpsi2[psi2];
        dpsi2dag_dt = Lpsi2dag[psi2dag];
      ]]>
    </integrate>
  </sequence>

  <output>
    <filename>Twin_beams.xsil</filename>
    <overwrite>yes</overwrite>

    <group>
      <sampling>
        <fourier_space> no</fourier_space>
        <lattice>       100</lattice>
        <moments>pow_dens pow_dens_mol</moments>
        <![CDATA[
          pow_dens = psi1dag*psi1;
          pow_dens_mol = psi2dag*psi2;
        ]]>
      </sampling>
    </group>

    <group>
      <sampling>
        <fourier_space> no</fourier_space>
        <lattice>0</lattice>
        <moments>p p_mol p1 p2</moments>
        <![CDATA[
          p = psi1dag*psi1;
          p_mol = psi2dag*psi2;
          int weight = (x>0.0);
          p1 = weight*p;
          p2 = (1.0-weight)*p;
        ]]>
      </sampling>
        
      <post_propagation>
        <fourier_space> no</fourier_space>
        <moments>p_av p_mol_av p1_av p2_av p1p1_av p2p2_av p1p2_av var var_av</moments>
        <![CDATA[
          p_av = real(p);
          p_mol_av = real(p_mol);
          p1_av = real(p1);
          p2_av = real(p2);
          p1p1_av = real(p1*p1);
          p2p2_av = real(p2*p2);
          p1p2_av = real(p1*p2);
          var = real(p1*p1+p2*p2-2*p1*p2);
          var_av = p1p1_av+p2p2_av-2*p1p2_av;
          // the normalized variance in the particle number difference fluctuations
          // should be calculated (in the matlab code) via
          // variance_norm=1+var_av/p_av
        ]]>
      </post_propagation>
    </group>
  </output>

</simulation>