% simple octave script to explore the P(I) curve


% first compare P(I) curve for a range of I values
sst=10; vmax=0.6*(1.066).^sst ;
par=.5; alpha=0.025;
r=0:10:350;  % solar radiation 
ri =r*par ;  % par correction
v1= alpha.*ri; 
% Evans and Parslow
f_EP = vmax.*v1 ./ sqrt( vmax^2 + (v1).^2) ;
% Old equation
f_old = vmax.*v1 ./ sqrt( vmax + (v1).^2) ;
% Brian Griffiths and the one used in model
f_BG = vmax*(1-exp(-1.*v1/vmax) ) ;
plot(r,f_EP,";EP;",r,f_old,";old;",r,f_BG,";BG;")

pause
% now using BG equation explore the difference between using 
% one daily averaged value or a diurnal cycle
% calculation assumes only half the day has light with linear increase
% and decrease.

imax=10;
for k=1:66
imax=imax+20;
 i=1:1:12;
 irad = imax*abs(6.5 .- i )/6.5 ; 
 imean = sum(irad)/24.;
  v1= par*alpha.*irad;
  v1m = par*alpha*imean;
  f_BG = vmax*(1-exp(-1.*v1/vmax) ) ;
  f_BGd(k) = sum(f_BG)/24. ;
  f_BGm(k) = vmax*(1-exp(-1*v1m/vmax) ); 
  iave(k) = imean;
end

plot(iave,f_BGd,'r',iave,f_BGm,'b', iave,f_BGd ./ f_BGm ,'x')
ratio = f_BGd ./ f_BGm ;

%results: At low light two calculations are similar and at high
% light the f_BG(daily) is half the f_BG(mean)
