function Results = NarcWiz
% NarcWiz - Narcissus analysis wizard
%
% Usage : NarcOut = NarcWiz      % NarcOut is case sensitive
%         report narcissus
%
% Takes the user through a narcissus analysis of the lens in the ZEMAX DDE server. NarcOut is a
% structure containing various items required for the report generator, so the function must be
% called exactly as shown above. The report will be generated in rich text format (.rtf) and
% displayed in MS Word.
%
% The general restrictions and assumptions are as follows :
% 1) The lens model to be analysed is for a 3rd generation thermal imager using a cooled detector
%    which has a cold stop. The cold stop is assumed to be at the surface before the image surface
%    and the detector window is assumed to come before the cold stop. Clear apertures of these
%    objects must be set correctly.
% 2) The inside of the dewar has emissivity of 1 all over and across the spectral band.
% 3) The side of the cold stop facing outward has a user-defined emissivity and the same temperature
%    as the inside of the dewar right out to the clear aperture of the dewar window.
% 4) Emissions and reflections of thermal radiation by non-optical surfaces are lambertian.
% 5) Lens must be single configuration. Save and analyse each configuration seperately, taking
%    special care to ensure that the apertures are as they would be for the system as-built. It is
%    particularly important that the clear aperture of the dewar window is correctly set from
%    mechanical drawings of the dewar.
% 6) Only rotationally symmetrical lenses are currently handled - mainly a
%    retriction of the ZEMAX functions used for the computation. The whole
%    approach to the problem would have to change significantly to cater
%    for non-symmetrical systems.
%
% The general functions performed are as follows :
% 1) Create a directory named after the lens to be analysed in which all results will be saved. The
%    directory name will have the .Narcissus extension.
% 2) Ensure that the system is correctly set up, has clear apertures and coatings defined. Currently this is
%    done by appealing to the users's sense of thoroughness.
% 3) Compute system transmission, YNI contributions and uniformity of image illumination.
% 4) Reverse the lens and create single bounce ghost reflection models for each optical surface.
% 5) Compute the Narcissus-Induced Delta T contributions for each surface.
% 6) Plot total NITD, NITD per surface and overall image irradiance uniformity per configuration.
% 7) Generate subjective images of image non-uniformity.
% 8) Rank surfaces in order of contribution to NITD and tabulate with YNI contributions.
% 9) Compute maximum slope of total NITD and overall image non-uniformity. This gives an indication
%    of the local visibility of the image non-uniformities.
% 10) If the user gives the command 'report narcissus', a report will be generated in MS Word,
%     containing all of the outputs.
%
% Possible return values for NarcWiz are
%  0 - Everything seems to be OK
% -1 - ZEMAX not running and could not start ZEMAX.
% -2 - The lens you wanted to analyse cannot be raytraced.
% -3 - Lens has fewer than 6 surfaces. Detector + cold stop + window + 2 additional lens surfaces is
%      six at least.
% -4 - Lens does not operate primarily in a thermal band.
% -5 - Lens is multi-config. Currently we only do single config. Save and analyse each config
%      seperately.
% -6 - User cancelled wizard.
% -7 - Double bounce ghost lenses were encountered - restart analysis.
% -8 - Lens has wrong field type. Change to "real image height".

% MZDDE - The ZEMAX DDE Toolbox for Matlab.
% Copyright (C) 2002-2004 Defencetek, CSIR
% Contact : dgriffith@csir.co.za
% 
% This file is part of MZDDE.
% 
%  MZDDE is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  MZDDE is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with MZDDE (COPYING.html); if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%



% $Revision: 1.4 $
% $Author: dgriffith $

Revision = '$Revision: 1.4 $';
Revision = Revision(11:(size(Revision,2)-1));
disp(['Narcissus Analysis Wizard. Version' Revision]);
RevDate = '$Date: 2005-04-22 09:35:38+02 $';
RevDate = RevDate(7:(size(RevDate,2)-1));
disp(['Copyright Defencetek, CSIR,' RevDate]);

