calcrn_vect.m

function [albs,arns,brns,crns,arnv,brnv,crnv,cras,crav]=calcrn_vect(rg,ratm,emiss,emisv,albe,albv,sigmav,sigmava)

albs=(albe-sigmav*albv)./((1-sigmav).^2+sigmav*albv.*albe-(sigmav*albv).^2);

i=or(albs>0.4,albs<0.05);

if i==1
albs=min(max(albs,0.05),0.4);
a=albs*sigmav^2;
b=-sigmav*(1+albs*albe);
c=albe-albs*(1-sigmav)^2;

f = @(x) a*x^2+b*x+c;

albv=fzero(f,albe);
albv=min(max(albv,0.05),0.3);
end

albe=sigmav*albv+(albs*(1-sigmav)^2)/(1-sigmav*albs*albv);

	v1=1-albv*albs.*sigmav;
    v1a=1-albv*albs.*sigmava;
	v2=1-emisv;
	v3=1-emiss;
	v4=1-sigmav;
    v4a=1-sigmava;
	v5=1-sigmav*v2*v3;
    v5a=1-sigmava*v2*v3;
    
    arns=-(v4a*emiss+emisv*emiss*sigmava)./v5a;
    brns=emisv*emiss*sigmava./v5a;
    crns=(rg*(1-albs).*v4./v1)+(v4a*emiss*ratm./v5a);
    cras=(v4a*emiss*ratm./v5a);
    
    arnv=brns;
    brnv=-sigmava.*(emisv+(emisv*emiss+v4a.*v3*emisv)./v5a);
    crnv=rg*(1-albv)*sigmav.*(1+(albs.*v4)./v1)+sigmava*emisv*ratm.*(1+(v4a.*v3)./v5a);
    crav=sigmava*emisv*ratm.*(1+(v4a.*v3)./v5a);

end