EPR.xmds

Script source:
EPR.xmds.gz

<?xml version="1.0"?>
<!-- Einstein-Podolsky-Rosen (EPR) simulation in atom-molecular BECs -->

<simulation>
  
  <name>EPR_10_undeplete</name>
  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>3000</paths>
  <seed>1 2</seed>
  <noises>4</noises>

  <use_mpi>yes</use_mpi>

  <globals>
    <![CDATA[
    // this is a comment
    /* this is also
    a comment */
    const double mass=1.433e-25;
    const double hbar=1.0545887e-34;

    // const double mol_peak_dens=4.6296e7;
    const double mol_peak_dens=3.7037e7;
    const double x0=30*1e-6;
    const double d0=x0;
    const double t0=2*mass*x0*x0/hbar;
    const double xi0=1;

    // const double psi1_initial_0=10;
    // const double phi1_initial_0=10;
    const double psi1_initial_0=0;
    const double phi1_initial_0=0;


    const double mol_peak_dens_dimensionless=mol_peak_dens*d0;

    const double t1=8e-4;
    // const double Delta_s =-2e4 ;
    const double Delta_s =-2.2698e4 ;
    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=0;
    const double loss2=0;

    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;

    const double delta=t0*Delta_s;
    const double q_0=sqrt(fabs(delta));

    ]]>
  </globals>


  <field>
    <name>main</name>
    <dimensions>   x    </dimensions>
    <lattice>     1809   </lattice>
    <domains>    (-9,9) </domains>
    <samples>1 1 1</samples>
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>psi1 psi1beta psi2 phi2 phi2beta psi2beta phi1 phi1beta</components>
      <fourier_space>no</fourier_space>
      <![CDATA[

      double theta=1.0-x*x/(xi0*xi0) ;
      double shape=theta;

      double psi1_initial=psi1_initial_0*exp(-x*x/(2*xi0*xi0/4));
      double phi1_initial=phi1_initial_0*exp(-x*x/(2*xi0*xi0/4));

      if (theta<0) shape=0;

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


      psi1 = pcomplex(psi1_initial,-q_0*x);
      psi1beta = pcomplex(psi1_initial,-q_0*x);
      phi1 = pcomplex(phi1_initial,q_0*x);
      phi1beta = pcomplex(phi1_initial,q_0*x);

      psi2 = pcomplex(psi2_initial,0.);
      psi2beta = pcomplex(psi2beta_initial,0.);
      phi2 = pcomplex(psi2_initial,0.);
      phi2beta = pcomplex(psi2beta_initial,0.);


      ]]>

    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>8e-4</interval>
      <lattice>400</lattice>
      <samples>100 100 100</samples>
      <k_operators>
	<constant>yes</constant>
	<operator_names>L1 L1beta Lphi1 Lphi1beta L2 L2beta Lphi2 Lphi2beta</operator_names>

	<![CDATA[
	L1= rcomplex(-loss1/2,-delta-kx*kx);
	L1beta= rcomplex(-loss1/2,-delta-kx*kx);
	Lphi1= rcomplex(-loss1/2,delta+kx*kx);
	Lphi1beta= rcomplex(-loss1/2,delta+kx*kx);
	L2= rcomplex(-loss2/2,-kx*kx/2);
	L2beta= rcomplex(-loss2/2,-kx*kx/2);
	Lphi2= rcomplex(-loss2/2,kx*kx/2);
	Lphi2beta= rcomplex(-loss2/2,kx*kx/2);
	]]>

      </k_operators>
      
      <vectors>main</vectors>
      <![CDATA[

      dpsi1_dt = L1[psi1] -i*(u11*conj(psi1beta)*psi1)*psi1 + chi*psi2*conj(psi1beta) + c_sqrt(chi*psi2-i*u11*psi1*psi1)*(n_1);

      dpsi1beta_dt = L1beta[psi1beta] -i*(u11*conj(psi1)*psi1beta)*psi1beta + chi*psi2beta*conj(psi1) + c_sqrt(chi*psi2beta-i*u11*psi1beta*psi1beta)*(n_3);

      dpsi2_dt = L2[psi2] + 0*(n_2);
      dpsi2beta_dt = L2beta[psi2beta] + 0*(n_4);



      dphi1_dt = Lphi1[phi1] +i*(u11*conj(phi1beta)*phi1)*phi1 + chi*phi2*conj(phi1beta) + c_sqrt(chi*phi2+i*u11*phi1*phi1)*(n_1) ;

      dphi1beta_dt = Lphi1beta[phi1beta] +i*(u11*conj(phi1)*phi1beta)*phi1beta + chi*phi2beta*conj(phi1) + c_sqrt(chi*phi2beta+i*u11*phi1beta*phi1beta)*(n_3);

      dphi2_dt = Lphi2[phi2] + 0*(n_2);
      dphi2beta_dt = Lphi2beta[phi2beta] + 0*(n_4);

      ]]>
    </integrate>

    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>8e-7</interval>
      <lattice>2</lattice>
      <samples>2 2 2</samples>
      <k_operators>
	<operator_names>L1 L1beta Lphi1 Lphi1beta L2 L2beta Lphi2 Lphi2beta</operator_names>

	<![CDATA[
	double delta=t0*Delta_s;

	L1= rcomplex(-loss1/2,-delta-kx*kx);
	L1beta= rcomplex(-loss1/2,-delta-kx*kx);
	Lphi1= rcomplex(-loss1/2,delta+kx*kx);
	Lphi1beta= rcomplex(-loss1/2,delta+kx*kx);
	L2= rcomplex(-loss2/2,-kx*kx/2);
	L2beta= rcomplex(-loss2/2,-kx*kx/2);
	Lphi2= rcomplex(-loss2/2,kx*kx/2);
	Lphi2beta= rcomplex(-loss2/2,kx*kx/2);
	]]>
      </k_operators>
      <vectors>main</vectors>
      <![CDATA[
      dpsi1_dt = L1[psi1];
      dpsi1beta_dt = L1beta[psi1beta];
      dphi1_dt = Lphi1[phi1];
      dphi1beta_dt = Lphi1beta[phi1beta];
      dpsi2_dt = L2[psi2];
      dpsi2beta_dt = L2beta[psi2beta];
      dphi2_dt = Lphi2[phi2];
      dphi2beta_dt = Lphi2beta[phi2beta];
      ]]>

    </integrate>

  </sequence>


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

    <group>
      <sampling>
	<fourier_space> no</fourier_space>
	<lattice>       201</lattice>
	<moments>pow_dens pow_dens_alt pow_dens_alt2 pow_dens_mol pow_dens_mol_alt</moments>
	<![CDATA[
	pow_dens = conj(psi1beta)*psi1;
	pow_dens_alt = phi1beta*conj(phi1);
	pow_dens_alt2 = conj(psi1beta)*conj(phi1);
	pow_dens_mol = conj(psi2beta)*psi2;
	pow_dens_mol_alt = phi2beta*conj(phi2);
	]]>
      </sampling>
    </group>

    <group>
      <sampling>
	<fourier_space> yes</fourier_space>
	<lattice>       1809</lattice>
	<moments>pow_dens_k pow_dens_k_alt XkXk XminuskXminusk XkXminusk YkYk YminuskYminusk YkYminusk</moments>

	<![CDATA[
	pow_dens_k = conj(psi1beta)*psi1;
	pow_dens_k_alt = phi1beta*conj(phi1); 

	// pow_dens_mol_k = conj(psi2beta)*psi2;

	// akak = psi1*psi1;

	// adagkadagk = conj(psi1beta)*conj(psi1beta);
	// akaminusk = psi1*conj(phi1);
	// adagkadagminusk = conj(psi1beta)*phi1beta;
	// adagkaminusk = conj(psi1beta)*conj(phi1);
	// adagminuskak = phi1beta*psi1; 

	XkXk = (psi1+conj(psi1beta))*(psi1+conj(psi1beta));
	XminuskXminusk = (conj(phi1)+phi1beta)*(conj(phi1)+phi1beta);
	XkXminusk = (psi1+conj(psi1beta))*(conj(phi1)+phi1beta);

	YkYk = -(psi1-conj(psi1beta))*(psi1-conj(psi1beta));
	YminuskYminusk = -(conj(phi1)-phi1beta)*(conj(phi1)-phi1beta);
	YkYminusk = -(psi1-conj(psi1beta))*(conj(phi1)-phi1beta);

	]]>
      </sampling>
    </group>

    <group>
      <sampling>
	<fourier_space> no</fourier_space>
	<lattice>0</lattice>
	<moments>p p_mol p1 p2</moments>
	<![CDATA[
	p = conj(psi1beta)*psi1;
	p_mol=conj(psi2beta)*psi2;
	int weight=(x>0);
	p1=weight*p;
	p2=(1-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;
	]]>
      </post_propagation>
    </group>


  </output>


</simulation>

Generated by GNU enscript 1.6.3.