1 files changed, 95 insertions, 0 deletions
diff --git a/xmds2/check_perturbative_approach/pp_I2.m b/xmds2/check_perturbative_approach/pp_I2.m
new file mode 100644
index 0000000..e77074c
--- /dev/null
+++ b/xmds2/check_perturbative_approach/pp_I2.m
@@ -0,0 +1,95 @@
+% We check results of  calculations via full solver vs the  case of strong E1
+% and E3. This two field create steady state coherences and populations which
+% we treat as constant (along z and t) during propagation. We keep track only of coherences
+% created by week fields E2 and E4.
+
+% In perturbative approach:
+% rho11, rho22, rho33, rho13, rho24, and rho44 dpepend  only on E1 and E2.
+% They are also treated ad time independent.
+% only r12 r14 r23 r34 treated fully since they are driven by weak fields E2 and E4
+% which are strongly z dependent as well as E1 and E2 and time.
+
+% Code below  checks result  of two  solvers. User  should double  check that
+% scripts are run with the same  command line and internal dimensions are the
+% same.
+
+% Notice that all fields feel the propagation but E1 and E2 will look like time 
+% independent in the perturbation approach.
+
+
+%% Loading fields data files
+% let's load exact solver solutions
+fname='./Nlevels_no_dopler_with_z_4wm_mg0.dat';
+fpDat = fopen(fname, 'r', 'ieee-le');
+if (fpDat < 0)
+  msg=strcat('Cannot open binary data file: ', fname);
+  disp(msg)
+  return
+end
+z_1Len = fread(fpDat, 1, 'uint32');
+z_Exact = zeros(1, z_1Len);
+z_Exact(:) = fread(fpDat, z_1Len, 'double');
+t_1Len = fread(fpDat, 1, 'uint32');
+t_Exact = zeros(1, t_1Len);
+t_Exact(:) = fread(fpDat, t_1Len, 'double');
+I1_out_1Len = fread(fpDat, 1, 'uint32');
+I1_Exact = fread(fpDat, [t_1Len, z_1Len], 'double');
+I2_out_1Len = fread(fpDat, 1, 'uint32');
+I2_Exact = fread(fpDat, [t_1Len, z_1Len], 'double');
+I3_out_1Len = fread(fpDat, 1, 'uint32');
+I3_Exact = fread(fpDat, [t_1Len, z_1Len], 'double');
+I4_out_1Len = fread(fpDat, 1, 'uint32');
+I4_Exact = fread(fpDat, [t_1Len, z_1Len], 'double');
+fclose(fpDat);
+clear fpDat z_1Len t_1Len I1_out_1Len I2_out_1Len I3_out_1Len I4_out_1Len 
+
+
+% let's load perturbative solver solutions
+fname  = './Nlevels_no_dopler_with_z_4wm_with_perturbations_mg0.dat';
+fpDat = fopen(fname, 'r', 'ieee-le');
+if (fpDat < 0)
+  msg=strcat('Cannot open binary data file: ', fname);
+  disp(msg)
+  return
+end
+z_1Len = fread(fpDat, 1, 'uint32');
+z_Perturbative = zeros(1, z_1Len);
+z_Perturbative(:) = fread(fpDat, z_1Len, 'double');
+t_1Len = fread(fpDat, 1, 'uint32');
+t_Perturbative = zeros(1, t_1Len);
+t_Perturbative(:) = fread(fpDat, t_1Len, 'double');
+I1_out_1Len = fread(fpDat, 1, 'uint32');
+I1_Perturbative = fread(fpDat, [t_1Len, z_1Len], 'double');
+I2_out_1Len = fread(fpDat, 1, 'uint32');
+I2_Perturbative = fread(fpDat, [t_1Len, z_1Len], 'double');
+I3_out_1Len = fread(fpDat, 1, 'uint32');
+I3_Perturbative = fread(fpDat, [t_1Len, z_1Len], 'double');
+I4_out_1Len = fread(fpDat, 1, 'uint32');
+I4_Perturbative = fread(fpDat, [t_1Len, z_1Len], 'double');
+fclose(fpDat);
+clear fpDat z_1Len t_1Len I1_out_1Len I2_out_1Len I3_out_1Len I4_out_1Len 
+
+%% extracting fields at the end of the cell
+I1_Exact_out=I1_Exact(:,end);
+I2_Exact_out=I2_Exact(:,end);
+I3_Exact_out=I3_Exact(:,end);
+I4_Exact_out=I4_Exact(:,end);
+
+I1_Perturbative_out=I1_Perturbative(:,end);
+I2_Perturbative_out=I2_Perturbative(:,end);
+I3_Perturbative_out=I3_Perturbative(:,end);
+I4_Perturbative_out=I4_Perturbative(:,end);
+
+%% time to uS
+t_Exact=t_Exact*1e6;
+t_Perturbative=t_Perturbative*1e6;
+
+%% comparison plot
+figure(1);
+plot(t_Exact, I2_Exact_out, '-k' , t_Perturbative, I2_Perturbative_out, '-r' )
+legend('exact', 'perturbative');
+xlabel('time (uS)');
+title('I_2 vs time');
+
+set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
+print('-color','fields_after_the_cell_I2.eps')