text:
function [tq,q,tu,u,t,p,Q,code]=hydrograph(uhname,rdname,method,...
                                           K,D,R,Iar,A,CN,Tc,units)
% SCS Unit Hydrograph Convolution (R11)
%
% [tq,q,tu,u,t,p,Q,code]=hydrograph(uhname,rdname,method,...
%                                   K,D,R,Iar,A,CN,Tc,units)
%
% Input
%   uhname  dimensionless unit hydrograph filename
%   rdname  dimensionless rainfall distribution filename
%   method  interpolation method
%   K       peak factor
%   D       storm duration (hr)
%   R       rainfall depth (mm | in)
%   Iar     initial abstraction ratio
%   A       basin area (km^2 | ac)
%   CN      curve number
%   Tc      time of concentration (min)
%   units   units code
%             0 = imperial
%             1 = metric
%
% Output
%   tq      runoff hydrograph time (hr)
%   q       runoff hydrograph flow rate (cms | cfs)
%   tu      unit hydrograph time (hr)
%   u       unit hydrograph flow rate (cms | cfs)
%   t       rainfall-runoff time (hr)
%   p       cumulative rainfall (mm | in)
%   Q       cumulative runoff (mm | in)
%   code    return code
%             0 = fail
%             1 = pass
%
% Example:
%
%  uhname='library\gamma.duh'; rdname='library\scsii-024.drd'; K=484;
%  method='linear'; D=24; R=200; Iar=0.2; A=3; CN=75; Tc=90; units=1;
%  [tq,q,tu,u,t,p,Q,code]=hydrograph(uhname,rdname,method,...
%                                    K,D,R,Iar,A,CN,Tc,units)
%
% See also help\hydographui.html.

% Version 2.03 Copyright(c)2008
% Tom Davis (tdavis@metzgerwillard.com)
%
% Last revision: 03/16/2008

if units
  R=R/25.4;                             % rainfall depth (in)
  cm2cf=35.3146667214886;               % cubic meters to cubic feet
  %    =1e6/(12^3*2.54^3);
  sk2ac=247.105381467165;               % square kilometers to acres
  %    =1e10/(43560*12^2*2.54^2);
  A=A*sk2ac;                            % basin area (ac)
end
A =A/640;                               % basin area (mi^2)
Tc=Tc/60;                               % time of concentration (hr)
d =Tc/7.5;                              % rain burst duration (hr)
Tp=5*d;                                 % time to peak (hr)
method1=lower(method);
method2=method1;
uhname=lower(uhname);
code=1;

if ~isempty(findstr('triangle',uhname))
  uhname='triangle';
elseif ~isempty(findstr('gamma',uhname))
  uhname='gamma';
end

switch uhname                           % dimensionless unit hydrograph
  case 'triangle'
    Tb=3872/(3*K);                      % time base
    % =2(5280^2/60^2/12)/K
    uh=[0,0;1,1;Tb,0];                  % triangular distribution
    method1='linear';
  case 'gamma'
    f =K*3/1936;
    % =K/(5280^2/60^2/12)
    a0=0.045+0.5*f+5.6*f^2+0.3*f^3;     % cubic estimate
    options=optimset('display','off','tolx',eps);
    g =inline('gamma(a)*exp(a)/a^a-1/f','a','f');
    a =fzero(g,a0,options,f);
    tu=(0:0.2:round(2500/K))';
    u =(tu.*exp(1-tu)).^a;              % gamma distribution
    u(end)=0;                           % force zero ordinate
  otherwise
    uh=load(uhname);                    % tabular distribution
    [m,n]=size(uh);
    [umax,ndx]=max(uh(:,2));            %#ok umax not used
    umin=min(min(uh));
    if m<3 | n~=2 | umin<0 | uh(1,:)~=[0 0] | uh(ndx,:)~=[1 1] %#ok (R11)
      msgbox('Unit Hydrograph is improperly formed.',...
        'File Error','error','modal');
      tq=[];q=[];tu=[];u=[];t=[];p=[];Q=[];code=0;
      return
    end
end

if ~strcmp(uhname,'gamma')
  uh=[uh;uh(end,1)+0.2,0];              % force zero ordinate
  Tu=uh(:,1);
  U =uh(:,2);
  % resample unit hydrograph
  tu=(0:0.2:Tu(end))';
  u =interp1(Tu,U,tu,method1);
end

tu=Tp*tu;                               % unit hydrograph time (hr)
u =K*A*u/Tp;                            % unit hydrograph flow rate (cfs)
u(u<0)=0;

rd=load(rdname);                        % dimensionless rainfall distribution
[m,n]=size(rd);
rmax=max(max(rd));
rmin=min(min(rd));
if m<2 | n~=2 | rmin<0 | rmax>1 | rd(1,:)~=[0 0] | rd(end,:)~=[1 1] %#ok (R11)
  msgbox('Rainfall Distribution is improperly formed.',...
    'File Error','error','modal');
  tq=[];q=[];tu=[];u=[];t=[];p=[];Q=[];code=0;
  return
end

T =D*rd(:,1);
P =R*rd(:,2);
% resample rainfall distribution
t =(0:d:T(end))';                       % rainfall time (hr)
p =interp1(T,P,t,method2);              % rainfall depth (in)

Q =zeros(size(t));
s =1000/CN-10;                          % retention (in)
i =find(p>Iar*s);
Q(i)=(p(i)-Iar*s).^2./(p(i)+(1-Iar)*s); % runoff depth (in)
dQ=[diff(Q);0];                         % incremental runoff depth (in)
q =conv(dQ,u);                          % runoff flow rate (cfs)
q(q<0)=0;
tq=(0:d:t(end)+tu(end))';               % runoff time (hr)

if units
  u=u/cm2cf; q=q/cm2cf;                 % cms
  p=p*25.4; Q=Q*25.4;                   % mm
end


text:
function [ dPriors ] = log_gmm_div_prior(W,Priors)
%LOG_GMM_DIV_PRIOR
%
%   input -----------------------------------------------------------------
%
%       o X         : (D x N),      Samples
%
%       o W         : (N x K),      Responsibility factor
%
%       o Priors    : (1 x K),      Weights
%
%   output ----------------------------------------------------------------
%
%       o dPriors   : (1 x K),      Derivative of weights
%
%
%

K       = length(Priors);
dPriors = zeros(1,K);

for k=1:K
                    %(N x 1)
    dPriors(k) = mean(W(:,k) ./ Priors(k));

end

 dPriors =  dPriors./sum(dPriors);

end

author: epfl-lasa


text:
function prob_test011 ( )

%*****************************************************************************80
%
%% TEST011 tests BETA, GAMMA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST011\n' );
  fprintf ( 1, '  BETA evaluates the Beta function;\n' );
  fprintf ( 1, '  GAMMA evaluates the Gamma function.\n' );

  a = 2.2;
  b = 3.7;

  beta1 = beta ( a, b );
  beta2 = gamma ( a ) * gamma ( b ) / gamma ( a + b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Argument A =                   %14f\n', a );
  fprintf ( 1, '  Argument B =                   %14f\n', b );
  fprintf ( 1, '  Beta(A,B) =                    %14f\n', beta1 );
  fprintf ( 1, '  (Expected value = 0.0454 )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Gamma(A)*Gamma(B)/Gamma(A+B) = %14f\n', beta2 );

  return
end

author: johannesgerer

repo: jburkardt-m