← Back to branch summary

~ubuntu-branches/ubuntu/trusty/psychtoolbox-3/trusty-proposed

« back to all changes in this revision

Viewing changes to Psychtoolbox/PsychDemos/IsomerizationsInEyeDemo.m

  • Committer: Package Import Robot
  • Author(s): Yaroslav Halchenko
  • Date: 2013-11-19 23:34:50 UTC
  • mfrom: (3.1.4 experimental)
  • Revision ID: package-import@ubuntu.com-20131119233450-f7nf92vb8qavjmk8
Tags: 3.0.11.20131017.dfsg1-3
Upload to unsable since fresh glew has arrived to sid!

Show diffs side-by-side

added added

removed removed

Lines of Context:
Psychtoolbox/PsychDemos/IsomerizationsInEyeDemo.m
6
6
% in watts/sr-m^2-wlinterval, or with a relative spectrum
7
7
% and a photopic troland value.
8
8
%
 
9
% NOTE, DHB, 7/19/13. This demo routine and its associated data routines 
 
10
% (DefaultPhotoreceptors, FillInPhotoreceptors, PrintPhotoreceptors)
 
11
% should be better integrated with the more recent code that
 
12
% implements the CIE physiological cone fundamentals, and the
 
13
% whole set of stuff should be better documented.  See also
 
14
%    IsomerizationsInDishDemo
 
15
%    CIEConeFundamentalsTest
 
16
%    ComputeCIEConeFundamentals
 
17
%    ComputeRawConeFundamentals
 
18
%    DefaultPhotoreceptors
 
19
%    FillInPhotoreceptors
 
20
%    PrintPhotoreceptors
 
21
%    RetIrradianceToIsoRecSec
 
22
% In particular, there should be some default for the 
 
23
% photoreceptors structure that gives one the CIE cone
 
24
% fundamentals in all their parametric glory, plus additional
 
25
% parameters that yield real energy/quantal sensitivites so
 
26
% that the resulting coordinates are isomerization rates in
 
27
% real units.  I think that we're close to having that, but
 
28
% better documentation and tidying is needed.
 
29
%
9
30
% 07/08/03 dhb  Wrote starting from IsomerizationsInDishDemo.
10
31
% 07/11/03 dhb  Grab data through subroutines.  Get rid of integration time.
11
32
% 07/15/03 dhb  Take eye size from function.
13
34
%               this is done, but the template may be useful someday.
14
35
% 03/20/12 dhb  Update cal file for PTB 3.
15
36
% 04/09/12 dhb  Add test of irradiance to troland conversion.
16
 
 
17
 
% Clear
18
 
clear all; close all;
19
 
 
20
 
% Set some photoreceptor properties.
21
 
% whatCalc = 'LivingDog';
22
 
whatCalc = 'LivingHumanFovea';
23
 
 
 
37
% 04/27/13 dhb  More extensive comments.
 
38
% 7/19/13  dhb  Print out photoreceptors structure using PrintPhotoreceptors.
 
39
%          dhb  Add monochromatic light option to the section that starts with trolands.
 
40
% 8/11/13  dhb  Add test of AborbtanceToAbsorbance.
 
41
%          dhb  Protect against case when absorbance is provided directly.
 
42
 
 
43
 
 
44
%% Clear
 
45
clear; close all;
 
46
 
 
47
%% Set photoreceptor properties.
 
48
%
 
49
% The photoreceptors structure gets filled with
 
50
% key parameters values (pupil size, eye length,
 
51
% pre-retinal absorbance, etc.)
 
52
%
 
53
% The routine DefaultPhotoreceptors is a high level
 
54
% call.  It fills in the 'source' fields and some
 
