|
gaugeatomlaser.xmds
Script source: gaugeatomlaser.xmds.gz
<?xmds version="1.0"?>
<!--Stochastic Atom laser simulation-->
<!-- $Id: gaugeatomlaser_body.part 999 2004-08-03 05:42:47Z cochrane $ -->
<!-- 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">
<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>
Generated by GNU enscript 1.6.3.
|