gaugeatomlaser.xmds

XMDS is a code generator that integrates equations. You write them down in human readable form in a

<?xmds version="1.0"?>
<!--Stochastic Atom laser simulation-->

<!-- $Id: gaugeatomlaser.xmds,v 1.13 2004/07/14 06:37:05 joehope Exp $ -->

<!--  Copyright (C) 2000-2004                                           -->
<!--                                                                    -->
<!--  Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane -->
<!--                                                                    -->
<!--  This file is part of xmds.                                        -->
<!--                                                                    -->
<!--  This program is free software; you can redistribute it and/or     -->
<!--  modify it under the terms of the GNU General Public License       -->
<!--  as published by the Free Software Foundation; either version 2    -->
<!--  of the License, or (at your option) any later version.            -->
<!--                                                                    -->
<!--  This program is distributed in the hope that it will be useful,   -->
<!--  but WITHOUT ANY WARRANTY; without even the implied warranty of    -->
<!--  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     -->
<!--  GNU General Public License for more details.                      -->
<!--                                                                    -->
<!--  You should have received a copy of the GNU General Public License -->
<!--  along with this program; if not, write to the Free Software       -->
<!--  Foundation, Inc., 59 Temple Place - Suite 330, Boston,            -->
<!--  MA  02111-1307, USA.                                              -->

 <simulation>

    <!-- Global system parameters and functionality -->
    <name>gaugeatomlaser</name> 

    <author>Unknown Author</author>
    <description>
      Gauge atom laser simulation
    </description>

    <stochastic>yes</stochastic> 
    <prop_dim>t</prop_dim>
    <error_check>yes</error_check>
    <paths>1</paths>
    <use_mpi>no</use_mpi> 
    <seed>1 7</seed>
    <noises>11</noises>
    <use_wisdom>yes</use_wisdom>

    <!-- Global variables for the simulation -->
    <globals>
      <![CDATA[
        const double noise = 1.0;
        /* physical constants  */
        const double r=5.0e3;
        const double rhos=4.0e7;
        const double omegax = 25;
        const double transversearea = 1.2e-11;
        const double g = 0*9.8;
        const double kappamax = 1.0e1;
        const double kick = -1.0e7;
        const double Utt = 0.005*1.96e-50;
        const double Utu = 0*2.9e-51;
        const double Uuu = 0.005*2.974797272874263e-51;
        const double inum = 375;
        const double losst1 = 7.0e-3;
        const double losst2 = 2.0e-19/transversearea;
        const double lossu1 = 7.0e-3;
        const double lossu2 = 4e-20/transversearea;
        const double losstu2 = 0*1.0e-19/transversearea;
        const double hbar = 1.05500000000e-34;
        const double M = 1.409539200000000e-25;
        const double offsetpotential = hbar*hbar*kick*kick/2/M-hbar*omegax/2;
        /* absorbing boundary constants */
        const double dpow = 1;
        const double absorbleft = 4.0e4/pow(2,dpow);
        const double absorbright = 4.0e4/pow(2,dpow);
        const double xleft = -7.0e-5;
        const double widthl = 3.0e-5;
        const double xright = 5.0e-5;
        const double widthr = 3.0e-5;
        
        /* derived constants */
        const double Uohtt = Utt/hbar/transversearea;
        const double Uohtu = Utu/hbar/transversearea;
        const double Uohuu = Uuu/hbar/transversearea;
        const fftw_complex gaa=rcomplex(-losst2,-Uohtt);
        const fftw_complex gbb=rcomplex(-lossu2,-Uohuu);
        const fftw_complex gab=rcomplex(losstu2,Uohtu);      
        const fftw_complex gaas=conj(gaa);
        const fftw_complex gbbs=conj(gbb);
        const fftw_complex gabs=conj(gab);   
        const fftw_complex csqrtgaa=c_sqrt(gaa)*noise;
        const fftw_complex csqrtgbb=c_sqrt(gbb)*noise;
        const fftw_complex csqrtgaas=c_sqrt(gaas)*noise;
        const fftw_complex csqrtgbbs=c_sqrt(gbbs)*noise;
        const double mu = (pow(M, 0.3333333333333333)*pow(1.5*inum*Uohtt*hbar*omegax,0.6666666666666666))/2.;
        const complex csqrtgabo2 = noise*c_sqrt(gab/2);
        const complex csqrtgabso2 = noise*c_sqrt(conj(gabs)/2);
        const double sqrtrhos = noise*sqrt(rhos);
        const fftw_complex Ka =  noise*(gab-2*gaa)/4/dx -losst1/2;
        const fftw_complex Kb =  noise*(gab-2*gbb)/4/dx -lossu1/2;
        const double modgaapgaare = noise*(mod(gaa)+gaa.re)/2/dx;
        const double losst2ondx = -noise*gaa.re/dx;
        const double modgbbpgbbre = noise*(mod(gbb)+gbb.re)/2/dx;
        const double lossu2ondx = -noise*gbb.re/dx;
        const double modgabon4 = noise*mod(gab)/4;
        const double quarter = noise/4;
        const double tworegab = 2*gab.re*noise;
        const double losst2noise = losst2*noise;
        const double lossu2noise = lossu2*noise;
        
        /* numerical shift constants  */
        const double ko=kick;
        const double muo=hbar*omegax/2+offsetpotential;
        const double Eo=muo;
        /*  If muo and Eo are ever different - will need to put time dependence back in functions  */
        const double muomeooh = (muo-Eo)/hbar;
      ]]>
    </globals>  

    <!-- Field to be integrated over -->
    <field>
      <dimensions>x</dimensions>
      <lattice>8192</lattice>
      <domains>(-1.0e-4, 8.0e-5)</domains>
      <samples>1 1 1 1</samples>
      <vector>
        <name>dconstants</name>
        <type>double</type>
        <components>Vt Vu damping ruo2 ruo4dx sqrtruo2</components>
        <fourier_space>no</fourier_space>
        <![CDATA[       
          double u=sqrt(M*omegax/M_PI/hbar)*exp(-x*x*M*omegax/hbar);
          ruo2 = r*u/2;
          ruo4dx = noise*r*u/4/dx;
          sqrtruo2 = noise*sqrt(ruo2);
          damping = x<xleft ? absorbleft*pow(1-cos(M_PI*(xleft-x)/widthl),dpow)
                            : ( x>xright ? absorbright*pow(1-cos(M_PI*(x-xright)/widthr),dpow): 0);   
          Vt = (0.5*M*omegax*omegax*x*x-muo+offsetpotential)/hbar;
          Vu = M*g*x/hbar;
        ]]>
      </vector>
      <vector>
        <name>cconstants</name>
        <type>complex</type>
        <components>kappa kappas</components>
        <fourier_space>no</fourier_space>
        <![CDATA[       
          kappa = kappamax*c_exp(-i*(kick-ko)*x);    
          /*  fftw_complex kappa = kappamax*c_exp(-i*((kick-ko)*x-muomeooh*t));   */
          kappas = conj(kappa);
        ]]>
      </vector>
      <vector>
        <name>main</name>
        <type>complex</type>
        <components>phita phitb phiua phiub theta</components>
        <fourier_space>no</fourier_space>
        <vectors> dconstants cconstants </vectors>
        <![CDATA[
          double realfn = sqrt(inum)*pow(M*omegax/M_PI/hbar,1/4.0)*exp(-x*x*M*omegax/2/hbar);
          phita = rcomplex(realfn,0);
          phitb = rcomplex(realfn,0);
          phiua = rcomplex(0,0);
          phiub = rcomplex(0,0);
          theta = rcomplex(0,0);
        ]]>
      </vector>                 
    </field>

    <!-- The sequence of integrations to perform -->
    <sequence>
      <integrate>
        <algorithm>SIIP</algorithm>
        <interval>1.0e-10</interval>
        <lattice>10</lattice>
        <samples>5 5 5 2</samples>

        <k_operators>
          <constant>yes</constant>
          <operator_names>Lta Ltb Lua Lub</operator_names>
          <![CDATA[
            Lta = rcomplex(Ka.re, Ka.im -hbar/M/2*kx*kx);
            Ltb = rcomplex(Ka.re,-Ka.im +hbar/M/2*kx*kx);
            Lua = rcomplex(Kb.re, Kb.im -hbar/M/2*(kx+ko)*(kx+ko) +Eo/hbar);
            Lub = rcomplex(Kb.re,-Kb.im +hbar/M/2*(kx-ko)*(kx-ko) -Eo/hbar);
          ]]>
        </k_operators>

        <vectors>main dconstants cconstants</vectors>
        <![CDATA[
          complex denst = phitb*phita;
          complex densu = phiub*phiua;
          double moddenst = mod(denst);
          double moddensu = mod(densu);
          complex difft = denst - moddenst;
          complex diffu = densu - moddensu;
          complex csqrtuaub = 2*noise*c_sqrt(phiua)*c_sqrt(phiub);
          complex csqrttatb = 2*noise*c_sqrt(phita)*c_sqrt(phitb);
          
          complex cnoisealpha = csqrtgabo2*c_sqrt(phita)*c_sqrt(phiua);
          complex cnoisebeta = csqrtgabso2*c_sqrt(phitb)*c_sqrt(phiub);
          
          complex invdenst=1/(rhos+denst);
          complex pumpdamp=ruo2*invdenst +ruo4dx*(3*rhos+denst)*invdenst*invdenst*invdenst -damping;
          complex pumpnoise=sqrtruo2*invdenst;
          
          dphita_dt = Lta[phita] + (-i*Vt +pumpdamp +gaa*moddenst -gab*moddensu +csqrtgaa*n_1)*phita
                                 + kappa*phiua -cnoisealpha*(n_5-i*n_7) +pumpnoise*(i*phita*n_9+sqrtrhos*(n_10+i*n_11));

          dphitb_dt = Ltb[phitb] + (i*Vt +pumpdamp +gaas*moddenst -gabs*moddensu +csqrtgaas*n_2)*phitb 
                                 + kappas*phiub -cnoisebeta*(n_6-i*n_8) +pumpnoise*(-i*phitb*n_9+sqrtrhos*(n_10-i*n_11));
                                 
          dphiua_dt = Lua[phiua] + (-i*Vu -damping +gbb*moddensu -gab*moddenst +csqrtgbb*n_3)*phiua
                                 - kappas*phita +cnoisealpha*(n_5+i*n_7);

          dphiub_dt = Lub[phiub] + (i*Vu -damping +gbbs*moddensu -gabs*moddenst +csqrtgbbs*n_4)*phiub
                                 - kappa*phitb +cnoisebeta*(n_6+i*n_8);
  
          dtheta_dt = difft*(csqrtgaa*n_1+csqrtgaas*n_2+difft*losst2noise)+denst*losst2ondx+modgaapgaare*moddenst
                      + diffu*(csqrtgbb*n_3+csqrtgbbs*n_4+diffu*lossu2noise+difft*tworegab)+densu*lossu2ondx+modgbbpgbbre*moddensu
                      + modgabon4*(mod(phitb*phiua)+mod(phita*phiub))*conj(phitb*phiub+phita*phiua)
                      + (difft+diffu)*conj(gab*phitb*phiub+gabs*phita*phiua)*quarter
                      + c_sqrt(phita)*c_sqrt(phiub)*((i*n_5-n_7)*csqrtgabo2*csqrtuaub.im-(i*n_6+n_8)*csqrtgabso2*csqrttatb.im)
                      + c_sqrt(phitb)*c_sqrt(phiua)*((i*n_6-n_8)*csqrtgabso2*csqrtuaub.im-(i*n_5+n_7)*csqrtgabo2*csqrttatb.im);
          /* dtheta_dt=rcomplex(0,0); */
        ]]>
      </integrate>
    </sequence>

    <!-- The output to generate -->
    <output format="ascii" precision="double">
      <filename>gaugeatomlaser.xsil</filename>
      <overwrite>yes</overwrite>
      <group>
        <sampling>
          <fourier_space>no</fourier_space>
          <lattice>0</lattice>
          <moments>integral</moments>
          <![CDATA[
            integral=theta;
          ]]>
        </sampling>
        <post_propagation>
        <fourier_space>no</fourier_space>
        <moments>weight</moments>
        <![CDATA[
          weight = exp(integral.re)*cos(integral.im);
        ]]>
        </post_propagation>
      </group>
      <group>
        <sampling>
          <fourier_space>   no</fourier_space>
          <lattice>         32</lattice>
          <moments>trapped untrapped</moments>
          <![CDATA[
            complex weight = c_exp(_mg0_raw[(_mg0_sample_pointer-1)*_mg0_raw_ncomponents]);
            trapped=phitb*phita*weight;
            untrapped=phiub*phiua*weight;
          ]]>
        </sampling>
      </group>
      <group>
       <sampling>
         <fourier_space>    no</fourier_space>
         <lattice>           0</lattice>
         <moments>ntr nuntr</moments>
         <![CDATA[
           complex weight = c_exp(_mg0_raw[(_mg0_sample_pointer-1)*_mg0_raw_ncomponents]);
           ntr=phitb*phita*weight;
           nuntr=phiub*phiua*weight;
         ]]>
         </sampling>
      </group>
      <group>
        <sampling>
          <fourier_space>    yes</fourier_space>
          <lattice>           128</lattice>
          <moments>ktrapped kuntrapped</moments>
          <![CDATA[
            complex weight = c_exp(_mg0_raw[(_mg0_sample_pointer-1)*_mg0_raw_ncomponents]);
            if(_i0>0)
            {   
            ktrapped = _active_main_main[_main_lattice0*_main_main_ncomponents - _main_main_index_pointer + 1]*phita*weight;
            kuntrapped = _active_main_main[_main_lattice0*_main_main_ncomponents - _main_main_index_pointer + 3]*phiua*weight;
            }
            else
            {
            ktrapped   = phitb*phita*weight;
            kuntrapped = phiub*phiua*weight;
            };
          ]]>
        </sampling>
      </group>
    </output>
</simulation>