function [ V ] = windprofile( r, rm, r0, vm)
%
% This function return a radial profile of azimuthal wind, V, given a
% radius vector r, a radius of maximum winds rm, and a radius of vanishing
% winds ro. (The latter two are scalars in the same units as r.) There are
% options for several different wind profiles.
%
% Copyright WindRiskTech, L.L.C., January, 2011
%-----------------------------------------------------------------------
%
wprofile=3; % Use holland (1) or emanuel (2) or er2011 (3) wind profile
%
%-----------------------------------------------------------------------
%
if wprofile == 1
    %
    % Holland 2010 wind profile model
    %
    bs=1.8;
    rn=300./20;     % Holland parameters
    xn=0.8;
    %
    q=max(size(r));
    %
    x=zeros(1,q);
    rh=r./rm;
    for i=1:q,
        if rh(i) <= 1
            x(i)=0.5;
        else
            x(i)=0.5+(rh(i)-1).*(xn-0.5)./(rn-1);
        end
    end
    %
    V=vm.*(rh.^-bs.*exp(1-rh.^-bs)).^x;
elseif wprofile == 2
    %
    % Shape parameters
    %
    b=0.25;
    nb=0.9;
    mb=1.6;  
    %
    mb2=2.*mb;
    fac1=(1-b).*(mb+nb);
    fac2=b.*(1+2.*mb);
    fac3=2.*(mb+nb);
    fac4=2.*mb+1; 
    %
    r1=min(r,r0);
    rat=r1./max(rm,1);
    V=vm.*((r0-r1)./(r0-rm)).*sqrt(rat.^mb2.*...
            (fac1./(nb+mb.*rat.^fac3)+fac2./(1+mb2.*rat.^fac4)));
elseif wprofile == 3     
    crat=1;
    f=5.0e-5;
    Mm=rm.*vm+0.5.*f.*rm.^2;
    rn=r./rm;
    M=Mm.*(2.*rn.^2./(2-crat+crat.*rn.^2)).^(1./(2.-crat));
    V=(M-0.5.*f.*r.^2)./(r+1e-8);
    V=max(V,0);
end        
end

