Merge branch 'testing' of qo.physics.wm.edu:/eit_filter_simulations into testing

Conflicts:
	config.m

7 files changed, 74 insertions, 49 deletions

diff --git a/config.m b/config.m
index 157cee8..d914e7e 100644
--- a/config.m
+++ b/config.m
@@ -3,10 +3,10 @@
 % Always run this file first.
 %
 
-decay_ground = [0.00001e6];	% Ground state decay
-decay_upper = [6e6];	% Excited state decay
-drive = [linspace(1e4,2e4,10)];	% Drive/control field
-drive = drive(1:9);	% Choose which drives you want
+decay_ground = [0.00001];	% Ground state decay
+decay_upper = [6];	% Excited state decay
+drive = [linspace(1e-2,1,100)];	% Drive/control field
+drive = drive;	% Choose which drives you want 
 detuning1 = [0];	% Single photon detuning
 dephase_bc = [0];	% Ground state dephasing
 
diff --git a/octave-workspace b/octave-workspace
new file mode 100644
index 0000000..8622b89
--- /dev/null
+++ b/octave-workspace
diff --git a/source/eit.xmds b/source/eit.xmds
index 81106d5..592c0fc 100644
--- a/source/eit.xmds
+++ b/source/eit.xmds
@@ -35,32 +35,38 @@
         real delta1;
         real dephase_bc;
         complex Pac, Pba;
+        complex probe_c, drive_c;
 
-        const double c = 3.e8;
-        const double lambda=794.7e-9; //wavelength in m
-        const double N=1e10*(1e6);    //number of particles per cubic m i.e. density
-        double eta;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+        const double c = 3.e8; // m/s^2
+        const double lambda=794.7e-9; // wavelength in m
+        const double N=1e9*(1e6);    //number of particles per cubic m 
+        const double Gamma_super = 6;
+        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)]
         
       ]]>
      </globals>
      <validation kind="run-time"/>
      <arguments>
-      <argument name="drive_arg" type="real" default_value="0.1e6" />
-      <argument name="decay_upper_arg" type="real" default_value="10.0e6" />
-      <argument name="decay_ground_arg" type="real" default_value="0.001e6" />
-      <argument name="domain_min" type="real" default_value="-1" />
-      <argument name="domain_max" type="real" default_value="1" />
+      <argument name="decay_upper_arg" type="real" default_value="1e1" />
+      <argument name="decay_ground_arg" type="real" default_value="1e-3" />
+      <argument name="drive_arg" type="real" default_value="1e-1" />
+      <argument name="probe_arg" type="real" default_value="1e-4" />
+      <!-- Expected detuning width needed is: 
+      decay_ground + (drive^2 + probe^2)/decay_upper 
+      -->
+      <argument name="detuning_min" type="real" default_value="-10" />
+      <argument name="detuning_max" type="real" default_value="10" />
       <argument name="delta1_arg" type="real" default_value="0" />
       <argument name="dephase_bc_arg" type="real" default_value="0" />
-      <argument name="final_time" type="real" default_value="10000" />
+      <argument name="start_time" type="real" default_value="0" />
+      <argument name="end_time" type="real" default_value="100" />
       <![CDATA[
+        
         decay_upper = decay_upper_arg;
         decay_ground = decay_ground_arg;
         delta1 = delta1_arg;
         dephase_bc = dephase_bc_arg;
-        r_ca = decay_upper - i*delta1;
-
-        eta = 3*lambda*lambda*N*decay_upper/8.0/M_PI;
+       
       ]]>
     </arguments>
     <bing />
@@ -81,8 +87,8 @@
       <dimension name="detuning" lattice="256"  domain="(domain_min, domain_max)" />
       -->
       
-      <dimension name="t"   lattice="300"   domain="(-2.0e-6, 4.0e-6)" />
-      
+      <dimension name="t"   lattice="1000"   domain="(start_time, end_time)"  />
+      <dimension name="detuning" lattice="100"  domain="(detuning_min, detuning_max)" />
     </transverse_dimensions>
   </geometry>
 
@@ -93,22 +99,23 @@
       // Initial (at starting 'z' position) electromagnetic field does not depend on detuning
       // as well as time
       drive = drive_arg;
-      probe = 0.0001e6; //*exp(-pow( ((t-0.0)/1e-6),2) );
+      probe = probe_arg; //*exp(-pow( ((t-0.0)/1e-6),2) );
       ]]>
     </initialisation>
   </vector>
 
   <!-- A 'vector' describes the variables that we will be evolving. -->
   <vector name="Gamma" type="complex">
-    <components> r_ab r_cb </components>
+    <components> r_ab r_cb r_ca</components>
     <initialisation>
       <![CDATA[
-      r_ab = decay_upper + i*(delta1);
-      r_cb = decay_ground + dephase_bc;
+      r_ab = decay_upper + i*(delta1 + detuning);
+      r_cb = decay_ground + dephase_bc + i*detuning;
+      r_ca = decay_upper - i*delta1;
       ]]>
     </initialisation>
   </vector>
