% figures    gr_sep09    REV. Jul 2012
%
clear all;
return

%% photoreceptors, contrast response, contrast threshold, Njnd

Lmn = 2;  Lmx = 500;   % min and max display luminance
L = linspace(Lmn,Lmx,1000);
cto = .01;  %  Barten
% figure(1);
for S = [10,50,100];        % fixed adaptation status
    P = L ./ (L+S);   % photor. response
%    h=semilogx(L,P); hold on; grid on; set(h,'LineWidth',1.5); title('response vs. luminance for S=10,50,100');
    cr = L .* S ./ (L+S).^2;   % (normalized) contrast response = dP/(dL/L)
%    h=semilogx(L,cr); hold on; grid on; set(h,'LineWidth',1.5); title('ncr vs. luminance for S=10,50,100');
    ct = (cto./4) * (1./cr);    % ct, scaled to get ct=cto for L=S;
%    h=loglog(L,ct); hold on; grid on; set(h,'LineWidth',1.5); title('ct vs. luminance for S=10,50,100'); axis([1 1e3 .009 .1]);
end;
% hold off;

figure(2);
Lamb = 0;  %  set  Lamb=0,1
L=L+Lamb; Lmn=Lmn+Lamb; Lmx=Lmx+Lamb;   %%%%%%%%%  beware... !
kk = [];  SS = [];
for S = [5:1:300] + Lamb;         % fixed adaptation status
   cr = L .* S ./ (L+S).^2;   % contrast response = dP/(dL/L)
   ct = (cto./4) * (1./cr);    % ct, scaled to get ct=cto for L=S;
   LL = Lmn;
   for k = 1:1000000
      ind = floor((LL - Lmn) / ((Lmx-Lmn)/1000)+1);
      LL = LL * (1+ct(ind));
      % LL = LL * (1+cto);   % ideal case, varied adaptation
      if LL>Lmx,   [k, LL, ind, ct(ind)];  break; end;
   end;
   kk = [kk k];   SS = [SS S];
end;
h=semilogx(SS,kk); grid on; set(h,'LineWidth',1.5);  % axis([1 1e3 240 360]); 
title('Njnd vs. S');   % title('Njnd vs. S  -- Lamb=0,1');
% hold on;   % to visually compare different Lamb values



%%  no. of JNDs in the range L1...L2, with Lamb
Lmn = .1;  Lmx = 500;   % min and max display luminance
L = linspace(Lmn,Lmx,1000);
Lamb = 0;  L=L+Lamb; Lmn=Lmn+Lamb; Lmx=Lmx+Lamb;   %%%%%%%%%  beware... !
cto = .01;  %  Barten
S = 50+Lamb;         % fixed adaptation state
cr = L .* S ./ (L+S).^2;   % contrast response = dP/(dL/L)
ct = (cto./4) * (1./cr);    % ct, scaled to get ct=cto for L=S;
L1 = 2+Lamb;  L2 = 10+Lamb;

LL = L1;
for k = 1:1000000
    ind = floor((LL - Lmn) / ((Lmx-Lmn)/1000)+1);
    % [k, LL, ind, ct(ind)]
    LL = LL * (1+ct(ind));
    if LL>L2,  break;  end;
end;
k   %  no. of JNDs in the range L1...L2, with  Lmn < L1 < L2 < Lmx

