diff options
| -rw-r--r-- | xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds | 208 |
1 files changed, 137 insertions, 71 deletions
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds index ae2fe75..7d33920 100644 --- a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds @@ -45,11 +45,23 @@ <![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 N=1e09*(1e6); //number of particles per cubic m i.e. density const double Gamma_super=6*(2*M_PI*1e6); // characteristic decay rate of upper level used for eta calculations expressed in [1/s] const double eta = 3*lambda*lambda*N*Gamma_super/8.0/M_PI; // 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 Temperature=1000.01; // cell temperature in [K] ! make sure it is not equal to zero! + 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 + const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); + // Maxwell distribution velocities range to take in account in [m/s] + const double V_maxwell_min = -10, V_maxwell_max = -V_maxwell_min; // TODO: express via mass and Temperature // repopulation rate (atoms flying in/out the laser beam) in [1/s] const double gt=0.01 *(2*M_PI*1e6); @@ -62,12 +74,16 @@ const double R31=0.5, R32=0.5; - complex E1c, E2c, E3c, E4c; // Complex conjugated Rabi frequencies + complex E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies complex r21, r31, r41, r32, r42, r43, r44; // density matrix elements + // 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 wariables--> <benchmark /> <arguments> <!-- Rabi frequency divided by 2 in [1/s] --> @@ -91,7 +107,7 @@ <propagation_dimension> z </propagation_dimension> <transverse_dimensions> <dimension name="t" lattice="1000" domain="(-14.0e-7, 4.0e-7)" /> - <dimension name="v" lattice="10" domain="(-1000, 1000)" /> + <dimension name="v" lattice="10" domain="(V_maxwell_min, V_maxwell_max)" /> </transverse_dimensions> </geometry> @@ -110,17 +126,44 @@ </initialisation> </vector> - <!-- - <computed_vector name="E_field_avgd" dimensions="z t" type="complex"> + <!--Maxwell distribution probability p(v)--> + <computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real"> + <components>probability_v</components> + <evaluation> + <![CDATA[ + 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[ + 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</dependencies> + <dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies> <![CDATA[ - // Calculate the current normalisation of the wave function. + 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"> @@ -155,55 +198,55 @@ <!-- Below initialization assumes strong E1 and E3 which were shining for long time before we even start to look at the problem --> <!-- Precalculated by Nate via Mathematica steady state solver --> - <dependencies>E_field</dependencies> + <dependencies>E_field_avgd</dependencies> <![CDATA[ - E1c = conj(E1); - E2c = conj(E2); - E3c = conj(E3); - E4c = conj(E4); + E1ac = conj(E1a); + E2ac = conj(E2a); + E3ac = conj(E3a); + E4ac = conj(E4a); // IMPORTANT: assumes no detunings - r11 = (gt*(4*mod2(E1) + (G3 + gt)*(G3 + 2*gt))*((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) + 4*mod2(E3)*G4*(R41 - - R42)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42))); - - r13 = -((E1c*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) - + 4*mod2(E3)*G4*(R41 - - R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42))); - - r22 = (gt*(4*mod2(E3) + (G4 + gt)*(G4 + 2*gt))*((G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) + 4*mod2(E1)*G3*(-R31 + - R32)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42))); - - r24 = -((E3c*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) - + 4*mod2(E1)*G3*(-R31 + - R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42))); - - r33 = (2*mod2(E1)*gt*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) + - 4*mod2(E3)*G4*(R41 - - R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)); - - r44 = (2*mod2(E3)*gt*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) + - 4*mod2(E1)*G3*(-R31 + - R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 + - 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) + - gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 + - 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)); + r11 = (gt*(4*mod2(E1a) + (G3 + gt)*(G3 + 2*gt))*((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) + 4*mod2(E3a)*G4*(R41 - + R42)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42))); + + r13 = -((E1ac*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) + + 4*mod2(E3a)*G4*(R41 - + R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42))); + + r22 = (gt*(4*mod2(E3a) + (G4 + gt)*(G4 + 2*gt))*((G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) + 4*mod2(E1a)*G3*(-R31 + + R32)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42))); + + r24 = -((E3ac*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) + + 4*mod2(E1a)*G3*(-R31 + + R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42))); + + r33 = (2*mod2(E1a)*gt*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) + + 4*mod2(E3a)*G4*(R41 - + R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)); + + r44 = (2*mod2(E3a)*gt*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) + + 4*mod2(E1a)*G3*(-R31 + + R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 + + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) + + gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 + + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)); ]]> </initialisation> @@ -215,16 +258,26 @@ <!--<integrate algorithm="ARK45" tolerance="1e-5" interval="1.5e-2">--> <!--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="1000" tolerance="1e-5" interval="1.5e-2"> + <integrate algorithm="RK4" steps="100" tolerance="1e-5" interval="1.5e-3"> <!--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>200 200</samples> + + <filters where="step end"> + <filter> + <dependencies>E_field_avgd</dependencies> + <![CDATA[ + //printf("z=%g, t=%g\n", z, t); + //printf("E_field_avgd filter in integrator\n", z, t); + ]]> + </filter> + </filters> + <samples>100 100 100</samples> <operators> <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t"> <integration_vectors>density_matrix</integration_vectors> - <dependencies>E_field</dependencies> + <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))--> <!-- @@ -238,10 +291,10 @@ --> </boundary_condition> <![CDATA[ - E1c = conj(E1); - E2c = conj(E2); - E3c = conj(E3); - E4c = conj(E4); + E1ac = conj(E1a); + E2ac = conj(E2a); + E3ac = conj(E3a); + E4ac = conj(E4a); r21=conj(r12); r31=conj(r13); @@ -252,17 +305,17 @@ //r44=1- r11 - r22 - r33; // Equations of motions according to Nate's mathematica code - dr11_dt = gt/2 - gt*r11 + i*(-(E1*r13) - E4*r14 + E1c*r31 + E4c*r41) + G3*r33*R31 + G4*r44*R41; - dr12_dt = -(gt*r12) + i*((-(delta1 + Kvec*v) + (delta2 + Kvec*v))*r12 - E2*r13 - E3*r14 + E1c*r32 + E4c*r42); - dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1c*r11) - E2c*r12 - (delta1 + Kvec*v)*r13 + E1c*r33 + E4c*r43); - dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4c*r11) - E3c*r12 - ((delta1 + Kvec*v) - (delta2 + Kvec*v) + (delta3 + Kvec*v))*r14 + E1c*r34 + E4c*r44); - dr22_dt = gt/2 - gt*r22 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42) + G3*r33*R32 + G4*r44*R42; - dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1c*r21) - E2c*r22 - (delta2 + Kvec*v)*r23 + E2c*r33 + E3c*r43); - dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4c*r21) - E3c*r22 - (delta3 + Kvec*v)*r24 + E2c*r34 + E3c*r44); - dr33_dt = i*(E1*r13 + E2*r23 - E1c*r31 - E2c*r32) - (G3 + gt)*r33; - - dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E4c*r31 - E3c*r32 + ((delta2 + Kvec*v) - (delta3 + Kvec*v))*r34); - dr44_dt = i*(E4*r14 + E3*r24 - E4c*r41 - E3c*r42) - (G4 + gt)*r44; + dr11_dt = gt/2 - gt*r11 + i*(-(E1a*r13) - E4a*r14 + E1ac*r31 + E4ac*r41) + G3*r33*R31 + G4*r44*R41; + dr12_dt = -(gt*r12) + i*((-(delta1 + Kvec*v) + (delta2 + Kvec*v))*r12 - E2a*r13 - E3a*r14 + E1ac*r32 + E4ac*r42); + dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1ac*r11) - E2ac*r12 - (delta1 + Kvec*v)*r13 + E1ac*r33 + E4ac*r43); + dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4ac*r11) - E3ac*r12 - ((delta1 + Kvec*v) - (delta2 + Kvec*v) + (delta3 + Kvec*v))*r14 + E1ac*r34 + E4ac*r44); + dr22_dt = gt/2 - gt*r22 + i*(-(E2a*r23) - E3a*r24 + E2ac*r32 + E3ac*r42) + G3*r33*R32 + G4*r44*R42; + dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1ac*r21) - E2ac*r22 - (delta2 + Kvec*v)*r23 + E2ac*r33 + E3ac*r43); + dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4ac*r21) - E3ac*r22 - (delta3 + Kvec*v)*r24 + E2ac*r34 + E3ac*r44); + dr33_dt = i*(E1a*r13 + E2a*r23 - E1ac*r31 - E2ac*r32) - (G3 + gt)*r33; + + dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1a*r14 + E2a*r24 - E4ac*r31 - E3ac*r32 + ((delta2 + Kvec*v) - (delta3 + Kvec*v))*r34); + dr44_dt = i*(E4a*r14 + E3a*r24 - E4ac*r41 - E3ac*r42) - (G4 + gt)*r44; ]]> </operator> <operator kind="ex" constant="yes"> @@ -331,6 +384,19 @@ ]]> </sampling> </group> + + <group> + <sampling basis="t(100) " initial_sample="yes"> + <dependencies>E_field_avgd</dependencies> + <moments>I1_out_avgd I2_out_avgd I3_out_avgd I4_out_avgd</moments> + <![CDATA[ + I1_out_avgd = mod2(E1a); + I2_out_avgd = mod2(E2a); + I3_out_avgd = mod2(E3a); + I4_out_avgd = mod2(E4a); + ]]> + </sampling> + </group> </output> </simulation> |