]> --Ndens=1e15 --gt=0.1e6 --Lcell=10.0e-2 --Temperature=1e-3 --Pwidth=0.4e-6 --delta1=0 --delta2=0 --delta3=0 --Ep1o=1e4 --Ep2o=1e2 --Ep3o=1e1 --Ep4o=1e1 --Em1o=1e4 --Em2o=1e2 --Em3o=1e1 --Em4o=1e1 --WLx=0 --WLy=0 --WLz=0 realistic_Rb_and_fields Eugeniy Mikhailov 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 than 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. &RbAtomConstantsFile; 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; // inner use variables double probability_v; // will be used as p(v) in Maxwell distribution ]]> 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; ]]> z Ep1 Em1 Ep2 Em2 Ep3 Em3 Ep4 Em4 probability_v probability_v_norm Maxwell_distribution_probabilities Ep1a Em1a Ep2a Em2a Ep3a Em3a Ep4a Em4a E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm &RbChosenRhoFile; // 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 ]]> &RbInitsFile; 100 density_matrix E_field_avgd &RbEquationsFile; Lt E_field density_matrix &RbPropEquationsFile; E_field_avgd Ip1_out Im1_out Ip2_out Im2_out Ip3_out Im3_out Ip4_out Im4_out