55
% values according to high-level descriptor (e.g.,
 
56
% ('CIE2Deg').  See help for that routine
 
57
% for available options.
 
58
%
 
59
% The routine FillInPhotoreceptors fetches the actual
 
60
% values for various fields, depending on the source.
 
61
%
 
62
% To get a feel for this, check what is in the photoreceptors
 
63
% structure after the first call, and then after the second.
 
64
whatCalc = 'CIE2Deg';
24
65
photoreceptors = DefaultPhotoreceptors(whatCalc);
 
66
photoreceptors.eyeLengthMM.source = 'LeGrand';
25
67
photoreceptors = FillInPhotoreceptors(photoreceptors);
26
68
 
27
 
% Define common wavelength sampling for this script.
 
69
%% Check AbsorbtancetoAbsorbance
 
70
%
 
71
% Simple check that this routine does what we expect, since
 
72
% we never use it anywhere else and this seems like as good
 
73
% a place to test it as any.
 
74
%
 
75
% We omit the normalization, because sometimes the wavelength
 
76
% sampling we use leads to a maximimum initial absorbance that
 
77
% is not unity, and letting AbsorbtanceToAbsorbance normalize
 
78
% causes disagreement.
 
79
testAbsorbance = photoreceptors.absorbance;
 
80
testAbsorbtance = photoreceptors.absorbtance;
 
81
checkAbsorbance = AbsorbtanceToAbsorbance(testAbsorbtance, photoreceptors.nomogram.S, photoreceptors.axialDensity.value,false);
 
82
diffs = testAbsorbance-checkAbsorbance;
 
83
if (max(abs(diffs(:))) > 1e-7)
 
84
    error('Cannot properly invert absorbance/absorbtance computations');
 
85
end
 
86
 
 
87
%% Define common wavelength sampling for this script.
 
88
% 
 
89
% S is [start delta nsamples] for the wavelengths in nm.
 
90
% This is standard PTB convention.
28
91
S = photoreceptors.nomogram.S;
29
92
 
30
 
% Get light spectrum.  You can choose various
31
 
% illustrative examples.
32
 
whichLight = 'fromMonitorRadiance';
33
 
%whichLight = 'fromTrolands';
34
 
 
35
 
% Here we'll start with a xenon arc lamp relative spectrum
36
 
switch (whichLight)
37
 
        % Convert from a specification of trolands and relative spectral power distribution.
38
 
        % Computation to retinal irradiance is done two ways just for fun.
39
 
        case 'fromTrolands',
40
 
                trolands = 100;
 
93
%% XYZ color matching functions
 
94
load T_xyz1931
 
95
T_xyz = SplineCmf(S_xyz1931,683*T_xyz1931,S);
 
96
T_Y = T_xyz(2,:);
 
97
        
 
98
%% Get light spectrum.  You can choose various illustrative examples.
 
99
%
 
100
% Available options:
 
101
%  'fromTrolands'
 
102
%  'fromMonitorRadiance'
 
103
%  'fromUniformQuantalSpd'
 
104
whichInputType = 'fromTrolands';
 
105
switch (whichInputType)
 
106
    
 
107
    % Start with troland value and a relative spectrum
 
108
    case 'fromTrolands'
 
109
        fprintf('Computing from troland value and relative spectrum\n');
 
110
        % Give troland value and type
 
111
        %
 
112
        % Type may be 'Photopic', 'Scotopic', or 'JuddVos'
 
113
        trolands = 1;
41
114
        trolandType = 'Photopic';
42
 
                load spd_xenonArc;
43
 
                irradianceWatts = TrolandsToRetIrradiance(spd_xenonArc,S_xenonArc,trolands, ...
 
115
        switch (trolandType)
 
116
            case 'Photopic'
 
117
                fprintf('Using photopic trolands\n');
 
118
            case 'Scotopic'
 
119
                fprintf('Using scotopic trolands\n');
 
120
            case 'JuddVos'
 
121
                fprintf('Using Judd-Vos luminosity function in troland calculations\n');
 
122
            otherwise
 
123
                fprintf('Unknown troland type specified');
 
124
        end
 
125
        
 
126
        %% Pupil.
 
127
        %
 
128
        % We do these computations for a fixed pupil size, ignoring the
 
129
        % pupilDiamter field of the photoreceptors structure.  That field
 
130
        % is set up to specify a source formula that estimates pupil diameter
 
131
        % from luminance.
 
132
        %
 
133
        % Since we are starting in trolands, the pupil size shouldn't actually
 
134
        % effect the calculations, except for finding the radiance that is
 
135
        % equivalent to the specified troland value.  
 
136
        % 
 
137
        % We remove the pupilDiameter.source field to make sure we aren't sending
 
138
        % mixed messages about how we want to handle pupil diameter.
 
139
        photoreceptors.pupilDiameter.value = 2;
 
140
        if (isfield(photoreceptors.pupilDiameter,'source'))
 
141
            photoreceptors.pupilDiameter = rmfield(photoreceptors.pupilDiameter,'source');
 
142
        end
 
143
        pupilAreaMm2 = pi*(photoreceptors.pupilDiameter.value/2)^2;
 
144
        
 
145
        % Specify relative spectrum to be used in
 
146
        % conversion to a full spectrum.
 
147
        %
 
148
        % Choices are:
 
149
        %  'Monochromatic'
 
150
        %  'XenonArc'
 
151
        %
 
152
        % If type is 'Monochromatic', must specify
 
153
        % wavelengthNm.
 
154
        spectrumType = 'Monochromatic';
 
155
        switch (spectrumType)
 
156
            case 'Monochromatic'
 
157
                monoWavelengthNm = 550;  
 
158
                wls = SToWls(S);
 
159
                monoWavelengthIndex = find(wls == monoWavelengthNm);
 
160
                if (isempty(monoWavelengthIndex))
 
161
                    error('No sample wavelength matches desired wavelength');
 
162
                end
 
163
                spd_fromTrolands = zeros(size(wls));
 
164
                spd_fromTrolands(monoWavelengthIndex) = 1;
 
165
                fprintf('Using monochromatic %d nm light as relative spectrum\n',monoWavelengthNm);
 
166
            case 'XenonArc'
 
167
                load spd_xenonArc;
 
168
                spd_fromTrolands = SplineSpd(S_xenonArc,spd_xenonArc,S);
 
169
                clear S_xenonArc spd_xenonArc
 
170
                fprintf('Using the spectrum of a xenon arc lamp as relative spectrum\n');
 
171
            otherwise
 
172
                error('Unknown spectrum type specified');
 
173
        end
 
174
        
 
175
        % Convert trolands back to spectral retinal irradiance.  This
 
176
        % depends on the pupil size and eye length specified.
 
177
                irradianceWattsPerUm2 = TrolandsToRetIrradiance(spd_fromTrolands,S,trolands, ...
44
178
                        trolandType,photoreceptors.species,photoreceptors.eyeLengthMM.value);
45
 
        irradianceTrolandsCheck = RetIrradianceToTrolands(irradianceWatts,S_xenonArc,trolandType, ...
 
179
        irradianceTrolandsCheck = RetIrradianceToTrolands(irradianceWattsPerUm2,S,trolandType, ...
46
180
            photoreceptors.species,photoreceptors.eyeLengthMM.value);
47
181
        trolandsCheck = sum(irradianceTrolandsCheck);
48
 
        fprintf('\nInput trolands is %0.1f, checked value is %0.1f\n\n',trolands,trolandsCheck);
49
 
                irradianceWatts = SplineSpd(S_xenonArc,irradianceWatts,S);
 
182
        fprintf('Input troland value is %0.1f, checked value is %0.1f\n',trolands,trolandsCheck);
 
183
                irradianceWattsPerUm2 = SplineSpd(S,irradianceWattsPerUm2,S);
50
184
        
51
185
                % Another way to do this calculation.  Pupil size should cancel out.  Should get
52
 
            % same answer as above.
53
 
                pupilSizeMM = 2;
54
 
                pupilAreaMM = pi*(pupilSizeMM/2)^2;
55
 
                luminance = TrolandsToLum(trolands,pupilAreaMM);
56
 
                radianceWatts = LumToRadiance(spd_xenonArc,S_xenonArc,luminance,trolandType);
57
 
                irradianceWattsCheck = RadianceToRetIrradiance(radianceWatts,S_xenonArc,pupilAreaMM,photoreceptors.eyeLengthMM.value);
 
186
            % same answer as above.  This has as a byproduct computing a stimulus radiance,
 
187
        % which is useful for some of the common printout below.
 
188
        luminanceCdM2FromTrolands = TrolandsToLum(trolands,pupilAreaMm2);
 
189
                radianceWattsPerM2Sr = LumToRadiance(spd_fromTrolands,S,luminanceCdM2FromTrolands,trolandType); 
 
190
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;
 
191
                irradianceWattsCheck = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value);
58
192
                figure(1); clf; hold on
59
 
                set(plot(SToWls(S),irradianceWatts,'r'),'LineWidth',2);
60
 
                set(plot(SToWls(S_xenonArc),irradianceWattsCheck,'k'),'LineWidth',2);
 
193
                set(plot(SToWls(S),irradianceWattsPerUm2,'r'),'LineWidth',2);
 
194
                set(plot(SToWls(S),irradianceWattsCheck,'k'),'LineWidth',2);
61
195
                set(title('Check of trolands to irradiance calculation'),'FontSize',14);
62
196
                set(xlabel('Wavelength (mm)'),'FontSize',14);
63
197
                set(ylabel('Irradiance'),'FontSize',14);
64
 
 
 
198
        
 
199
        % For case if monochromatic light, can check retinal irradiance against
 
200
        % the direct formulae provided in W&S, 2cd edition, p. 105, eqs 2.4.4 
 
201
        % Note that these equations are missing a factor of
 
202
        % the wavelength in the numerator for the quantal conversions, as
 
203
        % pointed out by Makous in his 1997 JOSA paper.
 
204
        %
 
205
        % This only works for 'Photopic' and 'Scotopic' trolands.
 
206
        %
 
207
        % There must be an eye length implicit in this calculation. Using
 
208
        % the LeGrand model eye length gives good agreement between these
 
209
        % and our more general calculations done in the main section below.
 
210
        %
 
211
        % Makous (1997), JOSA A, 14, p. 2331 gives retinal illuminances of
 
212
        % 5.44 quanta/[um2-sec] for 1 scotopic troland, and 14.65 quanta/[um2-sec]
 
213
        % for 1 photopic troland, with the calulations specified for 510 nm.
 
214
        % The calculations here, done in several different ways, yield
 
215
        % 5.442 (scotopic, agrees) and 26.85 (photopic, does not agree).  But at 
 
216
        % 550 nm, the code here yields 14.64 quanta/[um2-sec], which seems close
 
217
        % enough to the provided 14.65 to make me think that Makous' value is
 
218
        % actually for a wavelength close to 550 nm.  That would be a more typical
 
219
        % wavelength at which to do photopic calculations, despite what the text
 
220
        % in the paper says.
 
221
        %
 
222
        % Note that there are also errors in Tables 3 and 4 of the Makous
 
223
        % paper, and that corrected values appear in Tables 3 and 2 of
 
224
        % Makous (2004), Scotopic vision, In The Visual Neurosciences,
 
225
        % Werner and Chalupa (eds).  What does not seem to be specified
 
226
        % in either place is the wavelengths used in the calculations of
 
227
        % the two tables.
 
228
        if (strcmp(spectrumType,'Monochromatic'))
 
229
            switch (trolandType)
 
230
                case 'Photopic'
 
231
                    load T_xyz1931;
 
232
                    T_vLambda = SplineCmf(S_xyz1931,T_xyz1931(2,:),S);
 
233
                    clear T_xyz1931 S_xyz1931
 
234
                    magicFactorE = 5.261e-12;
 
235
                    magicFactorQ = 2.649e13;                  
 
236
                case 'Scotopic'
 
237
                    load T_rods;
 
238
                    T_vLambda = SplineCmf(S_rods,T_rods,S);
 
239
                    magicFactorE = 2.114e-12;
 
240
                    magicFactorQ = 1.064e13;
 
241
            end
 
242
        end
 
243
        if (strcmp(spectrumType,'Monochromatic'))
 
244
            switch (trolandType)
 
245
                case {'Photopic','Scotopic'}
 
246
                    irradianceDirectWattsPerMm2Check = magicFactorE*trolands/T_vLambda(monoWavelengthIndex);
 
247
                    irradianceDirectQuantaPerMm2SecCheck = 1e-9*monoWavelengthNm*magicFactorQ*trolands/T_vLambda(monoWavelengthIndex);
 
248
                    fprintf('Retinal irradiance in area units computed from from trolands via Wyszecki and Stiles formulae\n');
 
249
                    fprintf('\t%0.4g Watts/mm2\n\t%0.4g quanta/[mm2-sec]\n',sum(irradianceDirectWattsPerMm2Check),sum(irradianceDirectQuantaPerMm2SecCheck));
 
250
            end
 
251
        end
 
252
        
65
253
        % Start with radiance measurements, which we just
66
254
        % pull out of the Toolbox's default calibration file.
67
255
        case 'fromMonitorRadiance'
68
256
                % Load light radiance.  We'll use a monitor white.
69
257
                % The original units are watts/sr-m^2-wlinterval.
70
258
                cal = LoadCalFile('PTB3TestCal');
71
 
                radianceWatts = SplineSpd(cal.S_device,sum(cal.P_device,2),S);
 
259
                radianceWattsPerM2Sr = SplineSpd(cal.S_device,sum(cal.P_device,2),S);
72
260
                
73
261
                % Find pupil area, needed to get retinal irradiance.  We compute
74
 
                % pupil area based on the luminance of stimulus.  
75
 
                load T_xyz1931
76
 
                T_xyz = SplineCmf(S_xyz1931,683*T_xyz1931,S);
77
 
                theXYZ = T_xyz*radianceWatts; theLuminance = theXYZ(2);
78
 
                [nil,pupilAreaMM] = PupilDiameterFromLum(theLuminance,photoreceptors.pupilDiameter.source);
 
262
                % pupil area based on the luminance of stimulus according to the
 
263
        % algorithm specified in the photoreceptors structure.
 
264
                theXYZ = T_xyz*radianceWattsPerM2Sr; theLuminance = theXYZ(2);
 
265
                [nil,pupilAreaMm2] = PupilDiameterFromLum(theLuminance,photoreceptors.pupilDiameter.source);
 
266
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;
79
267
                
80
268
                % Convert radiance of source to retinal irradiance and convert to quantal units.
81
 
                irradianceWatts = RadianceToRetIrradiance(radianceWatts,S, ...
82
 
                        pupilAreaMM,photoreceptors.eyeLengthMM.value);
 
269
                irradianceWattsPerUm2 = RadianceToRetIrradiance(radianceWattsPerM2Sr,S, ...
 
270
                        pupilAreaMm2,photoreceptors.eyeLengthMM.value);
83
271
 
84
 
        % This is here to match a parameterization the Brian Wandell supplied
 
272
        % This light as well as some parameter tweaking are here to match a parameterization that Brian Wandell supplied
85
273
        % to match what his code to do these computations produces.  Note also
86
274
        % the mucking with the photoreceptors structure.  Wandell estimates
87
 
        % L, M, S isomerizations/cone-sec of 16.5, 12.68, 2.27.
 
275
        % L, M, S isomerizations/cone-sec of 16.5, 12.68, 2.27.  These are very close to the numbers
 
276
    % we get here.
88
277
        case 'fromUniformQuantalSpd',
89
278
                % Load corneal cone sensitivities in energy units, convert to quantal sensitivities
90
279
                % and set specified peak absorbtance.
98
287
                % load T_cones_ss10; T_cones = T_cones_ss10; S_cones = S_cones_ss10;
99
288
                % load T_cones_smj; T_cones = T_cones_smj; S_cones = S_cones_smj;
100
289
                % load T_cones_sp; T_cones = T_cones_sp; S_cones = S_cones_sp;
 
290
        
101
291
                peakIsomerizationEfficiency = [0.27 0.23 0.07]';
102
292
                T_cones = SplineCmf(S_cones,QuantaToEnergy(S_cones,T_cones')',S);
103
293
                T_cones(1,:) = T_cones(1,:)/max(T_cones(1,:));
106
296
                T_cones = diag(peakIsomerizationEfficiency)*T_cones;
107
297
                photoreceptors.isomerizationAbsorbtance = T_cones;
108
298
 
109
 
                % Get spectral luminous efficiency function
110
 
                load T_xyz1931; T_xyz = T_xyz1931; S_xyz = S_xyz1931;
111
 
                % load T_xyzJuddVos; T_xyz = T_xyzJuddVos; S_xyz = S_xyzJuddVos;
112
 
                T_Y = 683*SplineCmf(S_xyz,T_xyz(2,:),S);
113
 
 
114
299
                % Create a spectrally uniform spd (in quantal units), and convert
115
300
                % to energy units.
116
301
                uniformSpd = QuantaToEnergy(S,ones(S(3),1));
117
302
 
118
303
                % Normalize to radiance corresponding to 1 cd/m2.
119
304
                normConst = T_Y*uniformSpd;
120
 
                radianceWatts = uniformSpd/normConst;
 
305
                radianceWattsPerM2Sr = uniformSpd/normConst;
 
306
        photopicLuminanceCdM2 = T_Y*radianceWattsPerM2Sr;
121
307
 
122
 
                % Set pupil diameter for 1mm2 pupil area, photoreceptor diameter for 4mm2 collecting
 
308
                % Set pupil diameter for 1 mm2 pupil area, photoreceptor diameter for 4 mm2 collecting
123
309
                % area.  Set eye length to 17 mm.
124
310
                photoreceptors.pupilDiameter.value = 2*sqrt(1/pi);
125
 
                pupilAreaMM = pi*(photoreceptors.pupilDiameter.value/2)^2;
 
311
                pupilAreaMm2 = pi*(photoreceptors.pupilDiameter.value/2)^2;
126
312
                photoreceptors.ISdiameter.value = [2*sqrt(4/pi) 2*sqrt(4/pi) 2*sqrt(4/pi)]';
127
313
                photoreceptors.eyeLengthMM.value = 17;
128
 
                irradianceWatts = RadianceToRetIrradiance(radianceWatts,S,pupilAreaMM,photoreceptors.eyeLengthMM.value );
129
 
end
130
 
 
131
 
irradianceQuanta = EnergyToQuanta(S,irradianceWatts);
 
314
                irradianceWattsPerUm2 = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value );
 
315
end
 
316
 
 
317
%% Print out a whole bunch of quantities that are equivalent to the radiance, given
 
318
% other eye parameters.
 
319
radianceWattsPerCm2Sr = (10.^-4)*radianceWattsPerM2Sr;
 
320
radianceQuantaPerCm2SrSec = EnergyToQuanta(S,radianceWattsPerCm2Sr);
 
321
degPerMm = RetinalMMToDegrees(1,photoreceptors.eyeLengthMM.value);
 
322
irradianceWattsPerUm2Check = RadianceToRetIrradiance(radianceWattsPerM2Sr,S,pupilAreaMm2,photoreceptors.eyeLengthMM.value);
 
323
if (any(abs(irradianceWattsPerUm2 - irradianceWattsPerUm2Check) > 1e-10))
 
324
    error('Back computation of retinal irradiance from radiance does not check');
 
325
end
 
326
irradianceScotTrolands = RetIrradianceToTrolands(irradianceWattsPerUm2, S, 'Scotopic', [], num2str(photoreceptors.eyeLengthMM.value));
 
327
irradiancePhotTrolands = RetIrradianceToTrolands(irradianceWattsPerUm2, S, 'Photopic', [], num2str(photoreceptors.eyeLengthMM.value));
 
328
irradianceQuantaPerUm2Sec = EnergyToQuanta(S,irradianceWattsPerUm2);
 
329
irradianceWattsPerCm2 = (10.^8)*irradianceWattsPerUm2;
 
330
irradianceQuantaPerCm2Sec = (10.^8)*irradianceQuantaPerUm2Sec;
 
331
irradianceQuantaPerMm2Sec = (10.^-2)*irradianceQuantaPerCm2Sec;
 
332
irradianceQuantaPerUm2Sec = (10.^-6)*irradianceQuantaPerMm2Sec;
 
333
irradianceQuantaPerDeg2Sec = (degPerMm^2)*irradianceQuantaPerMm2Sec;
 
334
 
 
335
% Print out photoreceptor stucture information
 
336
fprintf('\n');
 
337
PrintPhotoreceptors(photoreceptors);
 
338
fprintf('\n');
 
339
 
 
340
% Radiometric iformation
 
341
fprintf('Luminance %0.3f cd/m2\n',photopicLuminanceCdM2);
 
342
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) watts/cm2\n',sum(irradianceWattsPerCm2),log10(sum(irradianceWattsPerCm2)));
 
343
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) watts/mm2\n',1e-2*sum(irradianceWattsPerCm2),log10(sum(irradianceWattsPerCm2)));
 
