%NAME: Hapke model
%PURPOSE:Compute the bi-directional reflectance(newer Hapke formula).
%DESCRIPTION:
%The function is coded from equation 12.55 on page 346 in Hapke's book(1993) Reflectance and Emittance Spectroscopy",and made some improvement base paper(2002).
%INPUTS:
%In_para.w-Single scattering albedo.
%In_para.imu-The cosine of incidence angle(s),in radian
%In_para.rmu-The cosine of reflection angle(s),in radian
%In_para.g-Phase angle,in radians.
%In_para.hs-Compaction parameter value(1986 formalism).Also,the width of the Shadow-Hiding Opposition Effect(SHOE).
%In_para.bs0-Empirical factor B0 equation 8.86 on page 228 in Hapke's Book(1993).Also,the hegiht of the Shadow-Hiding Opposition Effect(SHOE)
%In_para.hc-Parameter of the CBOE based on equation 34(Hapke,2002).Al the width of Coherent Backscatter Opposition Effect(CBOE
%In_para.bc0-The amplitude of the CBOE,and is constrained in[0,1](Hapke,
%In_para.stheta-Surface roughness value:a mean slope angle<theta>,(radians) Eq.12.5 in Hapke's book(1993).
%In_para.phf-Phase function type specified by a keyword.
%There are five legal input forms for p:
%1.Single(1-parameter)Henyey-Greenstein phase fuction.Equation 6.7&6.8 on page 108 in Hapke's book(1993).
%2.Double(2-parameter)Henyey-Greenstein phase fuction %Equation 6.18a on %page 121 in Hapke's book(1993),and %on page 529 in Hapke's paper(2002).
%3.3-parameter Henyey-Greenstein phase fuction.Equation(5)in Shepard's paper(2007).
%4.Value,scalar,is 1.0(isotropic)or some else value.Also,an vector same as w dimensionality.
%5.A second-order Legendre polynomial expansion.Equation 6.17 on page 121 in Hapke's book,and equation(1)in Pinty's paper(1989)
%%In_para.p-Parameters of the single particle phase function.
%OUTPUTS:
%Return value is the bi-directional reflectance of the surface[Single,Multiple,BDR].
%OPTIONAL INPUT PARAMETERS:
%N-Legendre expansion order for improvement of multiple scattering accuracy.
%KEYWORD PARAMETERS:
%B_err_fun-If it is set,uses the equation 8.79(p.224 in Hapke's book,1993) which is a error function method(more accurracy),or uses the eq.8.81 %(p.224 in Hapke's book,1993)which is about 3 with error
%H_integral-If it is set,uses the equation 8.57(integral,1%error),or uses the equation 8.55(two-stream,4%error).
%SHOE-The Shadow-Hiding Opposition Effect(SHOE).If it is set,the B_shoe(x)function will be called.(Eq.8.90,Hapke,1993)
%CBOE-The Coherent Backscatter Opposition Effect(CBOE).If it is set,The B_coboe(x)function will be called.(Eq.32,Hapke,2002)
%Multi_imprv-Improves the multiple scattering term for anisotropic scatterers.When the keyword is set,the anisotropic scattering method(Hapke,2002) will be implemented,or,the isotropic scattering method(Hapke,1993) will be implemented.
%pedantic-Flag,if set,returns NaN for w>1.
%MODIFICATION HISTORY:
%2008-Feb-28 Xu Yuanliu created the code that is modified and based on Buie's brdf-code(Marc W.Buie,Lowell Observatory,1997-2004 version)
%2008-May-11 Xu Yuanliu added coherent backscatter oppsoition effect,and improvement in multiple scattering term(Hapke's paper,2002)
%2008-Dec-18 Xu Yuanliu correct the mistake in the code about multiple scattering.
%2009-Feb-27 Xu Yuanliu make the surface roughness as shadow function
%2010-05-29 Li Pingheng translate the code from IDL to MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Main function
function Bi_rfl=hapke(imu,rmu,g,w,b,c,hs,bs0) %,stheta,hc,bc0
phf=4;
N=10;
SHOE=1;
B_err_fun=0;
H_integral=1;
Multi_imprv=1;
CBOE=0;
pedantic=1;
p=[b c];
gamma=sqrt(1.0-w);
if (~ pedantic)
z=find(w > 1);count=length(z);
if count > 0 gamma(z)=-sqrt(w(z)-1.);end
end
% 1-parameter Henyey-Greenstein phase fuction:
if phf==1
Pp=(1.0-p*p)/((1.0+2.0*p*cos(g)+p*p)^1.5);
end
% 2-parameter Henyey-Greenstein phase fuction:P_HG(g)
if phf==2
width=p(:,1)*ones(length(g),1);
fract=p(:,2)*ones(length(g),1);
fwd=(1.0-abs(width)*abs(width))/((1.0+2.0*abs(width)*cos(g)+abs(width)*abs(width))^1.5);
bkwd=(1.0-abs(width)*abs(width))/((1.0-2.0*abs(width)*cos(g)+abs(width)*abs(width))^1.5);
Pp=(1.0-fract)/2.0*bkwd+(1.0+fract)/2.0*fwd;
end
% 3-parameter Henyey-Greenstein phase fuction:
if phf==3
fore=p(:,1)*ones(length(g),1);
back=p(:,2)*ones(length(g),1);
fract=p(:,3)*ones(length(g),1);
fwd=(1.0-fore*fore)/((1.0+2.0*fore*cos(g)+fore*fore)^1.5);
bkwd=(1.0-back*back)/((1.0+2.0*back*cos(g)+back*back)^1.5);
Pp=(1.0-fract)/2.0*bkwd+(1.0+fract)/2.0*fwd;
end
% A second-order Legendre polynomial expansion P_L(g)
if phf==4
Pp=1+p(:,1)*cos(g)+p(:,2)*(3*cos(g)*cos(g)-1)*0.5;
end
Single_r=w*imu/(imu+rmu)*Pp/(4*pi);
if (SHOE)
bscat=zeros(length(g),1);
B_s=bscat;
m=find((g <pi/2)&(bs0 ~= 0));count=length(m);
if (count ~= 0)
B_s(m)=(1.0/hs(m))*(tan(g(m)/2.0));
in_dex=find(B_s == 0);countxx=length(in_dex);
if(countxx ~= 0) B_s(in_dex)=0.00001;end
if (~ B_err_fun)
bscat(m)=1.0/(1.0+B_s(m));%$%equ.8.81,p.224%
else bscat(m)=sqrt((4.0*pi)/B_s(m))*exp(1.0/B_s(m))*(ERF(sqrt(4.0/B_s(m)))-ERF(sqrt(1.0/B_s(m))))+exp(-3.0/B_s(m))-1.0;%%equ.8.79,p.224
end
end
bscat=1.0+bs0*bscat;
Single_r=Single_r*bscat;
end
if (H_integral)
r0=2.0/(1.0+gamma)-1.0;
Hsun=1.0/(1.0-(1.0-gamma)*imu*(r0+(1.0-0.5*r0-r0*imu)*log((1+imu)/imu)));
Hobs=1.0/(1.0-(1.0-gamma)*rmu*(r0+(1.0-0.5*r0-r0*rmu)*log((1+rmu)/rmu)));
else
Hsun=(1+2.0*imu)/(1+2*imu*gamma);
Hobs=(1+2.0*rmu)/(1+2*rmu*gamma);
end
if ((Multi_imprv) | (phf == 4))
Multi_r=w*imu/(imu+rmu)*(Hobs*Hsun-1)/(4*pi);
else
if (phf == 1)
icd=acos(imu);
rfl=acos(rmu);
p_imu=Multi_Pmu_hg1(N,icd,p);
p_rmu=Multi_Pmu_hg1(N,rfl,p);
p_mu=Multi_P_hg1(N,p);
end
if (phf == 2)
icd=acos(imu);
rfl=acos(rmu);
p_imu=Multi_Pmu_hg2(N,icd,p);
p_rmu=Multi_Pmu_hg2(N,rfl,p);
p_mu=Multi_P_hg2(N,p);
end
if (phf == 3)
p_imu=ones(nlen,1);
p_rmu=ones(nlen,1);
p_mu=ones(nlen,1);
end
Multi_imprv=p_imu*(Hobs-1.0)+p_rmu*(Hsun-1.0)+p_mu*(Hobs-1.0)*(Hsun-1.0);
Multi_r=w*imu/(imu+rmu)*Multi_imprv/(4*pi);
end
Bi_r=Single_r+Multi_r;
if (CBOE)
nlen=length(g);
B_cb=zeros(nlen,1);
B_c=B_cb;
m=find((g <pi/2)&(bc0 ~= 0.0));count=length(m);
if (count ~= 0)
B_c(m)=(1.0/hc(m))*(tan(g(m)/2.0));
in_dex=find(B_c == 0.0);countxx=length(in_dex);
if(countxx ~= 0) B_c(in_dex)=0.00001;end
B_cb(m)=(1.0+(1.0-exp(-B_c(m)))/B_c(m))/(2.0*(1.0+B_c(m))*(1.0+B_c(m)));
end
B_cb=1.0+bc0*B_cb;
Bi_r=Bi_r*B_cb;
end
% z=find(stheta ~= 0.0);count=length(z);
% if (count)~= 0
% sfun=Shadow_Hapke(imu,rmu,g,stheta);
% Bi_r(z)=Bi_r(z)*sfun(z);
% end
%Bi_rfl=[Single_r Multi_r Bi_r];
Bi_rfl=Bi_r;
hapke.zip_Hapke模型_hapke_soil
版权申诉
5星 · 超过95%的资源 116 浏览量
2022-07-15
17:09:29
上传
评论 2
收藏 3KB ZIP 举报
局外狗
- 粉丝: 65
- 资源: 1万+
评论3