function [ V ] = windprofiles( vm, rm, r, wp)
%
% This function return radial profiles of azimuthal wind, V, given matrices
% containing the maximum circular wind speed, vm (knots), the radii of maximum
% wind, rm (km), distances r (km) of each event from the points of interest, 
% and the wind profile wp (see below; must be equal to 1, 2 or 3). 
%
% Copyright WindRiskTech, L.L.C., 2011
%-----------------------------------------------------------------------
%
wprofile=wp; % Use holland (1) or emanuel (2) or er2011 (3) wind profile
wprofile=min(wprofile,3); % Forbid use of experimental profiles
%
vm=vm.*1852./3600;  % Convert maximum wind speed to m/s
%
if wprofile == 1
    %
    % Holland 2010 wind profile model
    %
    bs=1.8;
    rn=300/20;     % Holland parameters
    xn=0.8;
    %
    rh=r./rm;
    x=0.5+(rh-1).*(xn-0.5)./(rn-1);
    x=max(x,0.5);
    V=vm.*(rh.^-bs.*exp(1-rh.^-bs)).^x;   
    %
elseif wprofile == 2
    %
    % Outer radius (km)
    %
    r0=1000;
    %
    %  Re-scale radius of maximum winds by random number drawn from log-normal
    %  distribution
    %
    %r0=r0.*rfac;    
    %
    % 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; 
    %
    rat=r./max(rm,1);
    V=vm.*(max((r0-r),0)./(r0-rm)).*sqrt(rat.^mb2.* ...
        (fac1./(nb+mb.*rat.^fac3)+fac2./(1+mb2.*rat.^fac4)));
    %
elseif wprofile == 3     
    crat=1;
    f=5.0e-5;
    f=f.*1000;  % effectively converts radii from kms to meters
    %
    Mm=rm.*vm+0.5.*f.*rm.^2;
    rn=r./rm;
    if crat == 1
        M=Mm.*(2.*rn.^2./(1+rn.^2));
    else    
        M=Mm.*(2.*rn.^2./(2-crat+crat.*rn.^2)).^(1./(2.-crat));
    end
    %V=(M-0.5.*f.*r.^2)./(r+1e-8);
    V=M./(r+1e-8);  % (Add long tail to V to avoid discontinuity in vorticity at outer radius) 3/2013
    V=max(V,0);
%{    
elseif wprofile == 4 % (Experimental; not recommended)
    chi=80;  %  Chi parameter for outer profile
    r0=600;  %  Outer radius (but winds not forced to zero there in this implementation)
    vp=80;   %  Nominal potential intensity. But only chi/vp matters. 
    %
    %r0=rm.*sqrt(chi.*vm./vp).*((16.*vm./(1000.*f.*rm)).^(1/3)-1).^0.75;
    ra=(rm.*r0.^2.*vp./(chi*vm)).^(1/3);
    va=2.*vm.*(ra./rm)./(1+(ra./rm).^2);
    %u=heaviside(ra-r);
    u=max(sign(ra-r),0);
    u(isnan(u))=1; %If u=0, heaviside returns an NaN
    m1=2.*rm.*vm.*(r./rm).^2./(1+(r./rm).^2); % Inner (ER2011) profile
    m2=ra.*va.*(sqrt(1+2.*chi.*ra.*va.*(r-ra).*max((r0.^2-r.^2),0)./(vp.*r0.^4))); % Outer (EC2013) profile
    V=(u.*m1+(1-u).*m2)./max(r,0.1); % Join inner and outer profiles
elseif wprofile == 5  % (Experimental; not recommended)
    V=vm.*0;
    [~,st,~,~]=size(rm);
    for j=1:st
        rmb(1,1)=rm(1,j,1,1);
        vmb(1,1)=vm(1,j,1,1);
        [rnew,Vnew,~,~] = ECwind(rmb,vmb,0,0,2,1e-3,5e-5);
        V(:,j,:,:)=interp1(rnew,Vnew,r(:,j,:,:));
    end    
%}   
end        
%
V=V.*3600./1852;    % Convert wind speed to knots
%
end