344
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[cm2-sec]\n',sum(irradianceQuantaPerCm2Sec),log10(sum(irradianceQuantaPerCm2Sec)));
 
345
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[mm2-sec]\n',sum(irradianceQuantaPerMm2Sec),log10(sum(irradianceQuantaPerMm2Sec)));
 
346
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[um2-sec]\n',sum(irradianceQuantaPerUm2Sec),log10(sum(irradianceQuantaPerUm2Sec)));
 
347
fprintf('Stimulus retinal irradiance %0.4g (%0.1f log10) quanta/[deg2-sec]\n',sum(irradianceQuantaPerDeg2Sec),log10(sum(irradianceQuantaPerDeg2Sec)));
 
348
fprintf('\n');
 
349
        
 
350
%% Get retinal irradiance in quanta/[sec-um^2-wlinterval]
 
351
irradianceQuanta = EnergyToQuanta(S,irradianceWattsPerUm2);
132
352
figure(2); clf; set(gcf,'Position',[100 400 700 300]);
133
353
subplot(1,2,1); hold on
134
354
set(plot(SToWls(S),irradianceQuanta,'r'),'LineWidth',2);
136
356
set(xlabel('Wavelength (nm)'),'FontSize',12);
137
357
set(ylabel('Quanta/sec-um^2-wlinterval'),'FontSize',12);
138
358
 
