summaryrefslogtreecommitdiff
path: root/Nlevels_with_MOR.xmds
diff options
context:
space:
mode:
Diffstat (limited to 'Nlevels_with_MOR.xmds')
-rw-r--r--Nlevels_with_MOR.xmds487
1 files changed, 487 insertions, 0 deletions
diff --git a/Nlevels_with_MOR.xmds b/Nlevels_with_MOR.xmds
new file mode 100644
index 0000000..16b135f
--- /dev/null
+++ b/Nlevels_with_MOR.xmds
@@ -0,0 +1,487 @@
+<?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:
+-->
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 09:32:34 +0000