NarcIcon=1:64; NarcIcon=(NarcIcon'*NarcIcon)/64; % Make an icon for use with dialogs

Results.Status = zDDEInit; % Attempt to contact ZEMAX through DDE.
if (Results.Status == -1)        % Failed, try to start ZEMAX and ask user to select the file
    Results.Status = zDDEStart;
    if (Results.Status == 0)     % Success, now do the open file dialog
        uiwait(msgbox('ZEMAX is Starting. Please accept DDE commands from Matlab and open the lens file for Narcissus Analysis.', 'ZEMAX Startup', 'custom', NarcIcon, hot(64),'modal'))
        zOpenWindow('Ope'); zWindowMaximize(0);
    end
end
if (Results.Status == -1) return; end; % Give up

while (zGetRefresh == -2); end;       % Wait for user to accept DDE messaging and open the lens file

if (zGetRefresh ~= 0) 
    uiwait(msgbox('The lens you are attempting to analyse cannot be raytraced.', 'Fatal Lens Error', 'error', 'modal'));
    Results.Status = -2;
    return; % Give up
end
zPushLens(10);
while (zGetRefresh == -2); end;       % Wait for user to accept DDE messaging and open the lens file

% Get some lens data, including the wavelengths and the number of lens surfaces.
LenSys = zsGetSystem;
if (LenSys.numsurfs < 6)
    uiwait(msgbox('The lens you are attempting to analyse has less than 6 Surfaces.', 'Fatal Lens Error', 'error', 'modal'));
    Results.Status = -3;
    return;
end

% Check that the field type is set to "real image height"
FieldData = zGetField(0); FieldType = FieldData(1);
if (FieldType ~= 3)
    uiwait(msgbox('The lens you are attempting to analyse has the wrong field type. Change field type to "Real Image Height"', 'Fatal Lens Error', 'error', 'modal'));
    Results.Status = -8;
    return;
end;

PrimaryWave = zGetWave(0); NumberWaves = PrimaryWave(2); PrimaryWave = PrimaryWave(1);
PrimaryWaveLength = zGetWave(PrimaryWave); % Get the primary wavelength

if (PrimaryWaveLength(1) < 1) | (PrimaryWaveLength(1) > 14)
    uiwait(msgbox('The primary wavelength for this lens seems to be out of the normal thermal bands.', 'Fatal Lens Error', 'error', 'modal'));
    Results.Status = -4;
    return;
end

% By now we think there is a useable lens in the ZEMAX DDE server, and we start with analysis.

% Firstly find the location of the lens and create a directory for Narcissus analysis output
% The directory has the same name as the lens file, but has the extension .Narcissus
LensFilePath = zGetFile; 
[LensDir, LensFile] = fileparts(LensFilePath);
[Success, Message, MessageID] = mkdir(LensDir, [LensFile '.Narcissus']);
NarcDir = [LensDir '\' LensFile '.Narcissus\'];

% If directory already exists
if (strcmp(MessageID, 'MATLAB:MKDIR:DirectoryExists'))
     % uiwait(msgbox('The .Narcissus directory for this lens already exists. Old Data will be overwritten.', 'Directory Warning', 'custom', NarcIcon, hot(64), 'modal'));
    if (exist([NarcDir 'rho.txt'], 'file'))
      Answer = questdlg('Relative Illumination Data for this lens already exists. Recompute Relative Illumination ?', 'Recalculate Relative Illumination');
      if (strcmp(Answer,'Yes')), delete([NarcDir 'rho.txt']); end; if (strcmp(Answer, 'Cancel')), Results.Status = -6; return; end;
    end;
    if (exist([NarcDir 'nu.txt'], 'file'))
      Answer = questdlg('Vignetting Data for this lens already exists. Recompute Vignetting Data ?', 'Recalculate Veignetting Data');
      if (strcmp(Answer,'Yes')), delete([NarcDir 'nu*.txt']); delete([NarcDir 'cnu*.txt']); end;  if (strcmp(Answer, 'Cancel')), Results.Status = -6; return; end;  
    end;
    delete([NarcDir 'GH*.ZMX']);          % Delete all ghost lenses in this directory
end

% Here we will deal with multiconfiguration lenses. For now just dump the user if the lens is
% multi-config.

ConfigData = zGetConfig;
if (ConfigData(2) > 1)
    uiwait(msgbox('This wizard currently only handles single configuration systems.', 'Fatal Lens Error', 'error', 'modal'));
    Results.Status = -5;
    return;
end

uiwait(msgbox('Please ensure that all lens surfaces have correct clear apertures and coatings. Only then click OK in this message box.', 'Apertures/Coatings', 'custom', NarcIcon, hot(64), 'modal'));

% Get the latest lens
zGetRefresh;
% and save it as 'Primary.zmx' in the analysis directory
zSaveFile([NarcDir 'Primary.zmx']);

% Get some global variables, including the ambient scene temperature, the dewar termperature, the
% housing temperature, and the cold stop reflectivity.
Prompt = {'Scene Temperature (K)', 'Dewar Temperature (K)', 'Camera Housing Temperature (K)', 'Cold Stop Reflectivity (%)'};
DefAns = {'300',                   '120',                   '300',                            '1'                         };
Answer = inputdlg(Prompt, 'Narcissus Wizard : Temperature and Emmissivity Inputs', 1, DefAns);
if (~isempty(Answer))
  T_a = str2double(Answer{1});        % T_a = Scene ambient temperature
  T_d = str2double(Answer{2});        % T_d = Internal Dewar temperature in kelvins
  T_h = str2double(Answer{3});        % T_h = optical housing temperature
  r_c = str2double(Answer{4})/100.0;  % r_c = reflectivity of the outer surface of the cold stop
else
    Results.Status = -6;
    return;
end

% Compute the effective F/number from the size of the cold stop and it's distance from the focal
% plane array.
ImSurf = LenSys.numsurfs;
ColdStopSurf = ImSurf - 1;
SurfDataCode = zSurfDataCodes;
ColdStopSemiDia = zGetSurfaceData(ColdStopSurf, SurfDataCode.semidia);
ColdStopDist = zGetSurfaceData(ColdStopSurf, SurfDataCode.thickness);
F = 1 / (2 * sin(atan(ColdStopSemiDia/ColdStopDist)));

% Now compute axial transmission, image illumination uniformity and YNI contributions.
% To compute axial transmission, trace an axial ray.
PolTraceData = zGetPolTrace(PrimaryWave, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0);
tau_S = PolTraceData(2); % This is axial system transmission. An enhanced version could do better.

% Compute YNI contributions and image uniformity using zGetTextFile
% First the settings file must be located
% Settings for the vignetting and relative illumination data are particularly important
mDir = [fileparts(which('NarcWiz')) '\'];
SettingsFile = [LensDir '\' LensFile '.cfg']; % Get the settings from the originating file
zGetTextFile([NarcDir 'YNIf.txt'], 'Yni', SettingsFile, 1); % Compute and save YNI contributions
BusyMessage = ZEMAXBusy; pause(1); % Display a message telling the user that ZEMAX is busy
while (zDDEBusy), end; % Wait for command to complete
if (~exist([NarcDir 'rho.txt']))
  zGetTextFile([NarcDir 'rho.txt'], 'Rel', SettingsFile, 1); % Compute and save relative illumination
end;
while (zDDEBusy), end; % Wait for command to complete
delete(BusyMessage); % Pack away the message

% The next step is to reverse the lens and patch up
zWindowMaximize(0);
zOpenWindow('Rev'); % Reverse elements
uiwait(msgbox('Switch to ZEMAX and select all lens surfaces for reversal. When done click OK.', 'Reverse Warning', 'custom', NarcIcon, hot(64), 'modal'));
zGetRefresh; % Lens should be reversed now
% Fix up surface 0
zSetSurfaceData(0, SurfDataCode.thickness, ColdStopDist);
% Fix up surface 1 (the cold stop)
zDeleteSurface(1);
zSetSurfaceData(1, SurfDataCode.semidia, ColdStopSemiDia);
zSetAperture(1, 1, 0, ColdStopSemiDia, 0, 0, '');
% Set aperture characteristics
zSetSystemAper(0, 1, ColdStopSemiDia*2);
% Set the field type to object height
FieldData = zGetField(0);
zSetFieldType(1, FieldData(2));

% Fixup coatings
Coating = 'I.99';
LenSys = zsGetSystem;
for Surf = 1:(LenSys.numsurfs),
    if (~strcmp(zGetSurfaceData(Surf, SurfDataCode.glass),zGetSurfaceData(Surf-1, SurfDataCode.glass)))
        zSetSurfaceData(Surf, SurfDataCode.coating, Coating);
    else
        zSetSurfaceData(Surf, SurfDataCode.coating, '');
    end
end

% Remove the aperture from the image surface
zSetAperture(LenSys.numsurfs, 0, 0, 0, 0, 0, '');

% Get the radius of the cold stop
R_c = zGetAperture(1); R_c = R_c(3);

% Get the inner and outer window radii
R_wi = zGetAperture(2); R_wi = R_wi(3);
R_wo = zGetAperture(3); R_wo = R_wo(3);

% Get lens units
switch LenSys.unitcode,
    case 0, Units = 'mm';
    case 1, Units = 'cm';
    case 2, Units = 'in';
    case 3, Units = 'm';
end

zPushLens(10);
% Snatch a picture of the lens for the report
zGetMetaFile([NarcDir 'Primary Reversed.emf'], 'Lay', SettingsFile, 1);