Upload New File

evalSAS.m

+%
+% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
+%
+% Author: Nicolas Rodriguez (nicolas.bjorn.rodriguez@gmail.com)
+%
+function [p,omega]=evalSAS(dT,ps,S,evt_frac,type,varargin)
+% calculates the SAS functions and TTDs based on the RTD ps and SAS
+% function types and parameters
+
+% Omega: cumulative SAS function Om(Ps,t) or Om(S_T,t)
+% omega: absolute SAS function w(T,t)
+% p: TTD pQ(T,t) (pdf)
+% ps: RTD pS(T,t) (pdf)
+
+    switch type
+        case 'RS' % random sampling
+            omega=ones(size(ps));
+            p=ps;
+        case 'unif' % uniform SAS on a given range of S_T
+            S_T=S*cumsum(ps*dT);
+            S_Tmax=varargin{1};
+            Omega=(S_T<=S_Tmax).*S_T./S_Tmax+(S_T>S_Tmax);
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'pl' % power law
+            k=varargin{1};
+            Omega=cumsum(ps*dT).^k;
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'tvpl' % time-varying power law (parameterized with storage)
+            % k2>k1: wetter => lower k
+            k1=varargin{1};
+            k2=varargin{2};
+            Smin=varargin{3};
+            Smax=varargin{4};
+            wetness=max(min((S-Smin)/(Smax-Smin),1),0);
+            k=k1+(1-wetness)*(k2-k1);
+            Omega=cumsum(ps*dT).^k;
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+         case 'tvpl_Q' % time-varying power law (parameterized with discharge)
+            % k2>k1: wetter => lower k
+            k1=varargin{1};
+            k2=varargin{2};
+            Qmin=varargin{3}; % always 0
+            Qmax=varargin{4};
+            Q=varargin{5}; % Q(t)
+            wetness=max(min((Q-Qmin)/(Qmax-Qmin),1),0);
+            k=k1+(1-wetness)*(k2-k1);
+            Omega=cumsum(ps*dT).^k;
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'tvpl_dS' % time-varying power law (parameterized with storage variations)
+            % k2>k1: wetter => lower k
+            k1=varargin{1};
+            k2=varargin{2};
+            k=k1+(1-evt_frac)*(k2-k1);
+            Omega=cumsum(ps*dT).^k;
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'tvpl_S_dS' % time-varying power law (parameterized with storage and storage variations)
+            % k2>k1: events => lower k
+            % k2max > k2min: wetter => lower k2 (ISE)
+            k1=varargin{1};
+            k2min=varargin{2};
+            k2max=varargin{3};
+            Smin=varargin{4};
+            Smax=varargin{5};
+            wetness=max(min((S-Smin)/(Smax-Smin),1),0);
+            k2=k2min+(1-wetness)*(k2max-k2min);
+            k=k1+(1-evt_frac)*(k2-k1);
+            Omega=cumsum(ps*dT).^k;
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'bet' % beta distribution
+            a=varargin{1};
+            b=varargin{2};
+            Omega=betacdf(cumsum(ps*dT),a,b);
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'tvb' % time-varying beta distribution
+            amax=varargin{1};
+            afrac=varargin{2};
+            Smin=varargin{3};
+            Smax=varargin{4};
+            wetness=max(min((S-Smin)/(Smax-Smin),1),0);
+            a=amax*(1-afrac*wetness);
+            if a<=1
+                Omega=betacdf(cumsum(ps*dT),a,1);
+            elseif a>1
+                Omega=betacdf(cumsum(ps*dT),1,1-(a-1));
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'gam' % gamma distribution
+            S_T=varargin{10};
+            G1=varargin{11};
+            k=varargin{1};
+            theta=varargin{2};
+            if ~isempty(S_T)
+                Omega=interp1(S_T,G1,S*cumsum(ps*dT));
+            else
+                Omega=gamcdf(S*cumsum(ps*dT),k,theta);
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'tvg' % time-varying gamma distribution
+            k=varargin{1};
+            lambda=varargin{2};
+            dSc=varargin{3};
+            dS=varargin{4};
+            theta=max(lambda*(dS-dSc),1);
+            Omega=gamcdf(S*cumsum(ps*dT),k,theta);
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case '2g' % composite double gamma
+            S_T=varargin{10};
+            G1=varargin{11};
+            G2=varargin{12};
+            lambda_1=varargin{1};
+            lambda_2=varargin{2};
+            k1=varargin{3};
+            k2=varargin{4};
+            theta1=varargin{5};
+            theta2=varargin{6};
+            if ~isempty(S_T)
+                % interpolate the pre-calculated gamma distributions to the required values of S_T
+                Omega=lambda_1*interp1(S_T,G1,S*cumsum(ps*dT))+lambda_2*interp1(S_T,G2,S*cumsum(ps*dT));
+            else
+                Omega=lambda_1*gamcdf(S*cumsum(ps*dT),k1,theta1)+lambda_2*gamcdf(S*cumsum(ps*dT),k2,theta2);
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case '2norm' % composite double normal distribution
+            S_T=varargin{10};
+            lambda_1=varargin{1};
+            lambda_2=varargin{2};
+            mu1=varargin{3};
+            mu2=varargin{4};
+            theta1=varargin{5};
+            theta2=varargin{6};
+            Omega=lambda_1*normcdf(S*cumsum(ps*dT),mu1,theta1)+lambda_2*normcdf(S*cumsum(ps*dT),mu2,theta2);
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case '3g' % composite triple gamma distribution
+            S_T=varargin{10};
+            G1=varargin{11};
+            G2=varargin{12};
+            G3=varargin{13};
+            lambda_1=varargin{1};
+            lambda_2=varargin{2};
+            lambda_3=varargin{3};
+            k1=varargin{4};
+            k2=varargin{5};
+            k3=varargin{6};
+            theta1=varargin{7};
+            theta2=varargin{8};
+            theta3=varargin{9};
+            if ~isempty(S_T)
+                Omega=lambda_1*interp1(S_T,G1,S*cumsum(ps*dT))+lambda_2*interp1(S_T,G2,S*cumsum(ps*dT))+lambda_3*interp1(S_T,G3,S*cumsum(ps*dT));
+            else
+                Omega=lambda_1*gamcdf(S*cumsum(ps*dT),k1,theta1)+lambda_2*gamcdf(S*cumsum(ps*dT),k2,theta2)+lambda_3*gamcdf(S*cumsum(ps*dT),k3,theta3);
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+        case 'unif-gam' % composite unif+gamma distribution
+            S_T=varargin{10};
+            G1=varargin{11};
+            lambda_1=varargin{1}*evt_frac;
+            lambda_2=1-lambda_1;
+            S1=varargin{3};
+            k=varargin{4};
+            theta=varargin{5};
+            if ~isempty(S_T)
+                Omega=lambda_1*((S*cumsum(ps*dT)<=S1).*S.*cumsum(ps*dT)./S1+(S*cumsum(ps*dT)>S1))+lambda_2*interp1(S_T,G1,S*cumsum(ps*dT));
+            else
+                Omega=lambda_1*((S*cumsum(ps*dT)<=S1).*S.*cumsum(ps*dT)./S1+(S*cumsum(ps*dT)>S1))+lambda_2*gamcdf(S*cumsum(ps*dT),k,theta);
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+         case 'unif-2gam' % composite unif+2gamma distribution
+            S_T=varargin{10};
+            G1=varargin{11};
+            G2=varargin{12};
+            lambda_1=varargin{1}*evt_frac;
+            lambda_2=varargin{2};
+            lambda_3=1-lambda_2-lambda_1;
+            S1=varargin{4};
+            k1=varargin{5};
+            theta1=varargin{6};
+            k2=varargin{7};
+            theta2=varargin{8};
+            if ~isempty(S_T)
+                Omega=lambda_1*((S*cumsum(ps*dT)<=S1).*S.*cumsum(ps*dT)./S1+(S*cumsum(ps*dT)>S1))+lambda_2*interp1(S_T,G1,S*cumsum(ps*dT))+lambda_3*interp1(S_T,G2,S*cumsum(ps*dT));
+            else
+                Omega=lambda_1*((S*cumsum(ps*dT)<=S1).*S.*cumsum(ps*dT)./S1+(S*cumsum(ps*dT)>S1))+lambda_2*gamcdf(S*cumsum(ps*dT),k1,theta1)+lambda_3*gamcdf(S*cumsum(ps*dT),k2,theta2);
+            end
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+         case 'unif-2bet' % composite unif+2beta distribution
+            lambda_1=varargin{1}*evt_frac;
+            lambda_2=varargin{2};
+            lambda_3=1-lambda_2-lambda_1;
+            u1=varargin{4};
+            a1=varargin{5};
+            b1=varargin{6};
+            a2=varargin{7};
+            b2=varargin{8};
+            Ps=cumsum(ps*dT);
+            Omega=lambda_1*((Ps<=u1).*Ps./u1+(Ps>u1))+lambda_2*betacdf(Ps,a1,b1)+lambda_3*betacdf(Ps,a2,b2);
+            p=diff([0 Omega])/dT;
+            omega=(ps~=0).*p./ps;
+            omega(ps==0)=NaN;
+    end
+end
\ No newline at end of file