% 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);
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')