pro+sail的计算matlab版本

pro计算叶片透射率反射率； sail用于计算冠层反射率；

function [ refltran ] = fprospect（alphaknN）
%PRO Summary of this function goes here
%   Detailed explanation goes here
% 
    if （k<= 0）
        k=1.d+0;
    else
        inex = s13aaf_leaf（nNk）;
        k=（1.d+0-k）*exp（-k）+k^2*inex;
    end   
    
%N=1 R（1）=rho and T（1）=tau
rho90=frho（90nk）;
tau90=ftau（90nk）;
x=ftav（alphan）/ftav（90n）;
y=x*（ftav（90n）-1）+1-ftav（alphan）;
rho=x*rho90+y;
tau=x*tau90;

%N=N R（N）=refl and T（N）=tran
delta=sqrt（（tau90^2-rho90^2-1）^2-4*rho90^2）;
a=（1+rho90^2-tau90^2+delta）/（2*rho90）;
b=（1-rho90^2+tau90^2+delta）/（2*tau90）;
refl=rho*（a*b^（N-1）-a^-1*b^（1-N））+（tau*tau90-rho*rho90）*（b^（N-1）-b^（1-N））;
tran=tau*（a-a^-1）;
dmt=a*b^（N-1）-a^-1*b^（1-N）-rho90*（b^（N-1）-b^（1-N））;
refl=refl*1.0/dmt;
tran=tran*1.0/dmt;

    
end

function [tau] = ftav（teta n_leaf）

rd=pi/180.d+0;
tetard=teta*rd;
r2=n_leaf^2;
rp=r2+1.d+0;
rm=r2-1.d+0;
ax=（n_leaf+1.d+0）^2/2.d+0;
bx=-（r2-1.d+0）^2/4.d+0;
ds=sin（tetard）;

if（teta==0.d+0）
    tau=4.d+0*n_leaf/（n_leaf+1.d+0）^2;
    return ;
elseif（teta==90.d+0）
    b1=0.d+0;
else
    b1=sqrt（（ds^2-rp/2.d+0）^2+bx）;
end

b2=ds^2-rp/2.d+0;
b0=b1-b2;
ts=（bx^2/（6.d+0*b0^3）+bx/b0-b0/2.d+0）-（bx^2/（6.d+0*ax^3）+bx/ax-ax/2.d+0）;
tp1=-2.d+0*r2*（b0-ax）/rp^2;
tp2=-2.d+0*r2*rp*log（b0/ax）/rm^2;
tp3=r2*（1.d+0/b0-1.d+0/ax）/2.d+0;
tp4=16.d+0*r2^2*（r2^2+1.d+0）*log（（2.d+0*rp*b0-rm^2）/（2.d+0*rp*ax-rm^2））/（rp^3*rm^2）;
tp5=16.d+0*r2^3*（1.d+0/（2.d+0*rp*b0-rm^2）-1.d+0/（2.d+0*rp*ax-rm^2））/rp^3;
tp=tp1+tp2+tp3+tp4+tp5;
tau=（ts+tp）/（2.d+0*ds^2）;

end

function [inex] = s13aaf_leaf（n_leafvaik_leaf）

if （k_leaf>4.d+0） 
    
    if （k_leaf>=85.d+0） 
        inex=0.d+0;
    else
        xx=14.5d+0/（k_leaf+3.25d+0）-1.d+0;
        yy=（（（（（（（（（（（（（（（-1.62806570868460749d-12 ...
            *xx-8.95400579318284288d-13）*xx-4.08352702838151578d-12） ....
            *xx-1.45132988248537498d-11）*xx-8.35086918940757852d-11） ....
            *xx-2.13638678953766289d-10）*xx-1.10302431467069770d-9） ...
            *xx-3.67128915633455484d-9）*xx-1.66980544304104726d-8） ...
            *xx-6.11774386401295125d-8）*xx-2.70306163610271497d-7） ...
            *xx-1.05565006992891261d-6）*xx-4.72090467203711484d-6） ...
            *xx-1.95076375089955937d-5）*xx-9.16450482931221453d-5） ...
            *xx-4.05892130452128677d-4）*xx-2.14213055000334718d-3;
        yy=（（yy*xx-1.06374875116569657d-2）*xx-8.50699154984571871d-2）*xx ...
            +9.23755307807784058d-1;
        inex=exp（-k_leaf）*yy/k_leaf;
    end
else
    xx=0.5d+0*k_leaf-1.d+0;
    yy=（（（（（（（（（（（（（（（-3.60311230482612224d-13*xx+3.46348526554087424d-12）*xx-2.99627399604128973d-11） ...
        *xx+2.57747807106988589d-10）*xx-2.09330568435488303d-9） ...
        *xx+1.59501329936987818d-8）*xx-1.13717900285428895d-7） ...
        *xx+7.55292885309152956d-7）*xx-4.64980751480619431d-6） ...
        *xx+2.63830365675408129d-5）*xx-1.37089870978830576d-4） ...
        *xx+6.47686503728103400d-4）*xx-2

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4237  2013-08-13 20:47  fsail.m
     文件        3557  2013-08-13 20:48  fprospect.m