function [pstorm,rw] = psurf(vstorm,rmstorm,latstorm,wprofile,varargin)
%
%  This function uses gradient wind balance to estimate the distribution of
%  surface pressure with radius for single storms with maximum wind vstorm
%  (knots), radius of maximum winds rmstorm (km), latitude latstorm
%  (degrees), and desired wind profile wprofile. It accepts as optional
%  arguments, in order, secondary wind maximum (knots) and radius of
%  secondary wind maximum (km). It returns surface pressure pstorm (hPa)
%  and the radius vector it uses, rw (km).
%  July, 2022
%  Copyright WindRiskTech, LLC 2022
%--------------------------------------------------------------------------
rw=0:0.5:1000;    % Specification of radius vector (km)
p0=1013;          % Ambient surface pressure (hPa)
reduc=0.84;       % Factor to reduce gradient wind by to get surface wind 
Ts=300;           % Assumed surface air temperature (K)
%--------------------------------------------------------------------------
vstorm=vstorm/reduc;   %  Find gradient wind from 10 m wind
qrw=max(size(rw));
f=(pi/(6*3600))*sin(pi*abs(latstorm)/180); % Coriolis
dr=1000*(rw(2)-rw(1));
rw=max(rw,1e-6);
rwi=0.001./rw;
%
if nargin > 4
    vse=varargin{1};
    rmse=varargin{2};
    vse=vse/reduc;
    vcirc=windprofilem(vstorm,rmstorm,vse,rmse,rw,wprofile);
else
    vcirc=windprofiles(vstorm,rmstorm,rw,wprofile);
end 
vcirc=vcirc*1852/3600; % Convert to m/s
%
% Integrate gradient wind equation inward
%
ptemp=zeros(1,qrw);
for i=qrw-1:-1:1   
    ptemp(i)=ptemp(i+1)-0.5*dr*(vcirc(i).^2.*rwi(i)+f*vcirc(i)+ ...
        vcirc(i+1).^2.*rwi(i+1)+f*vcirc(i+1));
end    
pstorm=p0.*exp(ptemp./(287*Ts));
%
end