139
 
% Do the work in toolbox function
 
359
%% Do the work in toolbox function
140
360
[isoPerConeSec,absPerConeSec,photoreceptors] = ...
141
 
        RetIrradianceToIsoRecSec(irradianceWatts,S,photoreceptors);
 
361
        RetIrradianceToIsoRecSec(irradianceWattsPerUm2,S,photoreceptors);
142
362
 
143
363
% Make a plot showing the effective photoreceptor sensitivities in quantal
144
364
% units, expressed as probability of isomerization.
155
375
fprintf('***********************************************\n');
156
376
fprintf('Isomerization calculations for living human retina\n');
157
377
fprintf('\n');
158
 
fprintf('Calculations done using:\n');
159
 
fprintf('\t%s estimates for photoreceptor IS diameter\n',photoreceptors.ISdiameter.source);
160
 
fprintf('\t%s estimates for photoreceptor OS length\n',photoreceptors.OSlength.source);
161
 
fprintf('\t%s estimates for receptor specific density\n',photoreceptors.specificDensity.source);
162
 
fprintf('\t%s photopigment nomogram\n',photoreceptors.nomogram.source);
163
 
fprintf('\t%s estimates for lens density\n',photoreceptors.lensDensity.source);
164
 
fprintf('\t%s estimates for macular pigment density\n',photoreceptors.macularPigmentDensity.source);
165
 
