Replace calcrn_vect.m

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);
 
-   albs=(albe-sigmav*albv)./((1-sigmav).^2+sigmav*albv.*albe-(sigmav*albv).^2);
-    i=or(albs>0.4,albs<0.05);
+i=or(albs>0.4,albs<0.05);
 
-    if i==1
-    albs=max(min(albs,0.05),0.4);
-    a=albs*sigmav^2;
-    b=-sigmav*(1+albs*albe);
-    c=albe-albs*(1-sigmav)^2;
+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;
+f = @(x) a*x^2+b*x+c;
 
-    lbv=fzero(f,albe);
-    albv=min(max(albv,0.05),0.3);
-    end
+albv=fzero(f,albe);
+albv=min(max(albv,0.05),0.3);
+end
 
-    albe=sigmav*albv+(albs*(1-sigmav)^2)/(1-sigmav*albs*albv);
+albe=sigmav*albv+(albs*(1-sigmav)^2)/(1-sigmav*albs*albv);
 
 	v1=1-albv*albs.*sigmav;
     v1a=1-albv*albs.*sigmava;
@@ -37,4 +37,4 @@ function [albs,arns,brns,crns,arnv,brnv,crnv,cras,crav]=calcrn_vect(rg,ratm,emis
     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
+end
\ No newline at end of file