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
				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 = 0.00007949788511562656;
				const double lambda2 = 0.00007802412096860508;
				const double rt6 = 2.449489742783178;
				const double rt3 = 1.7320508075688772;
				const double rt2 = 1.4142135623730951;

				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;

				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
				// eta will be calculated in the <arguments> section
				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
				double eta1=0, eta2=0, eta3=0;
				
				//  ---------  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; 

				// eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
				eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;
				// !FIXME over simplification: we should use relevant levels linewidths
				eta1 = eta;
				eta2 = eta;
				eta3 = eta;
			]]>
    </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)-->
	<computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
		<components>probability_v</components>
		<evaluation>
			<![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 ); 
			]]>
		</evaluation>
	</computed_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-->
	<computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real">
		<components>probability_v_norm</components>
		<evaluation>
			<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;
			]]>
		</evaluation>
	</computed_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>
				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
		</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
				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;
			]]>
		</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
						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));
					]]>
        </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
				dE1_dz = 0.16666666666666666*i*(2.449489742783178*conj(r0309) + 4.242640687119286*conj(r0410) + 6.*conj(r0511))*eta1 - Lt[E1];
				dE2_dz = -0.4082482904638631*i*(conj(r0709) + conj(r0810))*eta1 - Lt[E2];
				dE3_dz = -0.3333333333333333*i*(1.7320508075688772*conj(r0614) + 3.*conj(r0715) + 4.242640687119286*conj(r0816))*eta2 - Lt[E3];
				dE4_dz = (i*(2.449489742783178*conj(r0113) + 3*conj(r0214) + 3*conj(r0315) + 2.449489742783178*conj(r0416))*eta2)/3. - Lt[E4];
				]]>
			</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: 
-->