From 0409c6db6e1fcfec987b93dc310030f5e145bb81 Mon Sep 17 00:00:00 2001
From: Gilles Boulet <gilles.boulet@cesbio.cnes.fr>
Date: Tue, 6 Jun 2023 04:26:58 +0000
Subject: [PATCH] Update calcrn_vect.m

---
 calcrn_vect.m | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/calcrn_vect.m b/calcrn_vect.m
index a105c9a..44107f8 100644
--- a/calcrn_vect.m
+++ b/calcrn_vect.m
@@ -1,7 +1,22 @@
 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);
+
+    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;
+
+    f = @(x) a*x^2+b*x+c;
+
+    lbv=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;
@@ -22,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
\ No newline at end of file
+end
-- 
GitLab