%==========================================================================
% Sodium_D_Line_DGS.m
%
% MatLab code for P589AEM HW4, problem 1: 
% Resolving the two Sodium-D Lines in a diffraction grating spectrometer.
%
%==========================================================================
clear all;
close all;

fprintf('\n');
fprintf('\n==================================================================');
fprintf('\n=== Resolving Sodium-D Lines: Diffraction Grating Spectrometer ===');
fprintf('\n===         UIUC Physics 598AEM Homework # 4, Problem 1        ===');
fprintf('\n==================================================================');

%=================================
% Wavelengths of 2 Sodium-D Lines:
%=================================
Lambda1 = 589.59; % (nm)
Lambda2 = 589.00; % (nm)
Lambda_Avg = 0.5*(Lambda1 + Lambda2);
Lambda_Dif =     (Lambda1 - Lambda2);

fprintf('\n');
fprintf('\n  Wavelength  of Sodium D line # 1  (nm) = %f',Lambda1);
fprintf('\n  Wavelength  of Sodium D line # 2  (nm) = %f',Lambda2);
fprintf('\n <Wavelength> of Sodium D lines     (nm) = %f',Lambda_Avg);
fprintf('\n  Wavelength Diff of Sodium D lines (nm) = %f',Lambda_Dif);

% Resolvance of Diffraction Grating Spectrometer:
m_order =    1;  % integer order of diffraction
N_lines = 6000;  % # of illuminated lines of diffraction grating

Resolvance_DGS = m_order*N_lines; % = Lambda/Sigma_Lambda

fprintf('\n');
fprintf('\n Order of Diffraction, m = %i',m_order);
fprintf('\n # Illuminated Lines,  N = %i',N_lines);
fprintf('\n Resolvance of DGS,    R = %i',Resolvance_DGS);

Sigma_Lambda1    = Lambda1/Resolvance_DGS; % DGS intrinsic Gaussian width of Lambda1
Sigma_Lambda2    = Lambda2/Resolvance_DGS; % DSQ intrinsic Gaussian width of Lambda2
Sigma_Lambda_Avg = 0.5*(Sigma_Lambda1 + Sigma_Lambda2); % Avg Sigma_Lambda

fprintf('\n');
fprintf('\n Sigma_Lambda1    (nm) = %f',Sigma_Lambda1);
fprintf('\n Sigma_Lambda2    (nm) = %f',Sigma_Lambda2);
fprintf('\n Sigma_Lambda_Avg (nm) = %f',Sigma_Lambda_Avg);

%==========================================================================
% Make plot of Intensity of Sodium-D Lines, as per Diff Grating Spect:
%==========================================================================
npts = 40000;

Lmda = zeros(1,npts);
I1   = zeros(1,npts);
I2   = zeros(1,npts);
Itot = zeros(1,npts);

Xlo = 587.3; % (nm)
Xhi = 591.3; % (nm)

delX = Xhi - Xlo;
  dX = delX/npts;

   X = Xlo;
  Io = 1.0; % light intensity (watts/m2)
  
  A1 = 1.0/(sqrt(2.0*pi)*Sigma_Lambda1);
  A2 = 1.0/(sqrt(2.0*pi)*Sigma_Lambda2);
  Xdif_old = 1.0e10;
for j = 1:npts;
    Lmda(j) = X;

        Y1  = 0.5*(((X - Lambda1)/Sigma_Lambda1)^2);
        Y2  = 0.5*(((X - Lambda2)/Sigma_Lambda2)^2);
      
      I1(j) = Io*A1*exp(-Y1);
      I2(j) = Io*A2*exp(-Y2);
      
    Itot(j) = I1(j) + I2(j);
    
    % Store wavelength index closest to Lambda_Avg:
    Xdif = abs(X - Lambda_Avg);
    if (Xdif < Xdif_old);
        Xdif_old = Xdif;
           j_mid = j;
    end
    
    X = X + dX;
end

fprintf('\n');
fprintf('\n Wavelength (nm) of I1-I2 mid-point  = %f',Lmda(j_mid));
fprintf('\n Total Intensity at Lambda_Avg point = %f',Itot(j_mid));

