% AnsiZ136DeloriTest
%
% ****************************************************************************
% IMPORTANT: Before using the AnsiZ136 routines, please see the notes on usage
% and responsibility in PsychAnsiZ136MPE/Contents.m (type "help PsychAnsiZ136MPE"
% at the Matlab prompt.
% ****************************************************************************
%
% Tests radiometric calculations and light safety calculations.
%
% Reads in tab delimted text file AnsiZ136MPEDeloriTestInput.txt that has a set of 
% rows in it.  Each row provides stimulus properties for a monochromtic
% light and values obtained from Delori's spreadsheet.
% (There are also some other test values entered by hand in the text file, for
% other radiometric quanties that I could not find computed in Delori's spreadsheet.)
%
% Delori's spreadsheet is described in footnote 49 of Delori et al. (2007,
% JOSA A, 24, pp.1250-1265).  That footnote indicates that the spreadsheet
% will be shared as long as recipients accept full responsibility for its use.
% The version of Delori's spreadsheet used to get these numbers is 0cLC_0.9.9.

% As of 5/1/13 (Notes from DHB):
%   I am unable to get complete agreement between Delori's calculations and
%   the PTB calculations.  But I get agreement for some wavelength, stimulus
%   size, and stimulus duration choices.  To get agreement, this routine
%   uses the effective pupil size formula from the Delori et al. 2007 paper
%   when converting the limit to power in the pupil.
%
%   A key question in going from the number in the 2007 standard and retinal
%   illuminant is what pupil size was assumed in the standard.  Delori et al.
%   (2007) give a formula for this, which (I think) goes from the radiant
%   power in light at the cornea overfilling the pupil to the pupil size. 
%   See Eq. 8.  That is implemented in this test program.
%
%   To get the degree of agreement I obtained, I also had to turn off
%   the limiting cone aperture calculation that is described in Table 2 of the
%   2007 standard.  This is controlled by a flag in the input file.
%
%   Conversions to retinal illuminance done wtth the pupil diameter and
%   eye length specified in the input file.  Delori et al. (2007) indicate
%   that the standard assumes an eye length of 17 mm.  The calculation uses
%   this to compute the allowable retinal illuminance, but if you think some
%   other length is better, you can fill it in.
%
%   In the cases where I don't get agreement with the spreadsheet, the discrepancy
%   arises in the value used for constant T2.  Delori's spreadsheet seems to use
%   10000 for a 2 degree stimulus, but as I read the standard (Table 6, p. 76), this
%   value has a maximum of 100 and a smaller value (21.83) for T2.  This matters when
%   the thermal limit is controlling, for some stimulus durations.  Currently this
%   happens for test conditions 3 and 6.  The code here
%   reproduces what I get when I do the computations by hand.  See files in this
%   directory
%     MPEByHand_580_2deg_100sec.doc
%     MPEByHand_700_2deg_20000sec.doc
%   I have emailed Delori (on 5/2/13) to ask him about this.  And also about the
%   limiting cone aperature calculation mentioned above.
%
% 2/27/13  dhb  Wrote it.

%% Clear and close
clear; close all

%% Load conditions
conditionStructs = ReadStructsFromText('AnsiZ136MPEDeloriTestInput.txt');

%% Loop over conditions and report
for i = 1:length(conditionStructs)
    % Set variables for this condition
    retinalIlluminanceUWattsPerCm2 = conditionStructs(i).retinalIlluminanceUWattsPerCm2;
    wavelengthNm = conditionStructs(i).wavelengthNm;
    stimulusDiameterDegrees = conditionStructs(i).stimulusDiamDegrees;
    stimulusDurationSeconds = conditionStructs(i).stimulusDurationSecs;
    eyeLengthMm = conditionStructs(i).eyeLengthMm;
    pupilDiameterMm = conditionStructs(i).pupilDiameterMm;
    ansiEyeLengthMm = conditionStructs(i).ansiEyeLengthMm;
    ansiPupilDiameterMm = conditionStructs(i).ansiPupilDiameterMm; 
    CONELIMITFLAG = conditionStructs(i).CONELIMITFLAG;
    fprintf('**********\nCondition %d\n\tInput retinal illuminance of %0.1f uWatts/cm2\n',i,retinalIlluminanceUWattsPerCm2);
    fprintf('\t\tWavelength %d nm\n',wavelengthNm);
    fprintf('\t\tStimulus diamter %0.1f degrees\n',stimulusDiameterDegrees);
    fprintf('\t\tStimulus duration %0.1f seconds\n',stimulusDurationSeconds);
    fprintf('\t\tEye length %0.1f mm\n',eyeLengthMm);
    fprintf('\t\tPupil diameter %0.1f mm\n',pupilDiameterMm);
    fprintf('\t\tAssuming ANSI standard eye length of %0.1f mm\n',ansiEyeLengthMm);
    fprintf('\t\tAssuming ANSI pupil diameter of %0.1f mm\n',ansiPupilDiameterMm);
    if (CONELIMITFLAG)
        fprintf('\t\tIncluding limiting cone angle calculation\n');
    else
        fprintf('\t\tExcluding limiting cone angle calculation\n');
    end
    
    % Get comparison values from Delori spreadsheet.  These need to be computed
    % by hand using the spreadsheet and then entered into the condition file.
    % We enter -1 in the spreadsheet for values not computed.
    deloriRadianceMWattsPerCm2Sr = conditionStructs(i).deloriRadianceMWattsPerCm2Sr;
    deloriCornealIrradianceUWattsPerCm2 = conditionStructs(i).deloriCornealIrradianceUWattsPerCm2;
    deloriPowerInPupilMWatts = conditionStructs(i).deloriPowerInPupilMWatts;
    deloriLog10PhotopicTrolands = conditionStructs(i).deloriLog10PhotopicTrolands;
    deloriLog10ScotopicTrolands = conditionStructs(i).deloriLog10ScotopicTrolands;
    deloriMPECb = conditionStructs(i).deloriMPECb;
    deloriMPECe = conditionStructs(i).deloriMPECe;
    deloriMPECt = conditionStructs(i).deloriMPECt;
    deloriMPEPupilFactor = conditionStructs(i).deloriMPEPupilFactor;
    deloriEffectivePupilDiameterMm = conditionStructs(i).deloriEffectivePupilDiameterMm;
    deloriMPETgamma = conditionStructs(i).deloriMPETgamma;
    deloriMPET2 = conditionStructs(i).deloriMPET2;
    deloriMPERetinalIrradianceWattsPerCm2 = conditionStructs(i).deloriMPERetinalIrradianceWattsPerCm2;
    deloriMPEPowerInPupilMWatts = conditionStructs(i).deloriMPEPowerInPupilMWatts;
    
    % In some cases we have other comparison values for other sources.
    % We enter -1 in the spreadsheet for values not computed.
    checkPhotopicLuminanceCdM2 = conditionStructs(i).checkPhotopicLuminanceCdM2;
    checkRetinalIlluminanceQuantaPerCm2Sec = conditionStructs(i).checkRetinalIlluminanceQuantaPerCm2Sec;
    
    % Do unit conversions and print with comparisons when such are available
    retinalIlluminanceWattsPerCm2 = 1e-6*retinalIlluminanceUWattsPerCm2;
    retinalIlluminanceWattsPerUm2 = 1e-8*retinalIlluminanceWattsPerCm2;
    retinalIlluminanceQuantaPerCm2Sec = EnergyToQuanta(wavelengthNm,retinalIlluminanceWattsPerCm2);
    photopicTrolands = RetIrradianceToTrolands(retinalIlluminanceWattsPerUm2, wavelengthNm, 'Photopic','Human',num2str(eyeLengthMm));
    scotopicTrolands = RetIrradianceToTrolands(retinalIlluminanceWattsPerUm2, wavelengthNm, 'Scotopic','Human',num2str(eyeLengthMm));
    photopicLuminanceCdM2 = TrolandsToLum(photopicTrolands,(pi/4)*pupilDiameterMm^2);
    radianceWattsPerM2Sr = RetIrradianceToRadiance(retinalIlluminanceWattsPerUm2,wavelengthNm,(pi/4)*pupilDiameterMm^2,eyeLengthMm);
    radianceMWattsPerCm2Sr = 1e3*1e-4*radianceWattsPerM2Sr;
    cornealIrradianceMWattsPerCm2 = RadianceAndDegrees2ToCornIrradiance(radianceMWattsPerCm2Sr,(pi/4)*stimulusDiameterDegrees^2);
    cornealIrradianceUWattsPerCm2 = 1e3*cornealIrradianceMWattsPerCm2;
    powerInPupilUWatts = cornealIrradianceUWattsPerCm2*(1e-2)*(pi/4)*pupilDiameterMm^2;
    powerInPupilMWatts = 1e-3*powerInPupilUWatts;
    pupilEnergyMJoules = powerInPupilMWatts*stimulusDurationSeconds;
    
    % Print out results, and comparisons if we have them
    AnsiZ136MPEPrintConditionalComparison('\tConverts to retinal illuminance of %0.1f log10 quanta/[cm2-sec]','%0.1f',retinalIlluminanceQuantaPerCm2Sec,checkRetinalIlluminanceQuantaPerCm2Sec,true);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to %0.2f log10 photopic trolands','%0.2f',photopicTrolands,10^deloriLog10PhotopicTrolands,true);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to %0.2f log10 scotopic trolands','%0.2f',scotopicTrolands,10^deloriLog10ScotopicTrolands,true);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to %0.1f cd/m2','%0.1f',photopicLuminanceCdM2,checkPhotopicLuminanceCdM2,false);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to radiance %0.1f mWatts/[cm2-sr]','%0.1f',radianceMWattsPerCm2Sr,deloriRadianceMWattsPerCm2Sr,false);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to corneal irradiance %0.1f uWatts/cm2','%0.1f',cornealIrradianceUWattsPerCm2,deloriCornealIrradianceUWattsPerCm2,false);
    AnsiZ136MPEPrintConditionalComparison('\tConverts to total radiant power in the pupil of %0.2g mW','%0.2g',powerInPupilMWatts,deloriPowerInPupilMWatts,false);
        
    % Compute MPE, with comparisons when available
    [MPELimitIntegratedRadiance_JoulesPerCm2Sr, ...
        MPELimitRadiance_WattsPerCm2Sr, ...
        MPELimitCornealIrradiance_WattsPerCm2, ...
        MPELimitCornealRadiantExposure_JoulesPerCm2] = ...
        AnsiZ136MPEComputeExtendedSourceLimit(stimulusDurationSeconds,stimulusDiameterDegrees,wavelengthNm,CONELIMITFLAG);
    MPELimitRadiance_WattsPerM2Sr =  1e4*MPELimitRadiance_WattsPerCm2Sr;
    
    % For comparison with Delori's calculations, we need to convert light at the cornea to retinal illuminance.  That in turn
    % requires an assumption about pupil size for whatever light is being measured.  We can use equation 8 of Delori et al. (2007, JOSA)
    % to match the assumption in Delori's calculations.
    MPEPupilFactor = AnsiZ136MPEComputePupilFactor(stimulusDurationSeconds,wavelengthNm);
    effectivePupilAreaCm2 = (pi/4)*((0.7)^2)/(MPEPupilFactor);
    effectivePupilDiameterMm = 10*sqrt(effectivePupilAreaCm2/(pi/4));
    if (abs((1e-2*(pi/4)*(effectivePupilDiameterMm^2))-effectivePupilAreaCm2) > 1e-7)
        error('Algegra boo-boo');
    end
    MPELimitPowerInPupilWatts = MPELimitCornealIrradiance_WattsPerCm2*effectivePupilAreaCm2;
    MPELimitPowerInPupilMWatts = 1e3*MPELimitPowerInPupilWatts;
    
    MPELimitRetinalIlluminanceWattsPerUm2 = RadianceToRetIrradiance(MPELimitRadiance_WattsPerM2Sr,wavelengthNm,(pi/4)*effectivePupilDiameterMm^2,ansiEyeLengthMm);
    MPELimitRetinalIlluminanceWattsPerCm2 = 1e8*MPELimitRetinalIlluminanceWattsPerUm2;
    MPELimitRetinalIlluminanceUWattsPerCm2 = 1e6*MPELimitRetinalIlluminanceWattsPerCm2;
    
    % We'll also compute some of the the factors that go into the MPE computation, for comparison with
    % what Delori's spreadsheet reports for the same numbers
    MPECb = AnsiZ136MPEComputeCb(wavelengthNm);
    MPECe = AnsiZ136MPEComputeCe(stimulusDiameterDegrees);
    MPRT2 = AnsiZ136MPEComputeT2(stimulusDiameterDegrees);
    MPECa = AnsiZ136MPEComputeCa(wavelengthNm);
    if (wavelengthNm >= 1050)
        MPECc = AnsiZ136MPEComputeCc(wavelengthNm);
    else
        MPECc = 1;
    end
    MPECt = MPECa*MPECa;
    MPET2 = AnsiZ136MPEComputeT2(stimulusDiameterDegrees);

    % Print out comparisons when available
    fprintf('\nMPE calculations\n');
    AnsiZ136MPEPrintConditionalComparison('\tUsing pupil factor %0.2f','%0.2f',MPEPupilFactor,deloriMPEPupilFactor,false);
    AnsiZ136MPEPrintConditionalComparison('\tEffective pupil diameter is %0.1f mm','%0.1f',effectivePupilDiameterMm,deloriEffectivePupilDiameterMm,false);
    AnsiZ136MPEPrintConditionalComparison('\tUsing Cb %0.2f','%0.2f',MPECb,deloriMPECb,false);
    AnsiZ136MPEPrintConditionalComparison('\tUsing Ce %0.2f','%0.2f',MPECe,deloriMPECe,false);
    AnsiZ136MPEPrintConditionalComparison('\tUsing Ct %0.2f','%0.2f',MPECt,deloriMPECt,false);
    AnsiZ136MPEPrintConditionalComparison('\tUsing T2 %0.2f','%0.2f',MPET2,deloriMPET2,false);
    AnsiZ136MPEPrintConditionalComparison('\tMPE power in pupil limit %0.3g mWatts','%0.3g',MPELimitPowerInPupilMWatts,deloriMPEPowerInPupilMWatts,false);
    AnsiZ136MPEPrintConditionalComparison('\tMPE retinal illuminance limit computed as %0.3g Watts/cm2','%0.3g',MPELimitRetinalIlluminanceWattsPerCm2,deloriMPERetinalIrradianceWattsPerCm2,false);
    fprintf('\tLimit - Stimulus log10 difference: %0.1f log10 units\n',log10(MPELimitRetinalIlluminanceWattsPerCm2)-log10(retinalIlluminanceWattsPerCm2));

    % Ready for next iteration
    fprintf('\n');

