path: root/xmds2/realistic_Rb/realistic_Rb.xmds

<?xml version="1.0"?>
<simulation xmds-version="2">
	<testing>
		<arguments> --Ndens=1e15 --Lcell=10.0e-2  --Temperature=1e-3  --Pwidth=0.4e-6  --delta1=0 --delta2=0 --delta3=0  --E1o=0 --E2o=1e2 --E3o=0 --E4o=0</arguments>
		<xsil_file name="realistic_Rb.xsil" expected="realistic_Rb_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5">
			<moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" />
		</xsil_file>
	</testing>

	<name>realistic_Rb</name>

	<author>Eugeniy Mikhailov</author>
	<description>
		License GPL.

		Solving simplified Rb atom model
		with fields propagation along spatial axis Z
		with Doppler broadening.


		We assume four-wave mixing condition when w3-w4=w2-w1 i.e. fields E3 and E4 drive the same
		resonance as fields E2 and E1.


		*              --------------- | F=1, 2P_3/2 >
		*                \        \    
		*                 \ E3_r   \    -------- | F=2, 2P_+1/2 >
		*                  \   E4_r \   /    \
		*                   \        \ / E2_l \ 
		*                    \        /        \ E1_l
		*   | F=2, 2S_1/2 > --------------      \
		*                               \        \
		*                                \        \
		*                               ------------- | F=1, 2S_1/2 >
		*            


		We are solving 
			dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher than i.
			Note that E is actually a Rabi frequency of electromagnetic field not the EM field
		in xmds terms it looks like
			dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space

		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
		normalize it to something else look drho/dt equation.
		No need to renormalizes eta as long as its express through
		the upper level decay rate in the same units as Rabi frequency.
	</description>

	<features>
		<globals>
			<![CDATA[
				// Some numerical constants
				const double pi = M_PI; 

				// atom related constants
				//read from Mathematica generated Constants.txt
//---------------- Constants.txt starts ------------------
const double ha0 = 2.1471788680034824e10;
const double ha1 = 2.558764384495815e9;
const double g1 = 3.612847284945266e7;
const double ha2 = 5.323020344462938e8;
const double hb2 = 7.85178251911697e7;
const double g2 = 3.8117309832741246e7;
const double lambda1 = 7.949788511562656e-7;
const double lambda2 = 7.802412096860509e-7;
const double eta1 = 5.450949336831401e-6;
const double eta2 = 5.5397657647874e-6;
const double rt6 = 2.449489742783178;
const double rt3 = 1.7320508075688772;
const double rt2 = 1.4142135623730951;
//---------------- Constants.txt ends   ------------------

				const double c=3.e8;
				const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
				const double lambda=794.7e-9; //wavelength in m
				// Fields k-vector
				const double Kvec = 2*M_PI/lambda;  
				// Simplified k-vectors
				const double Kvec1 = Kvec, Kvec2=Kvec, Kvec3=Kvec;

				//  ---------  Atom and cell properties -------------------------
				// range of Maxwell distribution atomic velocities
				const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] 
				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)

				// Average sqrt(v^2) in Maxwell distribution for one dimension
				// Maxwell related parameters will be calculated in <arguments> section
				double v_thermal_averaged=0;
				// Maxwell distribution velocities range to take in account in [m/s]
				double V_maxwell_min = 0, V_maxwell_max = 0;

				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
				const double gt=0.01 *(2*M_PI*1e6);

				// Larmor frequency 
				double WL=0;


				// inner use variables 
				double probability_v; // will be used as p(v) in Maxwell distribution

			]]>
		</globals>
		<validation kind="run-time"/> <!--allows to put ranges as variables-->
		<benchmark />
    <arguments>
			<!-- Rabi frequency divided by 2 in [1/s] -->
      <argument name="E1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
      <argument name="E2o" type="real" default_value="0.05*(2*M_PI*1e6)" />
      <argument name="E3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" />
      <argument name="E4o" type="real" default_value=".01*(2*M_PI*1e6)" />
			<!-- Fields detuning in [1/s] -->
      <argument name="delta1"  type="real" default_value="0.0" />
      <argument name="delta2"  type="real" default_value="0.0" />
      <argument name="delta3"  type="real" default_value="0.0" />
			<!--Pulse duration/width [s] -->
			<argument name="Pwidth"  type="real" default_value="0.1e-6" />
			<!--  Atom and cell properties -->
			<!--Cell length [m] -->
			<argument name="Lcell"  type="real" default_value="1.5e-2" />
			<!--Density of atoms [1/m^3] -->
			<argument name="Ndens"  type="real" default_value="1e15" />
			<!--Atoms temperature [K] -->
			<!--TODO: looks like Temperature > 10 K knocks solver, 
					 I am guessing detunings are too large and thus it became a stiff equation-->
			<!--! make sure it is not equal to zero!-->
			<argument name="Temperature"  type="real" default_value="5" />
			<!-- This will be executed after arguments/parameters are parsed -->
			<!-- Read the code Luke: took me a while of reading the xmds2 sources to find it -->
			<![CDATA[
				// Average sqrt(v^2) in Maxwell distribution for one dimension
				if (Temperature == 0)
					_LOG(_ERROR_LOG_LEVEL, "ERROR: Temperature should be >0 to provide range for Maxwell velocity distribution\n");
				v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
				// Maxwell distribution velocities range to take in account in [m/s]
				// there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
				V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
			]]>
    </arguments>
		<bing />
		<diagnostics /> 
		<fftw plan="estimate" threads="1" />
		<!--<fftw plan="patient" threads="1" />-->
		<!-- I don't see any speed up on 6 core CPU even if use threads="6" -->
		<!--<openmp />-->
		<auto_vectorise />
		<halt_non_finite />
	</features>

	<!-- 'z', 't', and 'v'  to have dimensions [m], [s], and [m/s]   -->
	<geometry>
		<propagation_dimension> z </propagation_dimension>
		<transverse_dimensions>
			<!-- IMPORTANT: looks like having a lot of points in time helps with convergence.
					 I suspect that time spacing should be small enough to catch
					 all pulse harmonics and more importantly 1/dt should be larger than
					 the largest detuning (including Doppler shifts).
					 Unfortunately calculation time is proportional to lattice size
					 so we cannot just blindly increase it.
					 Some rules of thumb:  
						* lattice="1000"   domain="(-1e-6, 1e-6)" 
							was good enough detunings up to 155 MHz (980 rad/s) notice that 1/dt=500 MHz
						* lattice="10000"   domain="(-1e-6, 1e-6)" 
							works for Doppler averaging in up to 400K for Rb when lasers are zero detuned
			 -->
			<dimension name="t"   lattice="10000"   domain="(-1e-6, 1e-6)" />
			<dimension name="v"   lattice="2"   domain="(V_maxwell_min, V_maxwell_max)" />
		</transverse_dimensions>
	</geometry>

	<!-- Rabi frequency --> 
	<vector name="E_field" type="complex" initial_space="t">
		<components>E1 E2 E3 E4</components>
		<initialisation>
			<![CDATA[
			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
			// as well as time
			E1=E1o;
			E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) );
			E3=E3o;
			E4=E4o;
			]]>
		</initialisation>
	</vector>

	<!--Maxwell distribution probability p(v)-->
	<vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
		<components>probability_v</components>
		<initialisation>
			<![CDATA[
			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
			probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); 
			]]>
		</initialisation>
	</vector>

	<!--Maxwell distribution norm sum(p(v))
			 Needed since we sum over the grid instead of true integral,
			 we also have finite cut off velocities-->
	<vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real">
		<components>probability_v_norm</components>
		<initialisation>
			<dependencies basis="v">Maxwell_distribution_probabilities</dependencies>
			<![CDATA[
			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
			probability_v_norm=probability_v;
			]]>
		</initialisation>
	</vector>


	<!-- Averaged across Maxwell distribution fields amplitudes -->
	<computed_vector name="E_field_avgd" dimensions="t" type="complex">
		<components>E1a E2a E3a E4a</components>
		<evaluation>
			<dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
			<![CDATA[
			double prob_v_normalized=probability_v/probability_v_norm;
			E1a=E1*prob_v_normalized;
			E2a=E2*prob_v_normalized;
			E3a=E3*prob_v_normalized;
			E4a=E4*prob_v_normalized;
			]]>
		</evaluation>
	</computed_vector>

	<vector name="density_matrix" type="complex" initial_space="t">
		<components>
			<!-- read from Mathematica generated RbChosenRho.txt -->