-  <vector name="population" type="complex" initial_space="t">
+  <vector name="population" type="complex" dimensions="t detuning">
     <components>
       Paa Pbb Pcc Pab Pca Pcb 
     </components>
@@ -127,7 +134,7 @@
     and what algorithm we want to use.
     -->
     <integrate algorithm="ARK45" tolerance="1e-5" interval="7e-2">
-      <samples>10 10</samples>
+      <samples>100 100</samples>
       <operators>
         <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
           <integration_vectors>population</integration_vectors>
@@ -143,18 +150,27 @@
 
           Pac = conj(Pca);
           Pba = conj(Pab);
+          probe_c = conj(probe);
+          drive_c = conj(drive);
+
+          
 
-          dPcc_dt = i*drive*Pac - i*drive*Pca + decay_upper*Paa - decay_ground*Pcc + decay_ground*Pbb;
+          
 
-          dPbb_dt = i*probe*Pab - i*probe*Pba + decay_upper*Paa - decay_ground*Pbb + decay_ground*Pcc;
+          dPcc_dt = i*drive_c*Pac - i*drive*Pca + decay_upper*Paa - decay_ground*Pcc + decay_ground*Pbb;
 
-          dPaa_dt = -i*probe*Pab + i*probe*Pba - i*drive*Pac + i*drive*Pca - 2*decay_upper*Paa;
+          dPbb_dt = i*probe_c*Pab - i*probe*Pba + decay_upper*Paa - decay_ground*Pbb + decay_ground*Pcc;
+
+          dPaa_dt = -i*probe_c*Pab + i*probe*Pba - i*drive_c*Pac + i*drive*Pca - 2*decay_upper*Paa;
 
           dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb;
 
-          dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb;
+          dPca_dt = -r_ca*Pca + i*drive_c*(Paa-Pcc) - i*probe_c*Pcb;
+
+          dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive_c*Pab;
+
+          //printf("%0.15f\n",detuning);
 
-          dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab;
           ]]>
         </operator>
 
@@ -168,8 +184,8 @@
         <integration_vectors>E_field</integration_vectors>
         <dependencies>population</dependencies>
           <![CDATA[
-          dprobe_dz = i*eta*conj(Pca) -Lt[probe] ;
-          ddrive_dz = i*eta*conj(Pba) -Lt[drive] ;
+          dprobe_dz = (i*eta*Pab -Lt[probe]) ;
+          ddrive_dz = (i*eta*Pac -Lt[drive]) ;
           ]]>
       </operators>
     </integrate>
@@ -177,7 +193,7 @@
 
   <!-- This part defines what data will be saved in the output file -->
   <output format="hdf5" >
-    <sampling_group basis="t(300)" initial_sample="yes">
+    <sampling_group basis="t detuning" initial_sample="yes">
       <dependencies>E_field</dependencies>
       <moments>I1_out I2_out</moments>
       <![CDATA[
@@ -185,7 +201,7 @@
       I2_out = mod2(drive);
       ]]>
     </sampling_group>
-    <sampling_group basis="t(100)" initial_sample="yes">
+    <sampling_group basis="t detuning" initial_sample="yes">
       <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments>
       <dependencies>population</dependencies>
       <![CDATA[
@@ -202,10 +218,11 @@
         PcbR = Pcb.Re();
         PcbI = Pcb.Im();
       ]]> -->
-      <moments>PabI</moments>
+      <moments>PabI PabR</moments>
       <dependencies>population</dependencies>
       <![CDATA[
         PabI = Pab.Im();
+        PabR = Pab.Re();
       ]]>
     </sampling_group>
   </output>
diff --git a/source/simulate.py b/source/simulate.py
index 47f7633..9639c6a 100644
--- a/source/simulate.py
+++ b/source/simulate.py
@@ -101,8 +101,8 @@ def simulate((index, decays_upper, decays_ground, drives, detunings1, dephasings
 						domain_min = -1*domain_max
 						final_time = 1/expected_width*time_multiplier
 						
-						options = '--drive_arg=%0.20f --decay_upper_arg=%f --decay_ground_arg=%f --delta1_arg=%f --dephase_bc_arg=%f --domain_min=%f --domain_max=%f --final_time=%f'
-						values = (drive, decay_upper, decay_ground, detuning1, dephase_bc, domain_min, domain_max, final_time)
+						options = '--drive_arg=%0.20f --decay_upper_arg=%f --decay_ground_arg=%f --delta1_arg=%f --dephase_bc_arg=%f --detuning_min=%f --detuning_max=%f'
+						values = (drive, decay_upper, decay_ground, detuning1, dephase_bc, domain_min, domain_max)
 						options = options % values
 						run(target, options)
 						
@@ -154,4 +154,4 @@ def main():
 	p.map(simulate, options)
 
 if __name__ == '__main__':
-	main()
\ No newline at end of file
+	main()
diff --git a/source/xmds/eit_no_z_no_doppler.xmds b/source/xmds/eit_no_z_no_doppler.xmds
index d0ef3f0..5df5d63 100644
--- a/source/xmds/eit_no_z_no_doppler.xmds
+++ b/source/xmds/eit_no_z_no_doppler.xmds
@@ -34,6 +34,7 @@
         real probe = .0001;
         real delta1;
         real dephase_bc;
+        complex Pac, Pba;
         
       ]]>
      </globals>
