function [LI, GCO2_kg_m2_hr, NCO2_kg_m2_hr]=Photosynthesis(CO2,Tl,fAPAR, TYPE)
%Goudriaan and van Laar, 1994. Listing 8.2, page 164
%INPUT: CO2 in ppm
%   Tl: Temperature of leaf in degree Celcius
%   fAPAR fraction of LI that is absorbed as PAR
%   Type: either 'C3' or 'C4' plants
%OUTPUT:
%   LI: light intercepted, for plotting on x-axis. Can also be PAR (then
%   fAPAR should be 1.0)
%   GCO2_kg_m2_hr: gross CO assimilation in kg/ha/hr
%   NCO2_kg_m2_hr: net CO assimilation in kg/ha/hr
    
    if nargin<4
        TYPE='C3';
    end
    if  strcmp(TYPE,'C3')
        GAMM25=50;
        EFF0=10;%17 C3 plant  ug CO2 J-1
        FMMT=[0 5   30   40 50;...
              0 0 1200 1200  0 ];%Max assimilation, ug CO2/m2/s, 1200=43.2 kgCO2ha/hr
        GMT=[ 0 2    25 39 50;...
              0 0 0.005  0  0];
    else
        GAMM25=5;
        EFF0=14;%14 C4 plant, 11 C3 plant  ug CO2 J-1    
        FMMT=[0 8   30   40 50;...
              0 0 2000 2000  0 ];%Max assimilation, ug CO2/m2/s
        GMT=[ 0 7   35 45 50;...
              0 0 0.02  0  0];
    end
    
    RD20=50;
    RS=100; RBV= 20; %Stomatic resistance param
    %Intercepted light (can be either radiation or PAR, will be returned for plotting purposes)
    LI=(1:600)';                                    %J/s/m2
    APAR= LI .* fAPAR;                              %J PAR/s/m2
    
    FMM=interp1(FMMT(1,:),FMMT(2,:),Tl);            %interpolate using Tl and FMMT
    GM=interp1(GMT(1,:),GMT(2,:),Tl);               %interpolate using Tl and GMT
    
    RD=RD20 .* exp(0.07*(Tl-20));                   %actual respiration
    GAMM= GAMM25.* exp(0.07 * (Tl -25));            %CO2 compensation point
    EFF=EFF0.*(CO2-GAMM)/(CO2 + 2.*GAMM);           %Net light conversion factor

    RM=1./GM;                                       %mesophyll resistance to CO2 
    GASLAW=(44/24).*(293./(273+Tl))*1000;           %Conversion factor 44 (g/mol) / 24 (mol CO2 / m3) to T to gC/l absorbed 
    FNC=(CO2-GAMM).*GASLAW./(RM + 1.6 * RS+1.6 * RBV);% Max CO2 by diffusion, ug CO2/m2/s
    
    FNMAX=min(FMM, FNC);                            %Net CO2 assimilation, ug CO2/m2/s
    FGMAX=FNMAX+RD;
    Fg=max(0,FGMAX.*(1 - exp(-EFF.*APAR./FGMAX)));  %Gross CO2 assimilation, ug CO2/m2/s
    Fn = Fg - RD;                                   %Net CO2 assimilation, ug CO2/m2/s
    
    convfact= 3600 .* 10000 .* 1E-9;                %s/hr * m2/ha * kg / ug
    GCO2_kg_m2_hr=Fg .* convfact;                   %kg CO2/ha/hr = ug CO2/m2/s * s/hr * m2/ha * kg / ug
    NCO2_kg_m2_hr=Fn .* convfact;                   %kg CO2/ha/hr = ug CO2/m2/s * s/hr * m2/ha * kg / ug
end