diff options
| -rw-r--r-- | xmds2/Nlevels_no_dopler_with_z/Makefile | 37 | ||||
| -rw-r--r-- | xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds | 228 | ||||
| -rw-r--r-- | xmds2/Nlevels_no_dopler_with_z/pp.m | 8 |
3 files changed, 273 insertions, 0 deletions
diff --git a/xmds2/Nlevels_no_dopler_with_z/Makefile b/xmds2/Nlevels_no_dopler_with_z/Makefile new file mode 100644 index 0000000..a89800a --- /dev/null +++ b/xmds2/Nlevels_no_dopler_with_z/Makefile @@ -0,0 +1,37 @@ +### -*- make -*- +### This file is part of the Debian xmds package +### Copyright (C) 2006 Rafael Laboissiere +### This file is relased under the GNU General Public License +### NO WARRANTIES! + +### This makefile can be used to build and run the XMDS examples + +XMDS_FILES = $(shell ls *.xmds) +RUN_FILES = $(patsubst %.xmds,%.run,$(XMDS_FILES)) +CC_FILES = $(patsubst %.xmds,%.cc,$(XMDS_FILES)) +XSIL_FILES = $(patsubst %.xmds,%.xsil,$(XMDS_FILES)) +M_FILES = $(patsubst %.xmds,%.m,$(XMDS_FILES)) + +XMDS = xmds2 +XSIL2GRAPHICS = xsil2graphics + +all: $(M_FILES) + +%.run: %.xmds + $(XMDS) $< + mv $(patsubst %.xmds,%,$<) $@ + +%.xsil: %.run + ./$< + +%.m: %.xsil + $(XSIL2GRAPHICS) $< + +plot: $(M_FILES) + octave pp.m + +clean: + rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf + +.PRECIOUS: %.run %.xsil %.m +.PHONY: all clean diff --git a/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds new file mode 100644 index 0000000..fc4da89 --- /dev/null +++ b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds @@ -0,0 +1,228 @@ +<?xml version="1.0"?> +<simulation xmds-version="2"> + + <name>Nlevels_no_dopler_with_z</name> + + <author>Eugeniy Mikhailov</author> + <description> + License GPL. + + Solving 4 level atom in N-field configuration, + with field propagation along spatial axis Z + no Doppler broadening + + For master equations look "Four-level 'N-scheme' in bare and quasi-dressed states pictures" + by T. Abi-Salloum, S. Meiselman, J.P. Davis and F.A. Narducci + Journal of Modern Optics, 56: 18, 1926 -- 1932, (2009). + + Present calculation matches Fig.3 (b) from the above paper. + Note that I need to double all Rabi frequencies to match the figure. + + * -------- |4> + * \ + * \ E3 -------- |3> + * \ / \ + * \ E2 / \ + * \ / \ E1 + * ------- |2> \ + * \ + * ------- |1> + * + + We moved to dimensionless units + t -> t*g ,time + z -> z*g/c , distance + rabi_frequency -> rabi_frequency/g + eta -> eta*c/g^2 , coupling constant + gij -> gij/g + Wij -> Wij/g + + where g is 1MHz rate + + We are solving + dE/dz+(1/c)*dE/dt=i*eta*rho_ij, where j llevel is higher then i + in xmds terms it looks like + dE_dz = i*eta*rhoij - L[E], here we moved t dependence to Fourier space + </description> + + <features> + <benchmark /> + <arguments> + <argument name="E1o" type="real" default_value="0.0025" /> + <argument name="E2o" type="real" default_value="2.5" /> + <argument name="E3o" type="real" default_value="3.0" /> + <argument name="delta1" type="real" default_value="0.0" /> + <argument name="delta2" type="real" default_value="0.0" /> + <argument name="delta3" type="real" default_value="0.0" /> + </arguments> + <bing /> + <fftw plan="patient" /> + <openmp /> + <auto_vectorise /> + <globals> + <![CDATA[ + const complex UpperDecayRate=10e6; // measured in s^-1 + const double c=3e8; + const double lambda=794.7e-9; //wavelength in m + const double N=1e9; //number of particles per cubic cm + const complex eta = 3*lambda*lambda*(N/1e-6)*UpperDecayRate/8.0/3.14*(c/UpperDecayRate/UpperDecayRate); //dimensionless coupling constant + + // repopulation rate (atoms flying in/out the laser beam) in MHz + const double gt=0.01/2; + // Natural linewidth of j's level in MHz + const double G3=2*2.7; + const double G4=2*3.0; + + // branching ratios + const double R41=0.5, R42=0.5; + const double R31=0.5, R32=0.5; + + //const double delta1=0; // E2 detuning in MHz + //const double delta2=0; // E2 detuning in MHz + //const double delta3=0; // E3 detuning in MHz + + //const complex E1o=0.005/2; // Rabi frequency in MHz + //const complex E2o=5.0/2; // Rabi frequency in MHz + //const complex E3o=6.0/2; // Rabi frequency in MHz + + complex E1c, E2c, E3c; // Complex conjugated Rabi frequencies + + complex r21, r31, r41, r32, r42, r43, r44; + ]]> + </globals> + </features> + + <geometry> + <propagation_dimension> z </propagation_dimension> + <transverse_dimensions> + <dimension name="t" lattice="128" domain="(-60, 60)" /> + </transverse_dimensions> + </geometry> + + <vector name="E_field" type="complex" dimensions="t"> + <components>E1 E2 E3</components> + <initialisation> + <![CDATA[ + // Initial (at starting 'z' position) electromagnetic field does not depend on detuning + // as well as time + E1=E1o; + E2=E2o; + E3=E3o; + ]]> + </initialisation> + </vector> + + <vector name="density_matrix" type="complex" dimensions="t"> + <!--<components>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34 r21 r31 r41 r32 r42 r43</components>--> + <!--<components>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34</components>--> + <components>r11 r22 r33 r12 r13 r14 r23 r24 r34</components> + <!-- + note one of the level population is redundant since + r11+r22+r33+r44=1 + so r11 is missing + --> + <initialisation> + <![CDATA[ + // Note: + // convergence is really slow if all populations concentrated at the bottom level |1> + // 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 + r11 = 1; r22 = 0; r33 = 0; + r12 = 0; r13 = 0; r14 = 0; + r23 = 0; r24 = 0; + r34 = 0; + ]]> + </initialisation> + </vector> + + <sequence> + <integrate algorithm="SIC" interval="4" steps="200"> + <samples>200</samples> + <operators> + <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t"> + <integration_vectors>density_matrix</integration_vectors> + <dependencies>E_field</dependencies> + <boundary_condition kind="left"> + <![CDATA[ + r11 = 1; r22 = 0; r33 = 0; + r12 = 0; r13 = 0; r14 = 0; + r23 = 0; r24 = 0; + r34 = 0; + ]]> + </boundary_condition> + <![CDATA[ + E1c = conj(E1); + E2c = conj(E2); + E3c = conj(E3); + + r21=conj(r12); + r31=conj(r13); + r41=conj(r14); + r32=conj(r23); + r42=conj(r24); + r43=conj(r34); + r44=1- r11 - r22 - r33; + + // Equations of motions according to Simon's mathematica code + dr11_dt = gt/2. - gt*r11 + E1*i*(-r13 + r31) + G3*r33*R31 + G4*r44*R41; + dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1*r32); + dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1*r11) - E2*r12 - delta1*r13 + E1*r33); + dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E3*r12) - (delta1 - delta2 + delta3)*r14 + E1*r34); + //dr21_dt = -(gt*r21) + i*((delta1 - delta2)*r21 - E1*r23 + E2*r31 + E3*r41); + dr22_dt = gt/2. - gt*r22 + i*(-(E2*r23) - E3*r24 + E2*r32 + E3*r42) + G3*r33*R32 + G4*r44*R42; + dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1*r21) - E2*r22 - delta2*r23 + E2*r33 + E3*r43); + dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E3*r22) - delta3*r24 + E2*r34 + E3*r44); + //dr31_dt = -((G3 + 2*gt)*r31)/2. + i*(E1*r11 + E2*r21 + delta1*r31 - E1*r33); + //dr32_dt = -((G3 + 2*gt)*r32)/2. + i*(E1*r12 + E2*r22 + delta2*r32 - E2*r33 - E3*r34); + dr33_dt = i*(E1*r13 + E2*r23 - E1*r31 - E2*r32) - (G3 + gt)*r33; + dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E3*r32 + (delta2 - delta3)*r34); + //dr41_dt = -((G4 + 2*gt)*r41)/2. + i*(E3*r21 + (delta1 - delta2 + delta3)*r41 - E1*r43); + //dr42_dt = -((G4 + 2*gt)*r42)/2. + i*(E3*r22 + delta3*r42 - E2*r43 - E3*r44); + //dr43_dt = -((G3 + G4 + 2*gt)*r43)/2. + i*(E3*r23 - E1*r41 - E2*r42 + (-delta2 + delta3)*r43); + //dr44_dt = E3*i*(r24 - r42) - (G4 + gt)*r44; + ]]> + </operator> + <operator kind="ex" constant="yes"> + <operator_names>Lt</operator_names> + <![CDATA[ + Lt = i*kt; + ]]> + </operator> + <integration_vectors>E_field</integration_vectors> + <![CDATA[ + dE1_dz = i*eta*r13 - Lt[E1]; + dE2_dz = i*eta*r23 - Lt[E2]; + dE3_dz = i*eta*r24 - Lt[E3]; + ]]> + </operators> + </integrate> + </sequence> + + + + + <!-- The output to generate --> + <output format="binary" filename="Nlevels_no_dopler_with_z.xsil"> + <group> + <sampling basis="t" initial_sample="yes"> + <dependencies>E_field</dependencies> + <moments>>I1_out I2_out I3_out</moments> + <![CDATA[ + I1_out = pow(abs(E1/E1o),2); + I2_out = pow(abs(E2/E2o),2); + I3_out = pow(abs(E3/E3o),2); + ]]> + </sampling> + </group> + </output> + +</simulation> + +<!-- +vim: ts=2 sw=2 foldmethod=indent: +--> diff --git a/xmds2/Nlevels_no_dopler_with_z/pp.m b/xmds2/Nlevels_no_dopler_with_z/pp.m new file mode 100644 index 0000000..30d4736 --- /dev/null +++ b/xmds2/Nlevels_no_dopler_with_z/pp.m @@ -0,0 +1,8 @@ +Nlevels_no_dopler_with_z + +figure(1) +imagesc(z_1, t_1, I1_out_1) +xlabel('z') +ylabel('t') +zlabel('I1_{out}') +title('I1_{out}') |