%% ------------------------------------------------------------------------
%    Determine range of spatiotemporal frequencies where penumbral-cone,
%              luminance contrast is below luminance threshold
%  ------------------------------------------------------------------------
% From "What can the eye see with melanopsin?" by Nugent, T.W. & Zele, A.J.
% Software S1.
%
% Instructions to run:
%   1. This script uses the luminance spatiotemporal threshold surface from
%      Kelly 1979, Eq. 5-8, doi: https://doi.org/10.1364/JOSA.69.001340
%   2. Ensure the file with the system spectra and haemoglobin absorption
%      (Software S2) is called "DatasetS2_SystemSpectra.xlsx" or change the
%      spectrumDataFileName to match the filename.
%   3. Calculate the maximum penumbral-cone, luminance contrast in the
%      silent substitution conditions (i.e., the melanopsin-isolated
%      direction). Type this number into line 29: penumbralConeContrast
%      Method for this calculation is in Nugent et. al. 2023,
%      step 18, doi: https://doi.org/10.1016/j.xpro.2023.102126
%   4. If you want to use this script to calculate penumbral cone
%      luminance contrast error, then add own spectra to
%      "spectrumDataFileName" path and change "condition1" and "condition2"
%      relative spectral powers below.
%
%% Clear workspace and set variables
clearvars; close all; clc;

% Name of system's spectrum file:
spectrumDataFileName = "DatasetS2_SystemSpectra.xlsx";

% This is the calculated penumbral cone intrusion for our given system
% (Accounting for larger individual observer penumbral cone luminance
% contrast).
penumbralConeContrast = 0.033;      % As a fraction, i.e., percent / 100

% What range of spatiotemporal frequencies for Kelly model
Fs_min = 0.08;              % cpd
Fs_max = 20;                % cpd
Ft_min = 0.1;               % Hz
Ft_max = 50;                % Hz
NumSpatialSamples = 46;     % samples
NumTemporalSamples = 54;    % samples

% What range of spatiotemporal frequencies for penumbral cone contrast
% calculations
fs_min = 0.08;              % cpd
fs_max = 0.5;               % cpd
ft_min = 0.1;               % Hz
ft_max = 2;                 % Hz
numSpatialSamples = 1e2;    % samples
numTemporalSamples = 1e2;   % samples

% If you want to calculate penumbral cone luminance errors
% Then enter the relative primary powers here (percentage of max power)

% This would be photoreceptor-low
condition1 = [1; 0; 0.794289206285317; 0; 0.786386436519884];
% This would be photoreceptor-high
condition2 = [0.394043887337721; 0.705478635178566; 0.306950337867556; 0.990538467088180; 0];

%% Create the spatiotemporal threshold surface from Kelly 1979
% This part is not needed for calculations actually

% Create the spatial and temporal vectors
% 1D vectors
alpha = 2*pi*(10.^(linspace(log10(Fs_min), log10(Fs_max), NumSpatialSamples))).';
ft = 10.^linspace(log10(Fs_min),log10(Ft_max), NumTemporalSamples);          % vel = ft ./ spatial
% 2D mesh
[ft, alpha] = meshgrid(ft, alpha);

% Equations 5-8 in Kelly
vel = 2 * pi * ft ./ alpha;
k = 6.1 + 7.3 .* (abs(log10(vel ./ 3)).^3);
amax = 45.9./(vel + 2);
g = k .* vel .* alpha.^2 .* exp(-2.*alpha./amax);
                            % Contrast = 1/g, fs = alpha/2/pi

% Plot the Kelly model
figure()
plot3(alpha/2/pi, ft, 1./g, 'k')
set(gca, 'XScale', 'log', 'YScale', 'log', 'ZScale', 'log', 'ZDir', 'reverse')
zlim([0.004 1])
xlabel('Spatial [cpd]'), ylabel('Temporal [Hz]')
title('Spatiotemporal threshold surface of luminance')

%% Calculate penumbral cone error of system

% Load in the photoreceptor spectral sensitivities from CIE S026/E:2018,
% doi: doi.org/10.25039/S026.2018
spectralSensitivities = readtable(spectrumDataFileName, 'Sheet', 'Photoreceptor Action Spectra');
photoreceptors = spectralSensitivities.Properties.VariableNames(2:end);
photoreceptorColours = [0 0 1; 0 1 0; 1 0 0; 0 0 0; 0 1 1];

% Load in the retinal vasculature absorption spectra from Prahl, url:
% http://omlc.org/spectra/hemoglobin/summary.html
Hb = readtable(spectrumDataFileName, 'Sheet', 'Penumbral Shadow Absorption');
d = 5;       % Diameter of Venules and Artioles [um]
c = 150;        % Concentration Hb [g/L]
m = 64500;      % Molar concentration [mol/L]
Hb.AHbO2 = d * log(10) * Hb.eHbO2/10000 * c / m;
Hb.AHb = d * log(10) * Hb.eHb/10000 * c / m;
alpha_total = 0.85 * Hb.AHbO2 + 0.15 * Hb.AHb;
Hb.haemoglobinTransmission = 10.^-alpha_total;

