diff --git a/human/humanUVSafety.m b/human/humanUVSafety.m index eb7f8ec8..681525bb 100644 --- a/human/humanUVSafety.m +++ b/human/humanUVSafety.m @@ -1,105 +1,259 @@ -function exposureMinutes = humanUVSafety(irradiance,wave) +function [val, level] = humanUVSafety(energy,wave,varargin) % Calculate the UV hazard safety for an irradiance % % Synopsis -% exposureMinutes = humanUVSafety(irradiance,wave) -% +% [val, level] = humanUVSafety(irradiance,wave,varargin) +% % Inputs: -% irradiance - in energy (watts/nm/m2) -% wave - wavelength samples of the irradiance +% energy - irradiance (skineye, eye) or radiance (bluehazard) (watts/nm/m2) +% wave - wavelength samples of the irradiance +% +% Optional key/val pairs +% method - 'skineye' or 'eye' (sections 4.3.1 and 4.3.2 of the PDF) +% duration - Stimulus duration in seconds used for 'eye' and 'bluehazard' +% methods +% +% Return: % -% Returns -% exposureMinutes - Maximum safe exposure time (minutes) per eight hour -% period +% Method: 'skineye' +% val - Maximum safe exposure time (minutes) per eight hour period +% level - hazardEnergy +% +% Method: 'eye' +% val - logical (true/false) +% level - the irradiance integrated over wavelength +% +% Method: 'bluehazard' +% val - logical (true/false) +% level - the dot product of the radiance with the blueLightHazard +% function. % % Description: % +% (Near UV is also called UV-A and is 315-400nm) +% % We caclulate for a constant light (not time-varying). % % 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 +% 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. +% Google Drive: +% https://drive.google.com/file/d/1JQC_R7QDGd58vpNByyKvrVlwEfTk5-Uu/view?usp=sharing +% +% The document explains the different functions in these sections. +% +% Actinic UV hazard exposure limit skin and eye (4.3.1) +% Near UV hazard limit for the eye (4.3.2) +% Retinal blue light hazard exposure (4.3.3) +% Retinal blue light small source (4.3.4) +% Retinal thermal hazard (4.3.5) +% Retinal theermal hazard for weak visual stimulus (4.3.6) +% Infrared exposure for the eye (4.3.7) +% Thermal hazard for the skin (4.3.8) +% +% Notes +% +% I was not sure if radiance -> irradiance is multiplied by 2pi or 1pi. +% So, we checked with Peter and the Internet. The answer: 1pi. +% +% From: Peter B. Catrysse +% Sent: Wednesday, August 21, 2019 11:02 AM +% To: Joyce Eileen Farrell +% Subject: RE: spectroradiometric measurements +% +% Hello Joyce, +% +% It is definitely pi. +% There is 2*pi solid angle for the hemisphere, but when you integrate +% you end up getting pi as factor. +% +% See you next week, +% +% Peter +% +% From: Joyce Eileen Farrell [mailto:jefarrel@stanford.edu] +% Sent: Tuesday, August 20, 2019 10:58 PM +% To: Peter Bert Catrysse +% Subject: Re: spectroradiometric measurements +% +% Hi Peter, +% +% Thanks so very much for giving me the conversion from radiance to +% irradiance. +% +% is it +% E = pi*L/R (where E is irradiance, L is radiance and R is reflectance) +% or +% E=2pi*L/R +% +% ---- +% +% See also this exchange online +% +% https://physics.stackexchange.com/questions/116596/convert-units-for-spectral-irradiance +% +% The person multiplies b6000 by pi, not 2pi +% +% Notes +% +% This calculation is for skin and eye (4.3.1). +% +% We are calculating the irradiance, not the irradiance at the retina. We +% could try to calculate irradiance at the retina, accounting for the pupil +% diameter and the tranmissivity of the lens. But the standard's committee +% baked that into the standard itself. (page 27). +% +% A related issue is the area subtended. In our case we are measuring +% irradiance as watts/m2, so the assumption is that any area is accounted +% for by the units, and that over the area that is illuminated by the light +% the intensity is uniform. % % 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). +% s_humanSafetyUVExposure, ieLuminance2Radiance +% -%} %% Check inputs -if ieNotDefined('irradiance'), error('Spectral irradiance required'); end -if ieNotDefined('wave') || length(wave) ~= length(irradiance) - error('wave representation missing or incorrect'); -end +varargin = ieParamFormat(varargin); -%% Load the standard function +p = inputParser; -fname = which('Actinic.mat'); -Actinic = ieReadSpectra(fname,wave); -% semilogy(wave,Actinic); xlabel('Wave'); grid on; +p.addRequired('irradiance',@isvector); +p.addRequired('wave',@isvector) +p.addParameter('method','skineye',@(x)(ismember(ieParamFormat(x),{'skineye','eye','bluehazard'}))); +p.addParameter('duration',1,@isnumeric); % Used for 'eye' method only +p.parse(energy,wave,varargin{:}); -%% The UV hazard safety formula -% -% sum (Actinic(lambda,t) irradiance(Lambda)) dLambda dTime -% -% Our stimuli are constant over time, so this simplifies to -% -% T * sum(Actinic(lambda) irradiance(lambda) dLambda -% -% where T is the total time. For simplicty we set T = 1, so we are -% calculating the total amount of time (in minutes) for a light with this -% spectral irradiance. -% -% Follow the units this way: -% -% Irradiance units: Watts/m2/nm -% Watts = Joules/sec -% -% Watts/m2/nm * (nm * sec) % Irradiance summed over nm and time -% Joules/sec/m2/nm * (nm * sec) -% Joules/m2 % Becomes Joules/area -% +method = ieParamFormat(p.Results.method); +duration = p.Results.duration; + +%% Set up wavelength sampling -% Summing over wavelength and time (1 sec) if numel(wave) == 1 dLambda = 10; disp('Assuming 10 nm bandwidth'); else dLambda = wave(2) - wave(1); end -hazardEnergy = (dot(Actinic,irradiance) * dLambda); - -% This is the formula from the standard to compute the maximum daily -% allowable exposure from the hazard energy. -%{ - The maximum permissible exposure time per 8 hours for ultraviolet - radiation incident upon the unprotected eye or skin shall be computed by: - t_max = 30/E_s (seconds) (Equation 4.2) +%% Load the standard function +switch method + case 'skineye' + fname = which('Actinic.mat'); + Actinic = ieReadSpectra(fname,wave); + % semilogy(wave,Actinic); xlabel('Wave'); grid on; + + %% The UV hazard safety formula + % + % sum (Actinic(lambda,t) irradiance(Lambda)) dLambda dTime + % + % Our stimuli are constant over time, so this simplifies to + % + % T * sum(Actinic(lambda) irradiance(lambda) dLambda + % + % where T is the total time. For simplicty we set T = 1, so we are + % calculating the total amount of time (in minutes) for a light with this + % spectral irradiance. + % + % Follow the units this way: + % + % Irradiance units: Watts/m2/nm + % Watts = Joules/sec + % + % Watts/m2/nm * (nm * sec) % Irradiance summed over nm and time + % Joules/sec/m2/nm * (nm * sec) + % Joules/m2 % Becomes Joules/area + % + % This standard is defined for 200-400 nm, but we do not have any + % data that go that short. + + % We only evaluate up to a wavelength limit of 400nm + wLimit = (wave <= 400); + hazardEnergy = (dot(Actinic(wLimit),energy(wLimit)) * dLambda); + + % This is the formula from the standard to compute the maximum daily + % allowable exposure from the hazard energy. + %{ + The maximum permissible exposure time per 8 hours for ultraviolet + radiation incident upon the unprotected eye or skin shall be computed by: - E_s is the effective ultraviolet irradiance (W/m^2). The formula for - E_s is defined in Equation 4.1. It is the inner product of the Actinic - function and the irradiance function, accounting for time and - wavelength sampling. -%} + t_max = 30/E_s (seconds) (Equation 4.2) -% We return the time in minutes. -exposureMinutes = (30/hazardEnergy)/60; + E_s is the effective ultraviolet irradiance (W/m^2). The formula for + E_s is defined in Equation 4.1. It is the inner product of the Actinic + function and the irradiance function, accounting for time and + wavelength sampling. + %} + + % We return the time in minutes. + val = (30/hazardEnergy)/60; + level = hazardEnergy; + + case 'eye' + %% The calculation in 4.3.2 is specialized for the eye. + % + % The formula they derive in that case is a maximum irradiance value that + % is calculated separately when the eye is exposed for less than 1,000 sec + % or more than 1,000 seconds. + % + % Suppose for an irradiance level that is constant over time we have + % + % IR = Integral (E(lambda) dlambda) (Watts/m2) + % + % where E(lambda) has units of Watts/nm/m2 + % + % ExpTime * IR has units of (sec * Watts/m2 = Joules/m^2) + % + % Then for times less than 1000 sec, the rule is + % + % ExpTime * IR < 10,000 + % + % For times greater than 1000 sec we require that + % + % IR < 10 + % + wLimit = (wave <= 400); + level = dLambda * sum(energy(wLimit)); + + % val is set to true for safe and false for not safe + val = false; + if duration <= 1000 % seconds + if level*duration < 10000, val = true; end + elseif level < 10 % Duration is long, so level must be low + val = true; + end + case 'bluehazard' + % Retinal blue light hazard, Section 4.3.3 + % + % Potential for a photochemically induced retinal injury + % resulting from radiation exposure at wavelengths primarily + % between 400 nm and 500 nm. This damage mechanism dominates + % over the thermal damage mechanism for times exceeding 10 + % seconds. (Page 15). + % + % The standard applies over the wavelength range from 300-700 nm. + % + fname = which('blueLightHazard.mat'); + blueHazard = ieReadSpectra(fname,wave); + % ieNewGraphWin; plot(wave,blueHazard) + + level = dLambda*dot(blueHazard,energy); + + % Equations 4.5a, 4.5b + val = false; + if duration <= 1e4 % seconds + if level*duration < 1e6, val = true; end + elseif level < 1e2 % Duration is long, so level must be low + val = true; + end + + otherwise + error('Unknown UV safety method %s\n',method); +end end \ No newline at end of file diff --git a/scripts/human/s_humanSafetyUVExposure.m b/scripts/human/s_humanSafetyUVExposure.m index a1012b9d..c4eb5d5b 100644 --- a/scripts/human/s_humanSafetyUVExposure.m +++ b/scripts/human/s_humanSafetyUVExposure.m @@ -34,71 +34,18 @@ % See also % -%{ -% Safety Notes -% -% Plot the three different functions and explain them here -% Make sure the formulae for hazards are implemented for -% -% Actinic UV hazard exposure limit skin and eye (4.3.1) -% Near UV hazard limit for the eye (4.3.2) -% Retinal blue light hazard exposure (4.3.3) -% Retinal blue light small source (4.3.4) -% Retinal thermal hazard (4.3.5) -% Retinal theermal hazard for weak visual stimulus (4.3.6) -% Infrared exposure for the eye (4.3.7) -% Thermal hazard for the skin (4.3.8) -%} - -%{ -We checked two ways if radiance -> irradiance is multiplied by 2pi or 1pi - -From: Peter B. Catrysse -Sent: Wednesday, August 21, 2019 11:02 AM -To: Joyce Eileen Farrell -Subject: RE: spectroradiometric measurements - -Hello Joyce, - - It is definitely pi. - There is 2*pi solid angle for the hemisphere, but when you integrate - you end up getting pi as factor. - -See you next week, - -Peter - -From: Joyce Eileen Farrell [mailto:jefarrel@stanford.edu] -Sent: Tuesday, August 20, 2019 10:58 PM -To: Peter Bert Catrysse -Subject: Re: spectroradiometric measurements - -Hi Peter, - -Thanks so very much for giving me the conversion from radiance to -irradiance. - -is it -E = pi*L/R (where E is irradiance, L is radiance and R is reflectance) -or -E=2pi*L/R - -See also this exchange - -https://physics.stackexchange.com/questions/116596/convert-units-for-spectral-irradiance - -The person multiplies b6000 by pi, not 2pi -%} +%% General parameters +wave = 300:700; %% 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; +ledRadiance = mean(radiance,2); +ledRadiance(wave > 500) = 0; % Only noise was measured above 500 nm + +lum405 = ieLuminanceFromEnergy(ledRadiance,wave); +plotRadiance(wave,ledRadiance); +irradiance = pi*ledRadiance; exposureMinutes = humanUVSafety(irradiance,wave); fprintf('Maximum exposure duration per eight hours: %f (min)\n',exposureMinutes) @@ -106,7 +53,6 @@ %% An example of the 385nm light in the OralEye camera fname = which('LED385nm.mat'); -wave = 300:700; radiance = ieReadSpectra(fname,wave); radiance = mean(radiance,2); plotRadiance(wave,radiance); @@ -120,7 +66,6 @@ %% The mean daylight we measured in California -wave = 300:700; [radiance,wave] = ieReadSpectra('DaylightPsychBldg.mat',wave); plotRadiance(wave,radiance); @@ -132,21 +77,71 @@ %% If you only know the luminance of an LED (monochromatic) and its bandwidth (s.d.) -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'}) +% The band width matters a lot for matching the curves. The luminance +% always matches correctly. +lum = lum405; % cd/m2, luminance of the 405 LED +thisWave = 405; % nm, center wavelength of the LED +bandwidth = 12; % nm, Gaussian standard deviation, FW at roughly 1/2 of the max +[estRadiance,estWave] = ieLuminance2Radiance(lum,thisWave,'sd',bandwidth); +plotRadiance(estWave,estRadiance); -ieLuminanceFromEnergy(estRadiance,wave) -ieLuminanceFromEnergy(radiance,wave) +% Check the luminance match +assert(ieLuminanceFromEnergy(estRadiance,estWave) - lum405 < 1e-10) -irradiance = pi*estRadiance; -exposureMinutes = humanUVSafety(irradiance,wave); +% Compare the curves +ieNewGraphWin; +plot(estWave,estRadiance,'o',wave,ledRadiance,'x'); +grid on; xlabel('Wave'); ylabel('Energy'); +legend({'estimated','measured'}); -% Conver the hazard energy into maximum daily allowable exposure in minutes +% Calculate the safety +irradiance = pi*estRadiance; +exposureMinutes = humanUVSafety(irradiance,estWave); 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)); -%% +%% Now check the other safety metric from Section 4.3.2 + +fname = which('LED385nm.mat'); +radiance = ieReadSpectra(fname,wave); +radiance = mean(radiance,2); + +irradiance = pi*radiance; +duration = 8; % Secs +safe = humanUVSafety(irradiance,wave,'method','eye','duration',duration); + +if safe + fprintf('Safe at a duration of %d secs\n',duration); +else + fprintf('Not safe at a duration of %d secs\n',duration); +end + +%% Blue light hazard calculation +fname = which('LED405nm.mat'); +radiance = ieReadSpectra(fname,wave); +radiance = mean(radiance,2); + +duration = 8*60*60; % Secs +safe = humanUVSafety(radiance,wave,'method','blue hazard','duration',duration); + +if safe + fprintf('Blue hazard: Safe at a duration of %d secs\n',duration); +else + fprintf('Blue hazard: Not safe at a duration of %d secs\n',duration); +end + +%% This is an extremely bright blackbody with a lot of short wave energy + +radiance = blackbody(wave,9000)*1e4; +plotRadiance(wave,radiance); + +duration = 1000; % Secs +safe = humanUVSafety(radiance,wave,'method','blue hazard','duration',duration); + +if safe + fprintf('Blue hazard: Safe at a duration of %d secs\n',duration); +else + fprintf('Blue hazard: Not safe at a duration of %d secs\n',duration); +end + +%% \ No newline at end of file