fprintf('\t%s method for pupil diameter calculation\n',photoreceptors.pupilDiameter.source);
166
 
fprintf('\t%s estimate (%g mm) for axial length of eye\n',photoreceptors.eyeLengthMM.source,photoreceptors.eyeLengthMM.value);
167
 
fprintf('\n');
168
378
fprintf('Photoreceptor Type             |\t       L\t       M\t     S\n');
169
379
fprintf('______________________________________________________________________________________\n');
170
380
fprintf('\n');
171
 
fprintf('Lambda max                     |\t%8.1f\t%8.1f\t%8.1f\t nm\n',photoreceptors.nomogram.lambdaMax);
 
381
if (isfield(photoreceptors.nomogram,'lambdaMax'))
 
382
    fprintf('Lambda max                     |\t%8.1f\t%8.1f\t%8.1f\t nm\n',photoreceptors.nomogram.lambdaMax);
 
383
end
172
384
fprintf('Outer Segment Length           |\t%8.1f\t%8.1f\t%8.1f\t um\n',photoreceptors.OSlength.value);
173
385
fprintf('Inner Segment Diameter         |\t%8.1f\t%8.1f\t%8.1f\t um\n',photoreceptors.ISdiameter.value);
174
386
fprintf('\n');
201
413
%     otherwise
202
414
% end
203
415
 
 
416