Debugging the LED analysis for human safety

color/ieLuminance2Radiance.m

@@ -1,65 +1,78 @@
-function [energy, photons] = ieLuminance2Radiance(lum,wave,varargin)
-% Convert the luminance of a monochromatic light to radiance (energy or
-% photons)
+function [energy, wave] = ieLuminance2Radiance(lum,thisWave,varargin)
+% From a monochromatic LED luminance, compute the spectral radiance (energy)
 %
 % Synopsis
-%   [energy, photons] = ieLuminance2Radiance(lum,wave)
+%   [energy, photons] = ieLuminance2Radiance(lum,thisWave,varargin)
 % 
 % Inputs
-%   lum:   Luminance of the light in cd/m2
-%   wave:  Wavelength of the monochromatic light
+%   lum:       Luminance of the monochromatic light in cd/m2
+%   thisWave:  Wavelength of the monochromatic light (between 350 and 720)
 %
 % Optional key/val pairs
-%   bin width:  How wide is the wavelength bin width (default is 10 nm)
+%   sd:  Standard deviation of the Gaussian spread of the spectral radiance
+%        (default is 10nm)
 %
 % Returns
-%   energy:   watts/sr/nm/m2
-%   photons:  photons/sr/nm/m2
+%   energy:   watts/sr/nm/m2  (a vector)
+%   wave:     wavelength samples (nm) of the energy vector
 %
 % Description
-%   The radiance of a monochromatic light (watts/sr/nm/m2) is scaled by the
-%   V(\lambda) function and some numerical constants to return the
-%   luminance in cd/m2.  This function inverts the scaling to return
-%   the radiance.
+%   Many LEDs have a spectral radiance close to a Gaussian centered at some
+%   wavelength. If we measure only the luminance of the LED, we can model
+%   the radiance energy (watts/sr/nm/m2) as a Gaussian centered at that
+%   wavelength.  We then scale the spectral radiance energy so that it has
+%   the desired luminance (cd/m2).
 %  
-%   The energy can be converted into photons as well.
-%
 % See also
-%   s_humanSafety, ieLuminanceFromEnergy, ieLuminanceFromPhotons
+%   s_humanSafety, ieLuminanceFromEnergy, ieLuminanceFromPhotons,
+%   Energy2Quanta
 
 % Examples:
 %{
-   energy = ieLuminance2Radiance(100,405);
-%}
-%{
-   lum = 100; wave = 405;
-   energy = ieLuminance2Radiance(100,wave,'bin width',10);
-   assert(lum == ieLuminanceFromEnergy(energy,wave))
+   % These are unit tests.  Make a spectral radiance (energy)
+   % with the peak wavelengty and Gaussian spread, as one might find for an
+   % LED.  Then check that the energy of the LED model has a luminance that
+   % equals the desired luminance. 
+
+   lum = 10; thisWave = 360;
+   [energy,wave] = ieLuminance2Radiance(lum,thisWave);
+   assert(lum - ieLuminanceFromEnergy(energy,wave) < 1e-10)
+   plotRadiance(wave,energy);
 %}
 %{
-
+   lum = 100; thisWave = 425;
+   energy = ieLuminance2Radiance(lum,thisWave,'sd',20);
+   assert(lum - ieLuminanceFromEnergy(energy,wave) < 1e-10)
+   plotRadiance(wave,energy);
 %}
 
-%% 
+%% Check input parameters
 p = inputParser;
 varargin = ieParamFormat(varargin);
 p.addRequired('lum',@isnumeric);
-p.addRequired('wave',@isnumeric);
-p.addParameter('binwidth',10,@isnumeric);
+p.addRequired('thisWave',@(x)(x >= 350 && x <= 720));
+p.addParameter('sd',10,@isnumeric);   % Standard deviation of the Gaussian radiance model
 
-p.parse(lum,wave,varargin{:});
-binwidth = p.Results.binwidth;
+p.parse(lum,thisWave,varargin{:});
+sd = p.Results.sd;
 
-%% For this wave, calculate the scale factor from radiance to luminance
+%% Model the radiance as a Gaussian centered at thisWave
 