<!-- ############### RbChosenRho.txt starts ################ -->
r0101 r0113 r0202 r0214 r0303 r0309 r0315 r0404 r0410 r0416 r0505 r0511 r0602 r0606 r0614 r0703 r0707 r0709 r0715 r0804 r0808 r0810 r0816 r0909 r0915 r1010 r1016 r1111 r1313 r1414 r1515 r1616 
<!-- ############### RbChosenRho.txt ends  ################# -->
		</components>
		<initialisation>
			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
			<![CDATA[
			// Note: 
			// convergence is really slow if all populations concentrated at the bottom level |1>
			// this is because if r11=1, everything else is 0 and then every small increment 
			// seems to be huge and adaptive solver makes smaller and smaller steps.
			// As quick and dirty fix I reshuffle initial population  
			// so some of the population sits at the  second ground level |2>
			// TODO: Fix above. Make the equation of motion for r11 
			//       and express other level, let's say r44
			//       through population normalization

			//read from Mathematica generated RbInits.txt
//---------------- RbInits.txt starts ------------------
r0101 = 0.125;
r0113 = 0;
r0202 = 0.125;
r0214 = 0;
r0303 = 0.125;
r0309 = 0;
r0315 = 0;
r0404 = 0.125;
r0410 = 0;
r0416 = 0;
r0505 = 0.125;
r0511 = 0;
r0602 = 0;
r0606 = 0.125;
r0614 = 0;
r0703 = 0;
r0707 = 0.125;
r0709 = 0;
r0715 = 0;
r0804 = 0;
r0808 = 0.125;
r0810 = 0;
r0816 = 0;
r0909 = 0;
r0915 = 0;
r1010 = 0;
r1016 = 0;
r1111 = 0;
r1313 = 0;
r1414 = 0;
r1515 = 0;
r1616 = 0;
//---------------- RbInits.txt ends  ------------------
			]]>
		</initialisation>
	</vector>

	<sequence>
		<!--For this set of conditions ARK45 is faster than ARK89-->
		<!--ARK45 is good for small detuning when all frequency like term are close to zero-->
		<integrate algorithm="ARK45" tolerance="1e-5" interval="Lcell"> 
		<!--<integrate algorithm="SI" steps="200" interval="Lcell"> -->
		<!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small-->
		<!--When ARK45 works it about 3 times faster than RK4 with 1000 steps-->
		<!--<integrate algorithm="RK4" steps="100" interval="1.5e-2">-->
		<!--SIC algorithm seems to be much slower and needs fine 'z'  step tuning and much finer time grid-->
		<!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04-->

		<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->

			<samples>100</samples>
			<!--<samples>100 100</samples>-->
			<!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4-->
			<!--<samples>100 100 100 100</samples>--> 
			<operators>
        <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
					<integration_vectors>density_matrix</integration_vectors>
          <dependencies>E_field_avgd</dependencies>
          <boundary_condition kind="left">
						<!--This set boundary condition at all 'z'  and left border of 't' (i.e. min(t))-->
						<!--
            <![CDATA[
							r11 = 0; r22 = 1; r33 = 0; r44 = 0;
							r12 = 0; r13 = 0; r14 = 0;
							r23 = 0; r24 = 0;
							r34 = 0;
							printf("z= %g, t= %g\n", z, t);
            ]]>
						-->
          </boundary_condition>
					<![CDATA[
						// Equations of motions according to Simon's mathematica code
						//read from Mathematica generated RbEquations.txt
//---------------- RbEquations.txt starts ------------------
dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (conj(r0113)*conj(E4a))/(4.*rt6));
dr0113_dt = (-(gt*r0113) - (gt + g2)*r0113)/2. - i*(WL*r0113 - ((2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0113 + (r0101*conj(E4a))/(4.*rt6) - (r1313*conj(E4a))/(4.*rt6));
dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (conj(r0214)*conj(E4a))/8.);
dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (conj(r0602)*conj(E3a))/(8.*rt3) + (r0202*conj(E4a))/8. - (r1414*conj(E4a))/8.);
dr0303_dt = gt/8. - gt*r0303 + (g1*r0909)/12. + (g1*r1010)/3. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1515)/4. - i*((r0309*E1a)/(4.*rt6) + (r0315*E4a)/8. - (conj(r0309)*conj(E1a))/(4.*rt6) - (conj(r0315)*conj(E4a))/8.);
dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*conj(E1a))/(4.*rt6) - (r0909*conj(E1a))/(4.*rt6) - (conj(r0703)*conj(E2a))/(4.*rt6) - (conj(r0915)*conj(E4a))/8.);
dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*conj(E1a))/(4.*rt6) - (conj(r0703)*conj(E3a))/8. + (r0303*conj(E4a))/8. - (r1515*conj(E4a))/8.);
dr0404_dt = gt/8. - gt*r0404 + (g1*r1010)/4. + (g1*r1111)/4. + (g2*r1414)/4. + (g2*r1515)/12. + (g2*r1616)/6. - i*((r0410*E1a)/(4.*rt2) + (r0416*E4a)/(4.*rt6) - (conj(r0410)*conj(E1a))/(4.*rt2) - (conj(r0416)*conj(E4a))/(4.*rt6));
dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*conj(E1a))/(4.*rt2) - (r1010*conj(E1a))/(4.*rt2) - (conj(r0804)*conj(E2a))/(4.*rt6) - (conj(r1016)*conj(E4a))/(4.*rt6));
dr0416_dt = (-(gt*r0416) - (gt + g2)*r0416)/2. - i*(-(WL*r0416)/2. - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0416 - (r1016*conj(E1a))/(4.*rt2) - (conj(r0804)*conj(E3a))/(4.*rt2) + (r0404*conj(E4a))/(4.*rt6) - (r1616*conj(E4a))/(4.*rt6));
dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (conj(r0511)*conj(E1a))/4.);
dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*conj(E1a))/4. - (r1111*conj(E1a))/4.);
dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (conj(r0214)*conj(E3a))/(8.*rt3));
dr0606_dt = gt/8. - gt*r0606 + (g1*r0909)/12. + (g1*r1010)/12. + (g2*r1313)/4. + (g2*r1414)/12. - i*(-(r0614*E3a)/(8.*rt3) + (conj(r0614)*conj(E3a))/(8.*rt3));
dr0614_dt = (-(gt*r0614) - (gt + g2)*r0614)/2. - i*((-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0614 - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0614 - (r0606*conj(E3a))/(8.*rt3) + (r1414*conj(E3a))/(8.*rt3) + (r0602*conj(E4a))/8.);
dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (conj(r0309)*conj(E2a))/(4.*rt6) + (conj(r0315)*conj(E3a))/8.);
dr0707_dt = gt/8. - gt*r0707 + (g1*r0909)/12. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1414)/3. + (g2*r1515)/4. - i*(-(r0709*E2a)/(4.*rt6) - (r0715*E3a)/8. + (conj(r0709)*conj(E2a))/(4.*rt6) + (conj(r0715)*conj(E3a))/8.);
dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*conj(E1a))/(4.*rt6) - (r0707*conj(E2a))/(4.*rt6) + (r0909*conj(E2a))/(4.*rt6) + (conj(r0915)*conj(E3a))/8.);
dr0715_dt = (-(gt*r0715) - (gt + g2)*r0715)/2. - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0715 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0715 + (r0915*conj(E2a))/(4.*rt6) - (r0707*conj(E3a))/8. + (r1515*conj(E3a))/8. + (r0703*conj(E4a))/8.);
dr0804_dt = -(gt*r0804) - i*((WL*r0804)/2. + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0804 + (r0810*E1a)/(4.*rt2) + (r0816*E4a)/(4.*rt6) + (conj(r0410)*conj(E2a))/(4.*rt6) + (conj(r0416)*conj(E3a))/(4.*rt2));
dr0808_dt = gt/8. - gt*r0808 + (g1*r1010)/12. + (g1*r1111)/12. + (g2*r1414)/12. + (g2*r1515)/4. + (g2*r1616)/2. - i*(-(r0810*E2a)/(4.*rt6) - (r0816*E3a)/(4.*rt2) + (conj(r0810)*conj(E2a))/(4.*rt6) + (conj(r0816)*conj(E3a))/(4.*rt2));
dr0810_dt = (-(gt*r0810) - (gt + g1)*r0810)/2. - i*((delta1 + v*Kvec1)*r0810 + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0810 + (r0804*conj(E1a))/(4.*rt2) - (r0808*conj(E2a))/(4.*rt6) + (r1010*conj(E2a))/(4.*rt6) + (conj(r1016)*conj(E3a))/(4.*rt2));
dr0816_dt = (-(gt*r0816) - (gt + g2)*r0816)/2. - i*((WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0816 - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0816 + (r1016*conj(E2a))/(4.*rt6) - (r0808*conj(E3a))/(4.*rt2) + (r1616*conj(E3a))/(4.*rt2) + (r0804*conj(E4a))/(4.*rt6));
dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (conj(r0309)*conj(E1a))/(4.*rt6) - (conj(r0709)*conj(E2a))/(4.*rt6));
dr0915_dt = (-((gt + g1)*r0915) - (gt + g2)*r0915)/2. - i*((-WL/6. - delta1 - v*Kvec1)*r0915 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0915 - (r0315*E1a)/(4.*rt6) + (r0715*E2a)/(4.*rt6) - (conj(r0709)*conj(E3a))/8. + (conj(r0309)*conj(E4a))/8.);
dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (conj(r0410)*conj(E1a))/(4.*rt2) - (conj(r0810)*conj(E2a))/(4.*rt6));
dr1016_dt = (-((gt + g1)*r1016) - (gt + g2)*r1016)/2. - i*(-((delta1 + v*Kvec1)*r1016) - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1016 - (r0416*E1a)/(4.*rt2) + (r0816*E2a)/(4.*rt6) - (conj(r0810)*conj(E3a))/(4.*rt2) + (conj(r0410)*conj(E4a))/(4.*rt6));
dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (conj(r0511)*conj(E1a))/4.);
dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (conj(r0113)*conj(E4a))/(4.*rt6));
dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (conj(r0614)*conj(E3a))/(8.*rt3) + (conj(r0214)*conj(E4a))/8.);
dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (conj(r0715)*conj(E3a))/8. + (conj(r0315)*conj(E4a))/8.);
dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (conj(r0816)*conj(E3a))/(4.*rt2) + (conj(r0416)*conj(E4a))/(4.*rt6));
//---------------- RbEquations.txt ends  ------------------
					]]>
        </operator>
				<!--
							 According to xmds2 docs operator kind="ip" should be faster
							 but our codes runs about 5% to 10% slower with it.
							 Maybe because we very close to the stiff condition so I use "ex" kind
							 <operator kind="ip" constant="yes">
					 -->
				<operator kind="ex" constant="yes" type="imaginary">
					<operator_names>Lt</operator_names>
					<![CDATA[
					Lt = -i/c*kt;
					]]>
				</operator>
        <integration_vectors>E_field</integration_vectors>
				<dependencies>density_matrix</dependencies>
				<![CDATA[
				//read from Mathematica generated RbPropEquations.txt
//---------------- RbPropEquations.txt starts ------------------
dE1_dz = 0.16666666666666666*i*Ndens*(2.449489742783178*conj(r0309) + 4.242640687119286*conj(r0410) + 6.*conj(r0511))*eta1 - Lt[E1];
dE2_dz = -0.4082482904638631*i*Ndens*(conj(r0709) + conj(r0810))*eta1 - Lt[E2];
dE3_dz = -0.3333333333333333*i*Ndens*(1.7320508075688772*conj(r0614) + 3.*conj(r0715) + 4.242640687119286*conj(r0816))*eta2 - Lt[E3];
dE4_dz = (i*Ndens*(2.449489742783178*conj(r0113) + 3*conj(r0214) + 3*conj(r0315) + 2.449489742783178*conj(r0416))*eta2)/3. - Lt[E4];
//---------------- RbPropEquations.txt ends  ------------------
				]]>
			</operators>
		</integrate>
	</sequence>



	<!-- The output to generate -->
	<output format="binary" filename="realistic_Rb.xsil">
		<group>
      <sampling basis="t(1000) " initial_sample="yes">
				<dependencies>E_field_avgd</dependencies>
				<moments>I1_out I2_out I3_out I4_out</moments>
				<![CDATA[
				I1_out = mod2(E1a);
				I2_out = mod2(E2a);
				I3_out = mod2(E3a);
				I4_out = mod2(E4a);
				]]>
			</sampling>
		</group>

		<!-- use the following two groups only for debuging 
				 otherwise they are quite useless and have to much information 
				 in 3D space (z,t,v) -->
		<!--
		<group>
      <sampling basis="t(100) v(10)" initial_sample="yes">
				<dependencies>E_field</dependencies>
				<moments>I1_out_v I2_out_v I3_out_v I4_out_v</moments>
				<![CDATA[
				// light field intensity distribution in velocity subgroups 
				I1_out_v = mod2(E1);
				I2_out_v = mod2(E2);
				I3_out_v = mod2(E3);
				I4_out_v = mod2(E4);
				]]>
			</sampling>
		</group>

		<group>
      <sampling basis="t(100) v(10)" initial_sample="yes">
				<dependencies>density_matrix</dependencies>
				<moments>
					r11_out_v r22_out_v r33_out_v r44_out_v 
					r12_re_out_v r12_im_out_v r13_re_out_v r13_im_out_v r14_re_out_v r14_im_out_v
					                      r23_re_out_v r23_im_out_v r24_re_out_v r24_im_out_v 
					                                            r34_re_out_v r34_im_out_v
				</moments>
				<![CDATA[
				// density matrix distribution in velocity subgroups 
				// populations output 
				r11_out_v = r11.Re();
				r22_out_v = r22.Re();
				r33_out_v = r33.Re();
				r44_out_v = r44.Re();
				// coherences output 
				r12_re_out_v = r12.Re();
				r12_im_out_v = r12.Im();
				r13_re_out_v = r13.Re();
				r13_im_out_v = r13.Im();
				r14_re_out_v = r14.Re();
				r14_im_out_v = r14.Im();
				r23_re_out_v = r23.Re();
				r23_im_out_v = r23.Im();
				r24_re_out_v = r24.Re();
				r24_im_out_v = r24.Im();
				r34_re_out_v = r34.Re();
				r34_im_out_v = r34.Im();
				]]>
			</sampling>
		</group>
		-->

	</output>

</simulation>
	
<!-- 
vim: ts=2 sw=2 foldmethod=indent: 
-->