% Load in the system spectra
spectrum = readtable(spectrumDataFileName, 'Sheet', 'Spectrum');
primaryNames = spectrum.Properties.VariableNames(2:end-1);
primaryColours = [97 0 255; 0 169 255; 129 255 0; 255 255 0; 180 0 0]./255;

% Plot the three spectral functions loaded in:
f = figure(); hold on
xlabel('Wavelength [nm]'), ylabel('Radiance [W/m2/sr]')
for ii = 1:size(photoreceptors, 2)
    area(spectralSensitivities.Wavelength, spectralSensitivities.(photoreceptors{ii}), 'FaceColor', photoreceptorColours(ii,:), 'FaceAlpha', 0.25, 'EdgeColor' ,'none')
end
plot(Hb.Wavelength, Hb.haemoglobinTransmission, 'k:')
xlim([400 700])
ax2 = axes(f, 'NextPlot', 'add', 'Color', 'none');
for ii = 1:size(primaryNames,2)
    area(spectrum.Wavelength, spectrum.(primaryNames{ii}), 'FaceColor', primaryColours(ii,:))
    xticks([])
    yticks([])
    xticklabels([])
    yticklabels([])
end
xlim([400 700])


%% Calculate penumbral cone intrusion in melanopsin isolated stimuli
dLambda = spectrum.Wavelength(2) - spectrum.Wavelength(1);

% Create an A-matrix to control photoreceptors
A = spectralSensitivities{:,2:end}.' * spectrum{:,2:end-1} .* dLambda;

% A-matrix under penumbral cones
A_penumbral = spectralSensitivities{:,2:end}.' * (spectrum{:,2:end-1} .* Hb.haemoglobinTransmission)  .* dLambda;

% Use A-matrix to calculate the filtered cone errors:
penumbralCone1 = A_penumbral * condition1;
penumbralCone2 = A_penumbral * condition2;
penumbralLuminanceContrast = (sum(penumbralCone2(2:3)) - sum(penumbralCone1(2:3))) ./ (sum(penumbralCone2(2:3)) + sum(penumbralCone1(2:3)));

fprintf('Calculated penumbral luminance contrast = %0.2f%%\n', 100 * penumbralLuminanceContrast);

%% Calculate spatiotemporal region free from penumbral cone intrusion
% Use the Kelly model to determine when luminance threshold is greater 
% than the penumbral-cone, luminance contrast in modulation

% From contrast to sensitivity
g = 1/penumbralConeContrast;

% Spatial frequencies to solve for:
alpha = 2*pi*(10.^(linspace(log10(fs_min), log10(fs_max), numSpatialSamples))).';

% Going to symbolically solve for ft, given g and alpha
syms a ft       % a = symbolic alpha

% Set up symbolic equations:
v = 2 * pi * ft / a;
k = (6.1 + 7.3 * abs(log10(v / 3))^3);
amax = 45.9 / (v + 2);
eqn = k * v * a^2 * exp(-2 * a / amax) == g;

% Substitute symbolic a for spatial frequencies in alpha
S = subs(eqn, a, alpha);
ft_limit = nan(size(alpha)).';
for ii = 1:length(alpha)
    try
        % There is a solution
        ft_limit(ii) = vpasolve(S(ii), ft);
    catch
        % There isn't a solution
        ft_limit(ii) = nan;
    end
end
ft_limit = double(ft_limit);

% alpha/2/pi and ft_max is the line intersecting the Kelly model at
% penumbralConeContrast. Below this is below penumbral cone error and vice
% versa.

%% Plot the range of spatiotemporal frequencies without penumbral cone
% error.
f = figure();
hold on
xPos = 10; yPos = 10; width = 4.5; height = 4.5;
set(f, 'Units', 'centimeters', 'Position', [xPos yPos width height])
fill([alpha(1); alpha; alpha(end); alpha(1)]./2./pi, [ft_max, ft_limit, ft_max, ft_max], [158 159 163]./255, 'FaceColor', [158 159 163]./255, 'FaceAlpha', 0.5)
fill([alpha(1); alpha; alpha(end); alpha(1)]./2./pi, [fs_min, ft_limit, fs_min, fs_min], [120 183 152]./255, 'FaceColor', [120 183 152]./255, 'FaceAlpha', 0.5)
plot(alpha/2/pi, ft_limit, 'k', 'LineWidth', 2)
text(0.16, 0.25, 'melanopsin-isolated', 'Color', 0.5*[120 183 152]./255,...
    'FontSize', 6, 'Units', 'data',...
    'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle')
text(0.3, 1.5, 'penumbral cones', 'Color', 0.5*[158 159 163]./255,...
    'FontSize', 6, 'Units', 'data',...
    'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle')
% This line is the max spatial frequencies melanopsin can see:
plot([0.4 0.4], [0.1 1.2], 'k--')
plot([0.4 0.4], [1.7 2], 'k--')
set(gca, 'XScale', 'log', 'YScale', 'log', 'FontSize', 6)
xlim([fs_min fs_max])
ylim([ft_min ft_max])
xlabel('Spatial [cpd]')
ylabel('Temporal [Hz]')
xticks([0.08 0.16 0.2 0.25 0.315 0.5])
yticks([0.1 0.25 0.5 1 2])
text(-0.22, 1.1, 'E', 'FontSize', 12, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'VerticalAlignment', 'top')