Modify xmds and nb file:

  change rf field to x direction
  make delta0 a function of z

1 files changed, 22 insertions, 14 deletions

diff --git a/xmds2/Genas_system/Genas_system.xmds b/xmds2/Genas_system/Genas_system.xmds
index e855e4f..ed059fc 100644
--- a/xmds2/Genas_system/Genas_system.xmds
+++ b/xmds2/Genas_system/Genas_system.xmds
@@ -10,12 +10,19 @@
 		Solving 4 level atom in 0->1 configuration,
 		with field propagation along spatial axis Z
 		no Doppler broadening
+		
+		Fields:
+			light field, default circularly polarized
+			rf field along x
+			static B-field along z
+		
+		M=0 upper state sublevel can be shifted independently
 
 		We are solving 
-			dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher then i.
+			dE/dz+(1/c)*dE/dt=i*eta*rho_ji,  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
+			dE_dz = i*eta*rhoji - 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.
@@ -50,13 +57,14 @@
 			<argument name="ExImo" type="real" default_value="0." />
 			<argument name="EyReo" type="real" default_value="0." />
 			<argument name="EyImo" type="real" default_value="(3.-0.001)*(2*M_PI*1.e6)" />
-			<!-- light detuning in [1/s] -->
+			<!--light detuning in [1/s] -->
 			<argument name="delta"  type="real" default_value="3.0*(2*M_PI*1e6)" />
 			<!--shift of upper-state M=0 sublevel-->
 			<argument name="delta0" type="real" default_value="1.*(2*M_PI*1e6)" />
 			<!--Static B-field Larmor frequency-->
 			<argument name="OL" type="real" default_value="4.*(2*M_PI*1e6)" />
-			<!--rf Rabi frequency-->
+			<!--normalized rf Rabi frequency:
+				Orf = (Larmor frequency due to peak rf field)/sqrt(2)-->
 			<argument name="Orf" type="real" default_value="0.1*(2*M_PI*1e6)" />
 			<!--rf frequency-->
 			<argument name="orf" type="real" default_value="4.*(2*M_PI*1e6)" />	
@@ -136,16 +144,16 @@
 						r43=conj(r34);
 
 						// Equations of motions according to Simon's mathematica code
-						dr11_dt = gt - gt*r11 + (-Ey - Ex*i)*r12 + i*(Ex + Ey*i)*r14 + i*(Exc + Eyc*i)*r21 + (-Eyc - Exc*i)*r41 + g0*(r22 + r33 + r44);
-						dr12_dt = (2*(Eyc - Exc*i)*r11 - i*((2*delta - (g0 + 2*gt)*i - 2*OL - 2*rf)*r12 - 2*(Exc + Eyc*i)*r22 + 2*(Exc - Eyc*i)*r42))/2.;
-						dr13_dt = (-g0/2. - gt - delta*i + delta0*i)*r13 + i*(Exc + Eyc*i)*r23 + (-Eyc - Exc*i)*r43;
-						dr14_dt = -(i*((2*delta - (g0 + 2*gt)*i + 2*OL + 2*rf)*r14 - 2*(Exc + Eyc*i)*r24 - 2*(Exc - Eyc*i)*(r11 - r44)))/2.;
-						dr22_dt = (Ey + Ex*i)*r12 + (Eyc - Exc*i)*r21 - (g0 + gt)*r22;
-						dr23_dt = (Ey + Ex*i)*r13 + i*(delta0 + i*(g0 + gt + i*OL) - rf)*r23;
-						dr24_dt = (Ey + Ex*i)*r14 + (Eyc + Exc*i)*r21 - (g0 + gt + 2*i*OL + 2*i*rf)*r24;
-						dr33_dt = -((g0 + gt)*r33);
-						dr34_dt = (Eyc + Exc*i)*r31 - i*(delta0 - (g0 + gt)*i + OL + rf)*r34;
-						dr44_dt = (Ey - Ex*i)*r14 + (Eyc + Exc*i)*r41 - (g0 + gt)*r44;
+						dr11_dt = gt - gt*r11 + (-Ey - Ex*i)*r12 + (-Ey + Ex*i)*r14 + (-Eyc + Exc*i)*r21 + g0*r22 + g0*r33 + (-Eyc - Exc*i)*r41 + g0*r44;
+						dr12_dt = (Eyc - Exc*i)*r11 + (-g0/2. - gt - delta*i + i*OL)*r12 + i*rf*r13 + (-Eyc + Exc*i)*r22 + (-Eyc - Exc*i)*r42;
+						dr13_dt = i*rf*r12 + (-g0/2. - gt - delta*i + i*delta0z(z))*r13 + i*rf*r14 + (-Eyc + Exc*i)*r23 + (-Eyc - Exc*i)*r43;
+						dr14_dt = (Eyc + Exc*i)*r11 + i*rf*r13 + (-g0/2. - gt - delta*i - i*OL)*r14 + (-Eyc + Exc*i)*r24 + (-Eyc - Exc*i)*r44;
+						dr22_dt = (Ey + Ex*i)*r12 + (Eyc - Exc*i)*r21 + (-g0 - gt)*r22 + i*rf*r23 - i*rf*r32;
+						dr23_dt = (Ey + Ex*i)*r13 + i*rf*r22 + (-g0 - gt - i*OL + i*delta0z(z))*r23 + i*rf*r24 - i*rf*r33;
+						dr24_dt = (Ey + Ex*i)*r14 + (Eyc + Exc*i)*r21 + i*rf*r23 + (-g0 - gt - 2*i*OL)*r24 - i*rf*r34;
+						dr33_dt = -(i*rf*r23) + i*rf*r32 + (-g0 - gt)*r33 + i*rf*r34 - i*rf*r43;
+						dr34_dt = -(i*rf*r24) + (Eyc + Exc*i)*r31 + i*rf*r33 + (-g0 - gt - i*OL - i*delta0z(z))*r34 - i*rf*r44;
+						dr44_dt = (Ey - Ex*i)*r14 - i*rf*r34 + (Eyc + Exc*i)*r41 + i*rf*r43 + (-g0 - gt)*r44;
 					]]>
 				</operator>
 				<operator kind="ex" constant="yes">