@@ -53,7 +54,7 @@
         decay_ground = decay_ground_arg;
         delta1 = delta1_arg;
         dephase_bc = dephase_bc_arg;
-        r_ca = decay - i*delta1;
+        r_ca = decay_upper - i*delta1;
       ]]>
     </arguments>
     <bing />
@@ -71,10 +72,7 @@
     <transverse_dimensions>
       
       <dimension name="detuning" lattice="256"  domain="(domain_min, domain_max)" />
-      
-      <!--
-      <dimension name="t"   lattice="1000"   domain="(-2.0e-6, 4.0e-6)" />
-      -->
+
     </transverse_dimensions>
   </geometry>
 
@@ -88,7 +86,7 @@
       ]]>
     </initialisation>
   </vector>
-  <vector name="population" type="complex" initial_space="t">
+  <vector name="population" type="complex" dimensions="detuning">
     <components>
       Paa Pbb Pcc Pab Pca Pcb 
     </components>
diff --git a/source/xmds/eit_with_z_no_doppler_mik.xmds b/source/xmds/eit_with_z_no_doppler_mik.xmds
index 2f8c09c..a63c169 100644
--- a/source/xmds/eit_with_z_no_doppler_mik.xmds
+++ b/source/xmds/eit_with_z_no_doppler_mik.xmds
@@ -87,7 +87,8 @@
 	<geometry>
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
-			<dimension name="t"   lattice="1000"   domain="(-2.0e-6, 4.0e-6)" />
+			<dimension name="t"   lattice="10"   domain="(-2.0e-6, 1e-5)" />
+			<dimension name="detuning" lattice="1" domain="(-1,1)" />
 		</transverse_dimensions>
 	</geometry>
 
@@ -140,7 +141,7 @@
 		<!--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>
+			<samples>10 10</samples>
 			<operators>
         <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
 					<integration_vectors>density_matrix</integration_vectors>
@@ -166,10 +167,13 @@
 						r43=conj(r34);
 						//r44=1- r11 - r22 - r33;
 
+
+						
+
 						// Equations of motions according to Simon's mathematica code
 						dr11_dt = gt/2. - gt*r11 + G3*r33*R31 + G4*r44*R41 + i*(-(E1*r13) + E1c*r31);
 						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(delta1*r13) - E1c*r11 - E2c*r12 + E1c*r33);
+						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-((delta1+detuning)*r13) - E1c*r11 - E2c*r12 + E1c*r33);
 						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-((delta1 - delta2 + delta3)*r14) - E3c*r12 + E1c*r34);
 						dr22_dt = gt/2. - gt*r22 + G3*r33*R32 + G4*r44*R42 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42);
 						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(delta2*r23) - E1c*r21 - E2c*r22 + E2c*r33 + E3c*r43);
@@ -202,7 +206,7 @@
 	<!-- The output to generate -->
 	<output format="hdf5">
 		<group>
-      <sampling basis="t(1000)" initial_sample="yes">
+      <sampling basis="t detuning" initial_sample="yes">
 				<dependencies>E_field</dependencies>
 				<moments>I1_out I2_out I3_out</moments>
 				<![CDATA[
@@ -214,7 +218,7 @@
 		</group>
 
 		<group>
-      <sampling basis="t(100)" initial_sample="yes">
+      <sampling basis="t detuning" initial_sample="yes">
 				<dependencies>density_matrix</dependencies>
 				<moments>
 					r11_out r22_out r33_out r44_out 
diff --git a/view/plotSimulations.m b/view/plotSimulations.m
index de5344c..6197c75 100644
--- a/view/plotSimulations.m
+++ b/view/plotSimulations.m
@@ -7,6 +7,12 @@ function plotSimulations(drives, y1, y2)
 	ylabel('(log scale)')
 	title(sprintf('%s and %s Vs Drive Intensity', inputname(2), inputname(3)))
 	legend(inputname(2), inputname(3))
+
+	% When plotting contrast, the y-axis upper limit is the maximum value (1)
+	% Since contrast approaches one quickly and becomes a practically straight line, it's hard to see
+	% Here we increase the y-axis upper limit to prevent the contrast line from overlapping with the axis
+	limits = ylim;
+	ylim([limits(1), limits(2)*10]);
 	hold off
 
 end
\ No newline at end of file