% Calculations to duplicate the CMSs from Mark's spreadsheet
% Using Matlab 7.3.1.1266 (R2006b) on my Macbook results agree to 10^(-15)

% If this is going to be called a whole lot, parameterize the data from
% the data files.

function [lms_bar_norm]=cie2006cmfs(age,fs)

if nargin < 2
    age = 32.; fc = 2.;    % age and field size
else
    bad=0;
    if age < 20 age=20; bad=1; end
    if age > 80 age=80; bad=1; end
    if fs < 1 fs=1;     bad=1; end
    if fs > 10 fs=10;   bad=1; end
    if bad
        fprintf(1,'Parameters out of range, adjusted to: age= %d  field=%d \n',age,fs);
    end
end

% read in the raw data. Note first two points of docul are interpolated.
rmd = load('/Users/drwpci/Documents/CIE_CMFs/RelativeMacularDensity.txt');
LMSa = load('/Users/drwpci/Documents/CIE_CMFs/Alms.txt');
docul = load('/Users/drwpci/Documents/CIE_CMFs/docul.txt');

wl = [390:5:780]';

% field size corrected macular density
corrected_rmd =rmd.*(0.485*exp(-fs./6.132))';

% age corrected lens/ocular media density 

if (age<=60)
    correct_lomd = docul(:,1) .* (1+0.02*(age-32)) + docul(:,2);
else
    correct_lomd = docul(:,1) .* (1.56+0.0667*(age-60)) + docul(:,2);
end

% corrected LMS (no age correction)
alpha_lms = 0. * LMSa;
alpha_lms(:,1) = 1-10.^(-(0.38+0.54*exp(-fs/1.333))*(10.^LMSa(:,1)));    % column X = LMSa(1)
alpha_lms(:,2) = 1-10.^(-(0.38+0.54*exp(-fs/1.333))*(10.^LMSa(:,2)));    % column Y = LMSa(2)
alpha_lms(:,3) = 1-10.^(-(0.3+0.45*exp(-fs/1.333))*(10.^LMSa(:,3)));     % column Z = LMSa(3)

% this fix is required because the above math fails for alpha_lms(3,:)==0
alpha_lms(find(wl>=620),3) = 0;

% Corrected to Corneal Incidence
lms_barq = alpha_lms .* repmat(10.^(-corrected_rmd-correct_lomd),1,3);

% Corrected to Energy Terms
lms_bar = lms_barq .* repmat(wl,1,3);

% normalized
lms_bar_norm = lms_bar ./ repmat(max(lms_bar),79,1);