full simulation of 5 level MOR to separate folder

1 files changed, 0 insertions, 481 deletions

diff --git a/Nlevels_with_MOR.xmds b/Nlevels_with_MOR.xmds
deleted file mode 100644
index 99d88a0..0000000
--- a/Nlevels_with_MOR.xmds
+++ /dev/null
@@ -1,481 +0,0 @@
-<?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 = Gamma_super;
-				complex gbc = (2*M_PI)*1e3; // 1/s untits
-				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; 
-
-				// 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="2"   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*exp(-pow( ((t-0.0)/Pwidth),2) );
-			EdR=EdRo*exp(-pow( ((t-0.0)/Pwidth),2) );
-			EpL=EpLo*exp(-pow( ((t-0.0)/Pwidth),2) );
-			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);
-
-					complex rba=conj(rab);
-					complex rac=conj(rca);
-					complex rbc=conj(rcb);
-					complex rBA=conj(rAB);
-					complex rAC=conj(rCA);
-					complex 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*(rca) +0*Lt[EdL] ;
-					dEdR_dz = i*eta*(rAB) +0*Lt[EdR] ;
-					dEpL_dz = i*eta*(rCA) +0*Lt[EpL] ;
-					dEpR_dz = i*eta*(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(2)" 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: 
--->