-sFactor = ieLuminanceFromEnergy(1,wave);
+wave   = 300:1:770;   % I chose a big range because we used for UV LEDs
+energy = zeros(numel(wave),1);
+energy(wave == thisWave) = 1;   % Set an impulse at the center wave
 
-energy = lum/sFactor * (binwidth/10);
+% Make the Gaussian and convolve the impulse
+g = fspecial('gaussian',[8*sd,1],sd);
+energy = conv(energy,g,'same');
 
-if nargout == 2
-    sFactor = ieLuminanceFromPhotons(1,wave);
-    photons = lum/sFactor * (binwidth/10);
-end
+% plotRadiance(wave,energy);
+
+% Figure out the luminance of the LED model radiance
+sFactor = ieLuminanceFromEnergy(energy,wave);
+
+% Scale the model radiance to the right luminance
+energy = energy * (lum/sFactor);
 
 end
 

color/ieLuminanceFromEnergy.m

@@ -77,6 +77,7 @@
 
 fName = fullfile(isetRootPath,'data','human','luminosity');
 V = ieReadSpectra(fName,wave);
+% ieNewGraphWin; plot(wave,V);  % The luminance curve
 
 % 683 is the standard factor for conversion when the energy are in Watts.
 % The wavelength difference accounts for the wavelength sampling.
@@ -85,8 +86,11 @@
 end
 
 % The luminance formula.  xwData and V are not always rows or columns, so
-% we use the dot() formula here rather than multiplying.
+% we use the dot() formula  rather than multiplying.
 lum = 683*dot(xwData,V) * binwidth;
 
+% Compare the luminance and energy data
+% ieNewGraphWin; semilogy(wave,V,'--',wave,xwData,'o'); 
+
 end
 

human/humanUVSafety.m

@@ -19,15 +19,27 @@
 %  The data for the Actinic safety function curves were taken from this
 %  paper
 %
-%  IEC 62471:2006 Photobiological Safety of Lamps and Lamp Systems. n.d.
-%  Accessed October 5, 2019. https://webstore.iec.ch/publication/7076 
-%  J.E. Farrell has a copy of this standard
+%    IEC 62471:2006 Photobiological Safety of Lamps and Lamp Systems. n.d.
+%    Accessed October 5, 2019. https://webstore.iec.ch/publication/7076 
+%    J.E. Farrell has a copy of this standard
 %
 % Notes:   Near UV is also called UV-A and is 315-400nm.
 %
 % See also
 %   s_humanSafetyUVExposure
 
+%%  Near-UV hazard exposure limit
+%{
+This calculation has no weighting function.  
+
+For times less than 1000 sec, add up the total irradiance
+from 315-400 without any hazard function (Equation 4.3a).  Call this E_UVA.
+Multiply by time (seconds).  The product should be less than 10,000.
+
+For times exceeding 1000 sec, add up the total irradiance, divide by the
+time, and the value must be less than 10 (Equation 4.3b).
+
+%}
 
 %% Check inputs
 
@@ -40,6 +52,7 @@
 
 fname = which('Actinic.mat');
 Actinic = ieReadSpectra(fname,wave);
+% semilogy(wave,Actinic); xlabel('Wave'); grid on; 
 
 %% The UV hazard safety formula
 %   

scripts/human/s_humanSafetyUVExposure.m

@@ -90,24 +90,13 @@
 The person multiplies b6000 by pi, not 2pi
 %}
 
-%%  Create a radiance and convert it to irradiance
-
-wave       = 300:700;
-radiance   = 10*blackbody(wave,9000);
-irradiance = radiance*pi;
-
-% ieNewGraphWin; plot(wave,irradiance);
-
-%% Run the UV hazard safety function
-
-exposureMinutes = humanUVSafety(irradiance,wave);
-fprintf('Safe exposure (hours) for 8 hour period is %.2f minutes\n',exposureMinutes);
-
 %% An example of a light measured in the lab
 
+wave = 300:770;
 fname = which('LED405nm.mat');
 radiance = ieReadSpectra(fname,wave);
 radiance = mean(radiance,2);
