path: root/Nlevels_with_MOR.xmds

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

	<name>Nlevels_with_MOR</name>

	<author>Eugeniy Mikhailov, M. Guidry</author>
	<description>
		License GPL.

		Solving split 3-level atom
		with field propagation along spatial axis Z
		with Doppler broadening.


		*                
		*                  -------- |a>  
		*                 /  / \  \
		*           EdL  /  /   \  \ 
		*               /  /     \  \   
		*   |c> ----------/---    \  \ EpR
		*	               / EpL     \  \
		*   |C> --------------      \  \
		*                       EdR  \  \
		*                             \  \
		*                        ------\------ |b>
		*			                        	\
		*                        ------------- |B>
		*
		*


		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[
				const double pi = M_PI; 
				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
				const double Kvec = 2*M_PI/lambda;  // k-vector
				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)]
				
				//  ---------  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 Ga=3.0 *(2*M_PI*1e6);
				//const double G4=3.0 *(2*M_PI*1e6);

				complex g = 10;
				complex gbc = 0.001;
				const complex Split = 0;
				const complex noise = 0;
				
				complex Gab, GAB, Gca, GCA, Gcb, GCB;

				// total decay of i-th level branching ratios. Rij branching of i-th level to j-th
				//const double Rab=0.5, Rac=0.5;


				complex EdLac, EdRac, EpLac, EpRac; 

				complex  rba, rac, rbc;  // density matrix elements

				complex rBA, rAC, rBC; // density matrix for MOR

				// 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="EdLo" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
      <argument name="EdRo" type="real" default_value="2*1.5*(2*M_PI*1e6)" />

      <argument name="EpLo" type="real" default_value="0.05*(2*M_PI*1e6)" /> 
      <argument name="EpRo" type="real" default_value="0.05*(2*M_PI*1e6)" />
 

			<!-- Fields detuning in [1/s] -->
      <argument name="delta_dL"  type="real" default_value="0.0" />
      <argument name="delta_pL"  type="real" default_value="0.0" />
      <argument name="delta_dR"  type="real" default_value="0.0" />
      <argument name="delta_pR"  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;
			]]>
    </arguments>
		<bing />
		<diagnostics /> 
		<fftw plan="estimate" 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="10"   domain="(V_maxwell_min, V_maxwell_max)" />
		</transverse_dimensions>
	</geometry>

	<!-- Rabi frequency --> 
	<vector name="E_field" type="complex" initial_space="t">
		<components>EdL EdR EpL EpR</components>
		<initialisation>
			<![CDATA[
			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
			// as well as time
			EdL = EdLo;
			EdR = EdRo;
			EpL = EpLo;
			//EpR = EpRo;

			EpR=EpRo*exp(-pow( ((t-0.0)/Pwidth),2) );

			//EpL = EpLo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5));	
			//EpR = EpRo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5));

			//E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) );
			]]>
		</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>EdLa EdRa EpLa EpRa</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;

			EdLa=EdL*prob_v_normalized;
			EdRa=EdR*prob_v_normalized;
			EpLa=EpL*prob_v_normalized;
			EpRa=EpR*prob_v_normalized;
			]]>
		</evaluation>
	</computed_vector>

	<!-- Averaged across Maxwell distribution density matrix components -->
	<computed_vector name="density_matrix_averaged" dimensions="t" type="complex">
		<components>rbb_av rBB_av rcc_av rCC_av raa_av rcb_av rab_av rca_av rCB_av rAB_av rCA_av</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;

			rbb_av=rbb*prob_v_normalized;
			rBB_av=rBB*prob_v_normalized;
			rcc_av=rcc*prob_v_normalized;
			rCC_av=rCC*prob_v_normalized;
			raa_av=raa*prob_v_normalized;

			rcb_av=rcb*prob_v_normalized;
			rab_av=rab*prob_v_normalized;
			rca_av=rca*prob_v_normalized;

			rCB_av = rCB*prob_v_normalized;
			rAB_av = rAB*prob_v_normalized;
			rCA_av = rCA*prob_v_normalized;

			]]>
		</evaluation>
	</computed_vector>


	<vector name="density_matrix" type="complex" initial_space="t">
		<components>rbb rBB rcc rCC raa rcb rab rca rCB rAB rCA</components>
		<initialisation>
			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
				<dependencies>E_field_avgd</dependencies>
				<![CDATA[
						EdLac = conj(EdLa);
						EdRac = conj(EdRa);
						EpLac = conj(EpLa);
						EpRac = conj(EpRa);

				rbb = 0.25; rcc = 0.25; rBB = 0.25; rCC = 0.25;
				raa = 0;
				rcb = 0; rab = 0; rca = 0;
				rCB = 0; rAB = 0; rCA = 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 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))-->
          </boundary_condition>
					<![CDATA[

					Gab=g+i*(delta_dL+delta_pR+0*noise);
					GAB=g+i*(delta_pL+delta_dR+0*noise);

					Gca=g-i*(delta_dL);
					GCA=g-i*(delta_pL);

					Gcb=gbc+i*(Split + delta_pR+0*noise);
					GCB=gbc+i*(-Split + delta_pL+0*noise);


                                rba=conj(rab);
                                rac=conj(rca);
                                rbc=conj(rcb);
                                rBA=conj(rAB);
                                rAC=conj(rCA);
                                rBC=conj(rCB);
					draa_dt = -i*EpRac*rab+i*EpRa*rba-i*EdLac*rac+i*EdLa*rca-2*g*raa
					          -i*EdRac*rAB+i*EdRa*rBA-i*EpLac*rAC+i*EpLa*rCA;

					drbb_dt =  i*EpRac*rab-i*EpRa*rba+g*raa-gbc*rbb+gbc*rcc;
					drBB_dt =  i*EdRac*rAB-i*EdRa*rBA+g*raa-gbc*rBB+gbc*rCC;

					drcc_dt =  i*EdLac*rac-i*EdLa*rca+g*raa-gbc*rcc+gbc*rbb;
					drCC_dt =  i*EpLac*rAC-i*EpLa*rCA+g*raa-gbc*rCC+gbc*rBB;

					drab_dt = -Gab*rab+i*EpRa*(rbb-raa)+i*EdLa*rcb;
					drca_dt = -Gca*rca+i*EdLac*(raa-rcc)-i*EpRac*rcb;
					drcb_dt = -Gcb*rcb-i*EpRa*rca+i*EdLac*rab;

					drAB_dt = -Gab*rAB+i*EdRa*(rBB-raa)+i*EpLa*rCB;
					drCA_dt = -Gca*rCA+i*EpLac*(raa-rCC)-i*EdRac*rCB;
					drCB_dt = -GCB*rCB-i*EdRa*rCA+i*EpLac*rAB;


					]]>
        </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[
					dEdL_dz = i*eta*conj(rac) +0*Lt[EdL] ;
					dEdR_dz = i*eta*conj(rBA) +0*Lt[EdR] ;
					dEpL_dz = i*eta*conj(rAC) +0*Lt[EpL] ;
					dEpR_dz = i*eta*conj(rab) +0*Lt[EpR] ;
 
          ]]>
			</operators>
		</integrate>
	</sequence>



	<!-- The output to generate -->
	<output format="binary" filename="Nlevels_with_MOR.xsil">
		<group>
      <sampling basis="t(1000) " initial_sample="yes">
				<dependencies>E_field_avgd</dependencies>
				<moments>IdL_out IpL_out IdR_out IpR_out</moments>
				<![CDATA[
				IdL_out = mod2(EdLa);
				IpL_out = mod2(EpLa);
				IdR_out = mod2(EdRa);
				IpR_out = mod2(EpRa);
				]]>
			</sampling>
		</group>

		<group>
      <sampling basis="t(100) v(10)" initial_sample="yes">
				<dependencies>density_matrix_averaged</dependencies>
				<moments>
					rbb_out rBB_out
					rcc_out rCC_out
					raa_out 
					rcb_re_out rcb_im_out
 					rCB_re_out rCB_im_out
					rab_re_out rab_im_out 
					rAB_re_out rAB_im_out
					rca_re_out rca_im_out 
					rCA_re_out rCA_im_out
				</moments>
				<![CDATA[
				// populations output 
				rbb_out = rbb_av.Re();
				rBB_out = rBB_av.Re();

				rcc_out = rcc_av.Re();	
				rCC_out = rCC_av.Re();
			
				raa_out = raa_av.Re();

				// coherences output 
				rcb_re_out = rcb_av.Re();
				rcb_im_out = rcb_av.Im();

				rCB_re_out = rCB_av.Re();
				rCB_im_out = rCB_av.Im();
				
				rab_re_out = rab_av.Re();
				rab_im_out = rab_av.Im();

				rAB_re_out = rAB_av.Re();
				rAB_im_out = rAB_av.Im();

				rca_re_out = rca_av.Re();
				rca_im_out = rca_av.Im();

				rCA_re_out = rCA_av.Re();
				rCA_im_out = rCA_av.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>
					rbb_out_v rcc_out_v raa_out_v r44_out_v 
					rbc_re_out_v rbc_im_out_v rba_re_out_v rba_im_out_v r14_re_out_v r14_im_out_v
					                      rca_re_out_v rca_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 
				rbb_out_v = rbb.Re();
				rcc_out_v = rcc.Re();
				raa_out_v = raa.Re();
				r44_out_v = r44.Re();
				// coherences output 
				rbc_re_out_v = rbc.Re();
				rbc_im_out_v = rbc.Im();
				rba_re_out_v = rba.Re();
				rba_im_out_v = rba.Im();
				r14_re_out_v = r14.Re();
				r14_im_out_v = r14.Im();
				rca_re_out_v = rca.Re();
				rca_im_out_v = rca.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: 
-->