1 files changed, 103 insertions, 0 deletions
diff --git a/xmds2/check_perturbative_approach/pp.m b/xmds2/check_perturbative_approach/pp.m
new file mode 100644
index 0000000..a80c9c4
--- /dev/null
+++ b/xmds2/check_perturbative_approach/pp.m
@@ -0,0 +1,103 @@
+% 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.
+
+% Code below  checks result  of two  solvers. User  shuold 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
+fpDat = fopen('../Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z_mg0.dat', 'r', 'ieee-le');
+if (fpDat < 0)
+  disp('Cannot open binary data file: Nlevels_no_dopler_with_z_mg0.dat')
+  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
+fpDat = fopen('../Nlevels_no_dopler_with_z_4wm_with_perturbations/Nlevels_no_dopler_with_z_mg0.dat', 'r', 'ieee-le');
+if (fpDat < 0)
+  disp('Cannot open binary data file: Nlevels_no_dopler_with_z_mg0.dat')
+  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);
+subplot(2,2,1);
+plot(t_Exact, I1_Exact_out, '-k' , t_Perturbative, I1_Perturbative_out, '-r' )
+legend('exact', 'perturbative');
+xlabel('time (uS)');
+title('I_1 vs time');
+
+subplot(2,2,2);
+plot(t_Exact, I2_Exact_out, '-k' , t_Perturbative, I2_Perturbative_out, '-r' )
+legend('exact', 'perturbative');
+xlabel('time (uS)');
+title('I_2 vs time');
+
+subplot(2,2,3);
+plot(t_Exact, I3_Exact_out, '-k' , t_Perturbative, I3_Perturbative_out, '-r' )
+legend('exact', 'perturbative');
+xlabel('time (uS)');
+title('I_3 vs time');
+
+subplot(2,2,4);
+plot(t_Exact, I4_Exact_out, '-k' , t_Perturbative, I4_Perturbative_out, '-r' )
+legend('exact', 'perturbative');
+xlabel('time (uS)');
+title('I_4 vs time');
+
+print('-color','fields_after_the_cell.eps')