zOpenWIndow('Sas');
uiwait(msgbox(['Switch to ZEMAX and save the lens as "Primary Reversed.zmx" in the directory ' NarcDir], 'Save Request', 'custom',NarcIcon, hot(64), 'modal'));
zOpenWindow('Gho');
uiwait(msgbox('Switch to ZEMAX and generate single bounce ghost lenses for all surfaces. Remember to specify Save Files and the correct coating', 'Ghost Lens Generation', 'custom', NarcIcon, hot(64), 'modal'));

% The next job is to perform a vignetting calculation and save text for each ghost lens
% Start by getting a list of all ghost lenses
GhostFiles = dir([NarcDir 'gh*.zmx']);
% Loop and load each ghost lens in turn

BusyMessage = ZEMAXBusy; pause(1);
nuFID = fopen([NarcDir 'nu.txt'], 'wt');
for i=1:size(GhostFiles,1)
    zLoadFile([NarcDir GhostFiles(i).name]);  % Load up the next ghost lens file
    LenSys = zsGetSystem;                     % Get system operating parameters esp. number of surfaces
    zPushLens(3);                             % Push the lens into the foreground
    GhostSurf = GhostFiles(i).name(3:5);      % Extract the ghosting surface from the filename
    NoNoSurf  = GhostFiles(i).name(6:8);      % This must be 000, otherwise user has generated double bounce lenses
    if (~strcmp(NoNoSurf, '000'))             % The user seems to have generated double bounce lenses
        uiwait(msgbox('Encountered double bounce ghost lenses. Restart the analysis.', 'Double Bounce Error', 'error', 'modal'));
        delete([NarcDir 'GH*.ZMX']);          % Delete all ghost lenses in this directory
        delete(BusyMessage);                  % Delete the busy message
        Results.Status = -7;                  % Give up
        return;
    end;
    nuDataFile = ['nu' GhostSurf '.txt'];     % Synthesize filename of vignetting data with cold stop
    nu_cDataFile = ['cnu' GhostSurf '.txt'];  % Filename of vignetting data without cold stop
    % Compute transmission for the axial ray
    PolTraceData = zGetPolTrace(PrimaryWave, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0);
    tau = PolTraceData(2); % This is axial system transmission. An enhanced version could do better.
    fprintf(nuFID, '%s %11.9f\n', nuDataFile, tau);       % Write filename and transmission data into file for later
    nuDataFile = [NarcDir nuDataFile];                    % Synthesize full pathname of vignetting data with cold stop
    nu_cDataFile = [NarcDir nu_cDataFile];                % Synthesize full pathname of vignetting data without cold stop
    if (~exist(nuDataFile, 'file')),                      % Check if the vignetting datafile already exists
      zGetTextFile(nuDataFile, 'Vig', SettingsFile, 1);   % if not, generate vignetting data for this ghost lens
    end;
    zSetAperture(LenSys.numsurfs - 1, 0, 0, 0, 0, 0, ''); % Now remove the cold stop aperture
    zPushLens(3);                                         % Push the lens into the foreground
    if (~exist(nu_cDataFile, 'file')),                    % Check if the vigentting data file exists
      zGetTextFile(nu_cDataFile, 'Vig', SettingsFile, 1); % if not, generate vignetting data for this ghost lens with cold stop removed
    end;
end
fclose(nuFID);
delete(BusyMessage);

% Here are the pixel resolutions at which to plot the qualitative narcissus map
xres = 640;
yres = 480;

% Here are the wavelengths over which we will integrate
deltalambda = 0.02;
lambda1 = 3.0;
lambda2 = 5.0;
lambda = lambda1:deltalambda:lambda2;
lambda = lambda';

% Compute Black body curves for ambient scene, dewar and housing
BB_a = planck(lambda, T_a);
BB_d = planck(lambda, T_d);
BB_h = planck(lambda, T_h);

% Read in the relative illumination (irradiance) rho
% Please check that the number of headerlines is correct

[yrho, rho] = textread([NarcDir 'rho.txt'], '%f %f', 'headerlines', 13);

% Now read in all the filename and transmittance data from nu.txt

[names, tau] = textread([NarcDir, 'nu.txt'], '%9c %f'); % Filenames of vignetting data with cold stop
names_c = cat(2,repmat('c',size(names,1),1),names);     % Filenames of vignetting data wihtout cold stop
% Now read in all of the vignetting curves
for i=1:size(names)
    [y, nu(:,i)] = textread([NarcDir names(i,:)], '%f %f', 'headerlines', 12);     % Read vignetting data with cold stop
    [y, nu_c(:,i)] = textread([NarcDir names_c(i,:)], '%f %f', 'headerlines', 12); % Read vignetting data without cold stop
end

% It is possible that rho has been computed on a less dense y field sampling. We must therefore resample rho.
rho = interp1(yrho, rho, y, 'spline');