end

%% See if we can match some conversions computed
% by Ed Pugh and conveyed to me by Brian Wandell.
%
% Ed starts with retinal illuminance of 10^15 in quanta/[cm2-sec]'
%   He gets 340 uW/cm2.  We get 342.
%   He gets ~590,000 trolands.  We get about this assuming 17 mm eye length.
%   He gets ~190,000 cd/m2 for a 2mm diameter pupil.  We also get about this.




%% Check ANSI light limit calculations against numbers in Eds document.  He doesn't say
% the durations or size he assumed, but he does say he got the numbers from Delori's
% spreadsheet.  Let's try making up a source size and duration and see what happens.


% When I plug these numbers (580 nm, 2 degree stimulus, 2 mm pupil diameter, 1 second exposure)
% into the version of Delori's spreadsheet I got via Ed Pugh (rev 1/10/08), it computes that
% the stimulus has:
%   a retinal irradiance of 340 uW/cm2 
%   a radiance of 31.28 mW/[cm2 sr]
%   a corneal irradiance of 29.93 uW/cm2,
%   total radiant energy in the pupil of 940.33 nJ.
%   retinal radiant exposure of 340 uJ/cm2.
% We compute these quantities using PTB routines.  From above we already have the retinal irradiance
% at about 340 Watts/cm2. The rest of the numbers also match up pretty well.  Rounding in the spreadsheet
% plus perhaps a different assumption about eye length can explain the differences, I think.

% The spreadsheet computes the exposure safety limit for this stimulus as:
%  radiant power in the pupil 2.96 mW 
%  radiant energy in the pupil of 2.96 mJ
%  retinal irradiance 1.07 W/cm2
%  retinal exposure of 1.07 J/cm2