% Photosynthesis.m

% Main Routine for calling and acting on the Leaf Photosynthesis subroutine LeafPs

% 12/21/2011

% Dennis Baldocchi
% University of California, Berkeley
% baldocchi@berkeley.edu

% debug notes
% fixed some problems with units
% adopted Collatz transition between wc and wj



clear all;

%cd,'d:\matlab code';

%cd, 'C:\Users\Baldocchi\Documents\MATLAB'

%   Initialize Meteorology
	met.parin=2000;  %photosynthetic active radiation, par, umol m-2 s-1
    met.co2=380;     % CO2 concentration, umol/mol
    met.P_kPa = 101.;  % air pressure
	met.P_Pa=met.P_kPa*1000.;
	met.U = 3.;
	met.Tair_C=25;  % air temperature
	met.Tair_K=273.15+met.Tair_C;
	met.eair_kPa=2;  % atmospheric vapor pressure
    met.eair_Pa=met.eair_kPa*1000.;
    met.Tleaf_K=300;
    
   ps.rstom = 50.;     % stomatal resistance, s/m
    
    bound_layer_res.heat = 20;
    bound_layer_res.vapor = 20;
    bound_layer_res.co2 = 20;
    
    param.vcopt=127.;
    
    % function [ps]=LeafPs(Iphoton,cca, tlk, rb_co2, P_kPa, eair_Pa,vcopt)
    
  % [ps]=LeafPs(met.parin,met.co2,met.Tleaf_K, bound_layer_res.co2,met.P_kPa,met.eair_Pa,param.vcopt);
   
     
    [ps]=LeafPs(100,400,300, 10,101.3,5000,100);
   

    
      for j=1:4
          
          param.vcopt=50*j;
       
    % light response curve
    
    clear A;
    clear par;
   
    
    for i=1:20
    met.parin=(i-1)*100 ;
    
    [ps]=LeafPs(met.parin,met.co2,met.Tleaf_K, bound_layer_res.co2,met.P_kPa,met.eair_Pa,param.vcopt);
    
    A(i)=ps.aphoto;
    par(i)=met.parin;
    
    wj(i)=ps.wj;
        
    end
    
    figure(1)
    
    plot(par,A)
     ylabel('A {\mu}mol m^{-2}s^{-1}')
     xlabel('PAR {\mu}mol m^{-2}s^{-1}')
     hold on
     legend
    
    
     
      end

      
      
       % co2 response curve
      
      clear A;
      
      for j=1:4
          
          param.vcopt=50*j;
       
    % light response curve
    
    clear A;
            
      met.parin=1000;
      met.Tleaf_K=298;
    
    for i=1:10
    met.co2=(i-1)*100 +50 ;
    
    [ps]=LeafPs(met.parin,met.co2,met.Tleaf_K, bound_layer_res.co2,met.P_kPa,met.eair_Pa, param.vcopt);
    
    A(i)=ps.aphoto;
    co2(i)=met.co2;

% plot wc and wj for different light levels
        
    end
    
    figure(2)
    
    plot(co2,A)
    axis([0,1000,0,40])
     ylabel('A {\mu}mol m^{-2}s^{-1}')
     xlabel('CO2')
     hold on

      end

      % temperature response
      
      clear A;
         
      met.parin=1000;
      met.co2=380;
      
    for j=1:4
          
          param.vcopt=50*j;
    for i=1:25
    met.Tleaf_K=273 + i*2 ;
    
    [ps]=LeafPs(met.parin,met.co2,met.Tleaf_K, bound_layer_res.co2,met.P_kPa,met.eair_Pa, param.vcopt);
    
    A(i)=ps.aphoto;
    Tleaf(i)=met.Tleaf_K;
        
    end
    
    
    figure(3)
   
    plot(Tleaf,A)
    axis([273,325,0,30])
    xlabel('Tleaf')
    ylabel('A {\mu}mol m^{-2}s^{-1}')
    hold on

    end 