% Now compute and sum all the irradiance terms

% First the contribution from the ambient scene, with and without the effect of rho
H_arho = tau_S * trapz(lambda, BB_a) * rho / (4 * F * F);
H_a = tau_S * trapz(lambda, BB_a) * ones(size(rho)) / (4 * F * F);
% Next comes the direct contribution from the cold dewar walls;
H_d = 2 * sqrt(1 - 1 / (4 * F * F)) * trapz(lambda, BB_d) * ones(size(H_a));
% The following contribution is the "cold return", or cold straylight from single bounce reflections
H_id = repmat(tau',size(nu,1),1).*nu*trapz(lambda, BB_d) / (4 * F * F);
sigmaH_id = sum(H_id,2);
% Secondly, we have the mixed return from the outside of the cold stop
H_idh = repmat(tau',size(nu,1),1).*(nu_c - nu)*trapz(lambda, (r_c * BB_h + (1 - r_c) * BB_d)) / (4 * F * F);
sigmaH_idh = sum(H_idh,2);
% Lastly, the warm return, assumed to come from the housing
H_ih = repmat(tau',size(nu,1),1).*(1-nu_c)*trapz(lambda, BB_h) / (4 * F * F);
sigmaH_ih = sum(H_ih,2);

% Here comes the final irradiance sum, with and without the effect of rho

H_rho = H_arho + H_d + sigmaH_id + sigmaH_idh + sigmaH_ih;
H = H_a + H_d + sigmaH_id + sigmaH_idh + sigmaH_ih;

% and here come the contributions per surface

H_i = repmat(H_a, 1, size(H_id,2)) + repmat(H_d, 1, size(H_id,2)) + H_id + H_idh + H_ih;

% And finally we get to NITD

% First compute planck curve for radiation one kelvin above ambient
BB_a1 = planck(lambda, T_a + 1);

% Compute Detector irradiance delta T
TD = (H_rho-H(1))./(tau_S * rho * (trapz(lambda, BB_a1) - trapz(lambda, BB_a))/(4 * F * F));

% Total NITD computed without the effect of rho
%NITD = (H-H(1))./(tau_S * ones(size(rho)) * (trapz(lambda, BB_a1) - trapz(lambda, BB_a))/(4 * F * F));
NITD = (H-H(1))/(tau_S *(trapz(lambda, BB_a1) - trapz(lambda, BB_a))/(4 * F * F));

% Compute the NITD contributions per surface 
NITD_i = (H_i - repmat(H_i(1,:), size(H_i,1),1))/(tau_S *(trapz(lambda, BB_a1) - trapz(lambda, BB_a))/(4 * F * F));

% Find worst offending surfaces

sigmaNITD_i = sum(NITD_i);
sigmaNITD_i = sigmaNITD_i';
Worst = flipdim(sortrows([str2num(names(:,3:5)) sigmaNITD_i],2),1); % Rank surfaces on worst total integrated NITD

% Compute the maximum gradient of NITD contribution per surface in kelvins
% per image diagonal
for Col = 1:size(NITD_i,2)
  slup = 2 * max(gradient(NITD_i(:,Col),mean(diff(y/max(y)))));    % maximum slope is most upward
  sldn = -2 * min(gradient(NITD_i(:,Col),mean(diff(y/max(y)))));   % minimum slope is most downward
  maxgradNITD_i(1,Col) = max([slup, sldn]);                        % maximum of either upward or downward slope
end;

WorstSlope = flipdim(sortrows([str2num(names(:,3:5)) maxgradNITD_i'],2),1); % Rank surfaces on maximum gradient of NITD contribution

% Read in the YNI contributions
[YNIsurf, YNIf] = textread([NarcDir, 'YNIf.txt'], '%f %f', 'headerlines', 11);
% Reverse order as for ghost analysis lenses
YNIf = flipud(YNIf);
Worst(:,3) = YNIf(Worst(:,1));
WorstSlope(:,3) = YNIf(WorstSlope(:,1));
% Display Worst offenders with rank and YNI contribution
Rank = (1:size(Worst,1))';
Worst = cat(2,Rank,Worst);
Rank = (1:size(WorstSlope,1))';
WorstSlope = cat(2,Rank,WorstSlope);
disp(Worst);         % Worst integrated NITD contributions
disp(WorstSlope);    % Worst NITD slope contributions

h1 = figure;
set(h1, 'Tag', 'NITDps');
plot(y,NITD_i);
grid;
xlabel('Field Position (mm)');
ylabel('NITD per Surface (K)');
legend(names(:,4:5));
title(['NITD per Surface for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K    r_c = ' num2str(r_c*100) '%']);

h2 = figure;
set(h2, 'Tag', 'TNITD');
plot(y,NITD);
grid;
xlabel('Field Position (mm)');
ylabel('Total NITD (K)');
title(['Total NITD for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K    r_c = ' num2str(r_c*100) '%']);

h3 = figure;
set(h3, 'Tag', 'IrrUni');
plot(y,TD);
grid;
xlabel('Field Position (mm)');
ylabel('TD (K)');
title(['Irradiance Uniformity for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K    r_c = ' num2str(r_c*100) '%']);


% Compute the NITD and TD maximum gradients in kelvins per image diagonal
slup = 2 * max(gradient(NITD,mean(diff(y/max(y)))));   % maximum slope is most upward
sldn = -2 * min(gradient(NITD,mean(diff(y/max(y)))));  % minimum slope is most downward
maxgradNITD = max([slup, sldn]);                       % find the steepest slope

slup = 2 * max(gradient(TD  ,mean(diff(y/max(y)))));   % maximum slope is most upward
sldn = -2 * min(gradient(TD  ,mean(diff(y/max(y)))));  % minimum slope is most downward
maxgradTD = max([slup, sldn]);                         % find the steepest slope

% Display the maximum gradient values
disp(['Maximum Gradient of NITD is ' num2str(maxgradNITD) ' K per Image Diagonal']);
disp(['Maximum Gradient of Image Uniformity is ' num2str(maxgradTD) ' K per Image Diagonal']);

% Now also make narcissus image map for "qualitative" impact of this total NITD and uniformity curve
% We find the distance of the pixel centre from the image centre and set the brightness according to
% the fraction of the image semi-diagonal.
% First find the length of the image semidiagonal in pixels.
semidiag = round(sqrt((xres/2)^2 + (yres/2)^2));

% Now compute the radial distance of the centre of each pixel from the centre of the image
r = max(y) * abs(repmat(1:1:xres, yres, 1) - xres/2 + j * (repmat((1:1:yres)',1, xres) - yres/2)) / semidiag;

% Finally, interpolate TD and NITD on this grid of radii
imTD = interp1(y, TD, r);
imNITD = interp1(y, NITD, r);
imH = interp1(y, H, r);
imH_rho = interp1(y, H_rho, r);

% Perform some scaling and shifting - currently this is not correct to yield a good representation
% SiTF must be taken into account.
% imH and imH_rho are probably the better indication at this stage.
minim = min(min(imTD));
imTD = imTD - minim;
imTD = imTD/max(max(imTD));

minim = min(min(imNITD));
imNITD = imNITD - minim;
imNITD = imNITD/max(max(imNITD));

% imH = 0.75 * imH/max(max(imH));
% imH_rho = 0.75 * imH_rho/max(max(imH_rho));

% Display the image non-uniformity representations
figure;
imshow(imTD);
title(['Total Image Uniformity for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K']);

figure;
imshow(imNITD);
title(['NITD Uniformity for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K']);

% figure;
% imshow(imH);
% title(['Image Irradiance for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K']);
% xlabel('Note : \rho is neglected');

% figure;
% imshow(imH_rho);
% title(['Image Irradiance for T_a = ' num2str(T_a) 'K    T_d = ' num2str(T_d) 'K    T_h = ' num2str(T_h) 'K']);
% xlabel('Note : \rho is included');

% Save the images for later reference

imwrite(imTD, [NarcDir 'FieldUniformity.jpg']);
imwrite(imNITD, [NarcDir 'NITDUniformity.jpg']);
% imwrite(imH, [NarcDir 'ImageIrradianceWithoutRho.jpg']);
% imwrite(imH_rho, [NarcDir 'ImageIrradianceWithRho.jpg']);

% Set up results to pass back
Results.Revision = Revision;        % Revision of NarcWiz which produced this output
Results.RevDate = RevDate;          % Date of Revision which produced this output
Results.LensFile = LensFilePath;    
Results.LensDir = LensDir;
Results.NarcDir = NarcDir;
Results.Worst = (cat(1, {'Rank', 'Surface Number', 'Contribution', 'YNI'}, num2cell(Worst)))';
Results.WorstSlope = (cat(1, {'Rank', 'Surface Number', 'Slope', 'YNI'}, num2cell(WorstSlope)))';
Results.T_a = T_a;    % Ambient scene temperature
Results.T_d = T_d;    % Internal dewar temperature
Results.T_h = T_h;    % Camera housing temperature
Results.r_c = r_c*100;% Reflectivity of the cold stop
Results.D_c = R_c*2;  % Diameter of the cold stop
Results.D_wo = R_wo*2;% Outer window clear aperture diameter
Results.D_wi = R_wi*2;% Inner window clear aperture diameter
Results.Units = Units;% Lens dimensional units
Results.maxgradNITD = maxgradNITD;      % Maximum gradient of total NITD
Results.maxgradTD = maxgradTD;          % Maximum gradient of total image non-uniformity
Results.maxgradNITD_i = maxgradNITD_i;  % Maximum gradient of surface contributions to NITD