figure(01);
plot(Lmda,I1,'r',Lmda,I2,'g',Lmda,Itot,'b');
axis tight;
grid on;
xlabel('Wavelength (nm)');
ylabel('Intensity');
legend('I(D1)','I(D2)','Itot');

figure(02);
semilogy(Lmda,I1,'r',Lmda,I2,'g',Lmda,Itot,'b');
axis tight;
grid on;
xlabel('Wavelength (nm)');
ylabel('Intensity');
legend('I(D1)','I(D2)','Itot');

%==========================================================================
% Make 3-D plot of Intensity of Sodium-D Lines, as per Diff Grating Spect:
%==========================================================================
 mpts =  4000;
 ipts =  1000; % n.b. 10x improves accuracy of determinine Sigma for Rayleigh Criteria - eats memory though...
%ipts =   100;

      lmda = zeros(ipts,mpts);
   siglmda = zeros(ipts,mpts);
      i1   = zeros(ipts,mpts);
      i2   = zeros(ipts,mpts);
      itot = zeros(ipts,mpts);

      itot_mid = zeros(1,ipts);
       sig_mid = zeros(1,ipts);
       i_ratio = zeros(1,ipts);
      itot_max = zeros(1,ipts);
       i1_max  = zeros(1,ipts);

 io = 1.0; % light intensity (watts/m2)

Rayleigh_Criterion = 0.735; % = Itot_min(@ saddle pt)/Itot_max
          deli_old = 1.0e10;

     xlo = 587.3; % (nm)
     xhi = 591.3; % (nm)

    delx = xhi - xlo;
      dx = delX/mpts;

 siglo   = 0.05; % (nm)
 sighi   = 0.55; % (nm)
 delsig  = sighi - siglo;
   dsig  = delsig/ipts;
  
idel_old = 1.0e10;

     sig = siglo;
for i = 1:ipts; % sigma lambda loop <= n.b. sig = *total* sigma!
    x = xlo;
  
    a = 1.0/(sqrt(2.0*pi)*sig);
  xdif_old = 1.0e10;
  max_itot = 0.0;
    for j = 1:mpts; % lambda loop
     siglmda(i,j) = sig;
        lmda(i,j) =   x;

            Y1  = 0.5*(((x - Lambda1)/sig)^2);
            Y2  = 0.5*(((x - Lambda2)/sig)^2);
      
        i1(i,j) = io*a*exp(-Y1);
        i2(i,j) = io*a*exp(-Y2);
      
            itot(i,j) = i1(i,j) + i2(i,j);
            
            % find max value of itot(i,j)
            if (itot(i,j) > max_itot);
                   max_itot = itot(i,j); % update new value of max_itot
                itot_max(i) = itot(i,j); % store the  biggest  max_itot
            end
    
        % Store wavelength index closest to Lambda_Avg:
        xdif = abs(x - Lambda_Avg);
        if (xdif < xdif_old);
            xdif_old = xdif;
               j_mid = j;
        end
    
        x = x + dx;
    end
    
    % Store total intensity at Lambda_Avg point for each sigma_lambda value:
     sig_mid(i) = sig;
    itot_mid(i) = itot(i,j_mid);
    
    % Determine sigma-tot point where itot_mid = itot_max:
    idel = abs(itot_mid(i) - itot_max(i));
    if (idel < idel_old);
        idel_old = idel;
        sig_mtch =  sig;
       itot_mtch = itot_max(i);
          i_mtch =    i; % index i for this point
    end
    
    % Calculate ratio itot_mid/itot_max for Rayleigh Criterion determination:
    i_ratio(i) = itot_mid(i)/itot_max(i);
     
     % store index of point where i_ratio = 0.735 (Rayleigh Criterion)
     deli = abs(i_ratio(i) - Rayleigh_Criterion);
     if (deli < deli_old);
         deli_old = deli;
           i_rayl = i;
     end
    
    sig = sig + dsig;
end

fprintf('\n');
fprintf('\n=======================================');
fprintf('\n=== Rayleigh Criterion Calculations ===');
fprintf('\n=======================================');