+lum405 = ieLuminanceFromEnergy(radiance,wave);
 plotRadiance(wave,radiance);
 irradiance = pi*radiance;
 
@@ -116,63 +105,48 @@
 
 %% An example of the 385nm light in the OralEye camera
 
-%fname = which('LED385nm.mat');
-
-fname = which('LED405nm.mat');
+fname = which('LED385nm.mat');
 wave = 300:700;
 radiance = ieReadSpectra(fname,wave);
 radiance = mean(radiance,2);
 plotRadiance(wave,radiance);
 
 irradiance = pi*radiance;
 exposureMinutes = humanUVSafety(irradiance,wave);
-
 fprintf('Maximum exposure duration per eight hours:  %f (min)\n',exposureMinutes)
 
 % Each exposure is very brief.  So you can safelty take this many exposures
-fprintf('For a 30 ms exposure, you can take %d exposures\n',round((exposureMinutes*60)/0.030));
+fprintf('For a 30 ms exposure, you can take %d exposures in an eight hour period.\n',round((exposureMinutes*60)/0.030));
 
-lum405 = ieLuminanceFromEnergy(radiance,wave);
-radiance2 = ieLuminance2Radiance(lum405,405,'bin width',dLambda/6); % watts/sr/nm/m2
-lst = (wave == 405);
+%%  The mean daylight we measured in California
 
-%% Start with a monochromatic light luminance
+wave       = 300:700;
+[radiance,wave] = ieReadSpectra('DaylightPsychBldg.mat',wave);
+plotRadiance(wave,radiance);
 
+% Convert radiance to irradiance
+irradiance = radiance*pi;
 
-% Suppose we know the luminance of a 380 nm light with a 10 nm bin width
-% lum = 31; wave = 385; dLambda = 3;
+exposureMinutes = humanUVSafety(irradiance,wave);
+fprintf('Safe exposure (hours) for 8 hour period is %.2f minutes (%.2f hours)\n',exposureMinutes,exposureMinutes/60);
 
-% We convert the luminance to energy
-dLambda = 4;
-radiance2 = ieLuminance2Radiance(lum405,405,'bin width',dLambda); % watts/sr/nm/m2
+%% If you only know the luminance of an LED (monochromatic) and its bandwidth (s.d.)
 
-irradiance = pi*radiance;
+lum = lum405;   % cd/m2, luminance of the 405 LED
+thisWave = 405; % nm, center wavelength of the LED
+bandwidth = 15; % nm, Gaussian standard deviation, FW at roughly 1/2 of the max
+[estRadiance,wave] = ieLuminance2Radiance(lum,thisWave,'sd',bandwidth); 
+plot(wave,radiance,'--',wave,estRadiance,'o');
+legend({'meas rad','est rad'})
+
+ieLuminanceFromEnergy(estRadiance,wave)
+ieLuminanceFromEnergy(radiance,wave)
 
+irradiance = pi*estRadiance;
 exposureMinutes = humanUVSafety(irradiance,wave);
 
 % Conver the hazard energy into maximum daily allowable exposure in minutes
 fprintf('Maximum exposure duration per eight hours:  %f (min)\n',exposureMinutes)
+fprintf('For a 30 ms exposure, you can take %d exposures in an eight hour period.\n',round((exposureMinutes*60)/0.030));
 
-
-%%  Plot the Actinic hazard function 
-
-ieNewGraphWin;
-mx = max(irradiance(:));
-dummy = ieScale(Actinic,1)*mx;
-plot(wave,irradiance,'k-',wave,dummy,'r--','linewidth',2);
-xlabel('Wave (nm)');
-ylabel('Irradiance (watts/m^2');
-grid on
-legend({'Irradiance','Normalized hazard'});
-
-%%  Near-UV hazard exposure limit
-%{
-This calculation has no weighting function.  
-
-For times less than 1000 sec, add up the total irradiance
-from 315-400 without any hazard function (Equation 4.3a).  Call this E_UVA.
-Multiply by time (seconds).  The product should be less than 10,000.
-
-For times exceeding 1000 sec, add up the total irradiance, divide by the
-time, and the value must be less than 10 (Equation 4.3b).
-%}
+%%