path: root/xmds2/realistic_Rb/realistic_Rb.xmds

<?xml version="1.0"?>
<simulation xmds-version="2">

	<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 then 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; 
				// proportional to splitting ratios  sqrt(6) , sqrt(3), sqrt(2)
				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);

				// Natural linewidth  of j's level in [1/s]
				const double g1 = 3.612847284945266e7;
				const double g2 = 3.8117309832741246e7;

				// levels energy
				const double ha0 = 2.1471788680034824e10;
				const double ha1 = 2.558764384495815e9;
				const double ha2 = 5.323020344462938e8;
				const double hb2 = 7.85178251911697e7;

				// Larmor frequency 
				double WL=0;



				complex  E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies

				// density matrix elements which calculated via Hermitian property r_ij=conj(r_ji)
				complex  
						r1301, 
						r1402, 
						r0903, 
						r1503, 
						r1004, 
						r1604, 
						r1105, 
						r0206, 
						r1406, 
						r0307, 
						r0907, 
						r1507, 
						r0408, 
						r1008, 
						r1608, 
						r1509, 
						r1610; 


				// 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="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="100"   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>

	<!-- Averaged across Maxwell distribution density matrix components -->
	<computed_vector name="density_matrix_averaged" dimensions="t" type="complex">
		<components>
				r0101a
				r0113a
				r0202a
				r0214a
				r0303a
				r0309a
				r0315a
				r0404a
				r0410a
				r0416a
				r0505a
				r0511a
				r0602a
				r0606a
				r0614a
				r0703a
				r0707a
				r0709a
				r0715a
				r0804a
				r0808a
				r0810a
				r0816a
				r0909a
				r0915a
				r1010a
				r1016a
				r1111a
				r1313a
				r1414a
				r1515a
				r1616a
		</components>
		<evaluation>
			<dependencies basis="v">density_matrix Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
			<![CDATA[
			double prob_v_normalized=probability_v/probability_v_norm;

				r0101a = r0101*prob_v_normalized;
				r0113a = r0113*prob_v_normalized;
				r0202a = r0202*prob_v_normalized;
				r0214a = r0214*prob_v_normalized;
				r0303a = r0303*prob_v_normalized;
				r0309a = r0309*prob_v_normalized;
				r0315a = r0315*prob_v_normalized;
				r0404a = r0404*prob_v_normalized;
				r0410a = r0410*prob_v_normalized;
				r0416a = r0416*prob_v_normalized;
				r0505a = r0505*prob_v_normalized;
				r0511a = r0511*prob_v_normalized;
				r0602a = r0602*prob_v_normalized;
				r0606a = r0606*prob_v_normalized;
				r0614a = r0614*prob_v_normalized;
				r0703a = r0703*prob_v_normalized;
				r0707a = r0707*prob_v_normalized;
				r0709a = r0709*prob_v_normalized;
				r0715a = r0715*prob_v_normalized;
				r0804a = r0804*prob_v_normalized;
				r0808a = r0808*prob_v_normalized;
				r0810a = r0810*prob_v_normalized;
				r0816a = r0816*prob_v_normalized;
				r0909a = r0909*prob_v_normalized;
				r0915a = r0915*prob_v_normalized;
				r1010a = r1010*prob_v_normalized;
				r1016a = r1016*prob_v_normalized;
				r1111a = r1111*prob_v_normalized;
				r1313a = r1313*prob_v_normalized;
				r1414a = r1414*prob_v_normalized;
				r1515a = r1515*prob_v_normalized;
				r1616a = r1616*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
				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 then 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[
						E1ac = conj(E1a);
						E2ac = conj(E2a);
						E3ac = conj(E3a);
						E4ac = conj(E4a);

						// Density matrix is Hermitian so we use r_ij=conj(r_ji)

						r1301 = conj(r0113);
						r1402 = conj(r0214);
						r0903 = conj(r0309);
						r1503 = conj(r0315);
						r1004 = conj(r0410);
						r1604 = conj(r0416);
						r1105 = conj(r0511);
						r0206 = conj(r0602);
						r1406 = conj(r0614);
						r0307 = conj(r0703);
						r0907 = conj(r0709);
						r1507 = conj(r0715);
						r0408 = conj(r0804);
						r1008 = conj(r0810);
						r1608 = conj(r0816);
						r1509 = conj(r0915);
						r1610 = conj(r1016);

						// Equations of motions according to Simon's mathematica code
						dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (r1301*E4ac)/(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*E4ac)/(4.*rt6) - (r1313*E4ac)/(4.*rt6));
						dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (r1402*E4ac)/8.);
						dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (r0206*E3ac)/(8.*rt3) + (r0202*E4ac)/8. - (r1414*E4ac)/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. - (r0903*E1ac)/(4.*rt6) - (r1503*E4ac)/8.);
						dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*E1ac)/(4.*rt6) - (r0909*E1ac)/(4.*rt6) - (r0307*E2ac)/(4.*rt6) - (r1509*E4ac)/8.);
						dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*E1ac)/(4.*rt6) - (r0307*E3ac)/8. + (r0303*E4ac)/8. - (r1515*E4ac)/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) - (r1004*E1ac)/(4.*rt2) - (r1604*E4ac)/(4.*rt6));
						dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*E1ac)/(4.*rt2) - (r1010*E1ac)/(4.*rt2) - (r0408*E2ac)/(4.*rt6) - (r1610*E4ac)/(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*E1ac)/(4.*rt2) - (r0408*E3ac)/(4.*rt2) + (r0404*E4ac)/(4.*rt6) - (r1616*E4ac)/(4.*rt6));
						dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (r1105*E1ac)/4.);
						dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*E1ac)/4. - (r1111*E1ac)/4.);
						dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (r1402*E3ac)/(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) + (r1406*E3ac)/(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*E3ac)/(8.*rt3) + (r1414*E3ac)/(8.*rt3) + (r0602*E4ac)/8.);
						dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (r0903*E2ac)/(4.*rt6) + (r1503*E3ac)/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. + (r0907*E2ac)/(4.*rt6) + (r1507*E3ac)/8.);
						dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*E1ac)/(4.*rt6) - (r0707*E2ac)/(4.*rt6) + (r0909*E2ac)/(4.*rt6) + (r1509*E3ac)/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*E2ac)/(4.*rt6) - (r0707*E3ac)/8. + (r1515*E3ac)/8. + (r0703*E4ac)/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) + (r1004*E2ac)/(4.*rt6) + (r1604*E3ac)/(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) + (r1008*E2ac)/(4.*rt6) + (r1608*E3ac)/(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*E1ac)/(4.*rt2) - (r0808*E2ac)/(4.*rt6) + (r1010*E2ac)/(4.*rt6) + (r1610*E3ac)/(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*E2ac)/(4.*rt6) - (r0808*E3ac)/(4.*rt2) + (r1616*E3ac)/(4.*rt2) + (r0804*E4ac)/(4.*rt6));
						dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (r0903*E1ac)/(4.*rt6) - (r0907*E2ac)/(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) - (r0907*E3ac)/8. + (r0903*E4ac)/8.);
						dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (r1004*E1ac)/(4.*rt2) - (r1008*E2ac)/(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) - (r1008*E3ac)/(4.*rt2) + (r1004*E4ac)/(4.*rt6));
						dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (r1105*E1ac)/4.);
						dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (r1301*E4ac)/(4.*rt6));
						dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (r1406*E3ac)/(8.*rt3) + (r1402*E4ac)/8.);
						dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (r1507*E3ac)/8. + (r1503*E4ac)/8.);
						dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (r1608*E3ac)/(4.*rt2) + (r1604*E4ac)/(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[
					dE1_dz = i*0.16666666666666666*eta1*(2.449489742783178*r0309 + 4.242640687119286*r0410 + 6.*r0511) - Lt[E1];
					dE2_dz = -i*0.8164965809277261*eta1*(r0709 + r0810) - Lt[E2];
					dE3_dz = -i*1.*eta2*(1.7320508075688772*r0614 + 3.*r0715 + 4.242640687119286*r0816) - Lt[E3];
					dE4_dz = i*(4*eta2*(2.449489742783178*r0113 + 3*r0214 + 3*r0315 + 2.449489742783178*r0416))/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>

		<!--
		<group>
      <sampling basis="t(100) v(10)" initial_sample="yes">
				<dependencies>density_matrix_averaged</dependencies>
				<moments>
					r11_out r22_out r33_out r44_out 
					r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out
					                      r23_re_out r23_im_out r24_re_out r24_im_out 
					                                            r34_re_out r34_im_out
				</moments>
				<![CDATA[
				// populations output 
				r11_out = r11a.Re();
				r22_out = r22a.Re();
				r33_out = r33a.Re();
				r44_out = r44a.Re();
				// coherences output 
				r12_re_out = r12a.Re();
				r12_im_out = r12a.Im();
				r13_re_out = r13a.Re();
				r13_im_out = r13a.Im();
				r14_re_out = r14a.Re();
				r14_im_out = r14a.Im();
				r23_re_out = r23a.Re();
				r23_im_out = r23a.Im();
				r24_re_out = r24a.Re();
				r24_im_out = r24a.Im();
				r34_re_out = r34a.Re();
				r34_im_out = r34a.Im();
				]]>
			</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: 
-->