fprintf('\n');
fprintf('\n Sigma-tot point (nm)   where [Itot-min(@ Saddle Pt.)=Itot-max]: %f', sig_mtch);
fprintf('\n Itot-max  point (W/m2) where [Itot-min(@ Saddle Pt.)=Itot-max]: %f',itot_mtch);

fprintf('\n');
fprintf('\n Rayleigh Criterion: Ratio Itot_mid/Itot_max = 0.735, can *just* resolve the two Gaussian peaks...');
fprintf('\n Sigma-tot (nm) at Rayleigh Criterion saddle point = %f', siglmda(i_rayl));
fprintf('\n Itot-max(W/m2) at Rayleigh Criterion saddle point = %f',itot_max(i_rayl));
fprintf('\n Itot-mid(W/m2) at Rayleigh Criterion saddle point = %f',itot_mid(i_rayl));
fprintf('\n Ratio Itot-mid/Itot-max at Rayleigh Criterion pnt = %f', i_ratio(i_rayl));

xg0(1) =  sig_mtch;
yg0(1) = itot_mtch;

xg1(1) =  sig_mtch;
xg1(2) =  sig_mtch;
yg1(1) = 0.0;
yg1(2) = 8.0;

xg2(1) = siglo;
xg2(2) = sighi;
yg2(1) = itot_mtch;
yg2(2) = itot_mtch;

%xg2(1) = siglo;
%xg2(2) = sighi;
%yg2(1) = itot_max(i_rayl);
%yg2(1) = itot_max(i_rayl);

xg3(1) = siglo;
xg3(2) = sighi;
yg3(1) = Rayleigh_Criterion;
yg3(2) = Rayleigh_Criterion;

xg4(1) = siglmda(i_rayl);
xg4(2) = siglmda(i_rayl);
yg4(1) = 0.0;
yg4(2) = 1.0;

% n.b. accuracy/resolution {= 0.005 nm} on Sigma-tot determination is
% affected by # steps {= 1000} in sigma loop from 0.05 < sigma-tot < 0.55 !!!

% Now calculate upper limit on Sigma_Lambda_CCD due to Gaussian electronic noise linear CCD array:
sig_ccd = sqrt((siglmda(i_rayl))^2 - (Sigma_Lambda_Avg)^2);

fprintf('\n');
fprintf('\n Sodium D Lines: Sigma_Lambda_Avg (nm)      = %f',Sigma_Lambda_Avg);
fprintf('\n Sigma-CCD due to CCD electronic noise (nm) = %f',sig_ccd);

figure(11);
surf(lmda,siglmda,itot);
axis tight;
view(35.0,60.0);
shading interp;
xlabel('Wavelength (nm)');
ylabel('Sigma-Lambda-tot (nm)');
zlabel('Itot (W/m2)');
title ('DGS Itot vs. Lambda vs. Sigma-Lambda-tot');

figure(12);
plot(sig_mid,itot_mid,'b',sig_mid,itot_max,'m',xg0,yg0,'bo',xg0,yg0,'m*',xg1,yg1,'k',xg2,yg2,'k');
axis tight;
grid on;
xlabel('Sigma-Lambda-tot (nm)');
ylabel('Itot (W/m2)');
title ('Itot-mid (at Lambda-avg saddle point) and Itot-max vs. Sigma-Lambda-tot');
legend('Itot-mid','Itot-max');

figure(13);
plot(sig_mid,i_ratio,'b',xg3,yg3,'k',xg4,yg4,'k',sig_mid(i_rayl),i_ratio(i_rayl),'m*',sig_mid(i_rayl),i_ratio(i_rayl),'bo');
axis tight;
grid on;
xlabel('Sigma-Lambda-tot (nm)');
ylabel('Iratio');
title ('Rayleigh Criterion: Intensity Ratio (Itot-mid/I1max) vs. Sigma-Lambda-tot');

%==========================================================================
beep();
fprintf('\n');
fprintf('\n Sodium_D_Line_DGS calculations finished \n');
%==========================================================================