rm = 4e4;   % radius of max winds
vm = 40.0;  % max gradient wind speed
b = 1.3;    % Holland (1980) b parameter
lat = 15;   % latitude
rho = 1.1;  % air density
dp = vm^2 * rho * 2.71828 / b;  % environment-centre pressure difference.
pe = 1e5;   % Environment pressure
f = 2*7.292e-5*sin(lat*pi/180);  % Coriolis parameter
Ut = -5;    % Translation speed
ds = 5e3;   % Grid spacing
lp = 30;    % Grid size is 2*lp + 1

[x,y] = meshgrid( [-lp:lp]*ds, [-lp:lp]*ds );
lam = atan2(y,x);
r = sqrt(x.^2 + y.^2); r(find(r==0)) = nan;

% other profiles here...
%beta = [vm rm 50000 92200]
% end...
xkm = x*1e-3; ykm = y*1e-3; xx = xkm(1,:); yy = ykm(:,1);
cx = rm*1e-3*sin(0:0.01:6.28);  % Circle for rmw plots.
cy = rm*1e-3*cos(0:0.01:6.28);

% Define parametric gradient wind (V) and vorticity (Z) profiles
% using Holland (1980) bogus with core modified to quadratic.
rmrb = (rm./r).^b;
V = sqrt( b*dp/rho * exp(-rmrb) .* rmrb + (0.5*f*r).^2 ) - 0.5*f*r;
Z = -f + sqrt( 4*b*dp/rho * exp(-rmrb) .* rmrb + (f*r).^2 ) ...
       .* ( b*dp/rho * rmrb .* (2 + b*(rmrb - 1)) .* exp(-rmrb) + (f*r).^2 )...
       ./ ( r .* ( 4*b*dp/rho * exp(-rmrb) .* rmrb + (f*r).^2 ) );

if 1 % Quadratic core
   Vm = max(V(:));
   rs = r( find( V==Vm ) ); rs = rs(1);
   rmrsb = (rm/rs)^b;
   Vs = sqrt( b*dp/rho * exp(-rmrsb) * rmrsb + (0.5*f*rs)^2 ) - 0.5*f*rs;
   icore = find(r < rs);
   V(icore) = Vs * (r(icore)/rs) .* (2 - (r(icore)/rs));
   Z(icore) = Vs/rs * (4 - 3*r(icore)/rs);
   'FIRST CASE*************'
else % Fit cubic at rm, matching derivatives
   E = exp(1);
   Vm = sqrt( b*dp/(rho*E) + (0.5*f*rm)^2 ) - 0.5*f*rm;
   dVm = -f/2 + ( E*(f^2)*rm*sqrt((4*b*dp/rho)/E + (f*rm)^2) )/ ...
                ( 2*(4*b*dp/rho + E*(f*rm)^2) );
   d2Vm = ( b*dp*(-4*b^3*dp/rho - (-2 + b^2)*E*(f*rm)^2) ) / ...
          ( E*rho*sqrt((4*b*dp)/(E*rho) + (f*rm)^2) * ...
            (4*b*dp*rm^2/rho + E*(f*rm^2)^2) );
   aa = ( d2Vm/2 - ( dVm - Vm/rm )/rm ) / rm;
   bb = ( d2Vm - 6*aa*rm ) / 2;
   cc = dVm - 3*aa*rm^2 - 2*bb*rm;
   icore = find(r < rm);
 % xx = r .* ( r .* ( r*aa + bb ) + cc );
   V(icore) = r(icore) .* ( r(icore) .* ( r(icore)*aa + bb ) + cc );
   Z(icore) = r(icore) .* ( r(icore) * 4*aa + 3*bb ) + 2*cc;
   pm = pe - dp*(1 - exp(-1));
   p(icore) = pm + ( r(icore).^6 - rm^6 ) * (aa^2)/6 ...
                 + ( r(icore).^5 - rm^5 ) * (2*aa*bb)/5 ...
                 + ( r(icore).^4 - rm^4 ) * (aa*(2*cc+f)+bb^2)/4 ...
                 + ( r(icore).^3 - rm^3 ) * (bb*(2*cc+f))/3 ...
                 + ( r(icore).^2 - rm^2 ) * (cc*(cc+f))/2;
   'ELSE CASE*********'
end
Vm
%use the other profile; ie equate V to vw...
%[kp,vw]=willoughby_pv(beta,r,f)
%V=vw
%
Vx = -V.*y./r;  % Cartesian winds
Vy =  V.*x./r;

I2 = (f + 2*V./r).*(f + Z);  % Inertial stability squared

K = 50;    % Diffusivity
C = 0.002; % Drag coefficient

% Calculate alpha, beta, gamma, chi, eta and psi in Kepert (2001). 
% The III's are to get the solution in the case where V/r > I.
al = (2 * V ./ r + f) / (2*K);
be = (f + Z) / (2*K);
gam = V ./ (2*K*r);
III = find(gam > sqrt(al.*be));
fixup = find(isnan(al) | isnan(be) | isnan(gam));
chi = (C/K)*V ./ sqrt(sqrt(al.*be)      ); chi( fixup ) = nan;
eta = (C/K)*V ./ sqrt(sqrt(al.*be) + gam); eta( fixup ) = nan;
psi = (C/K)*V ./ sqrt( abs( sqrt(al.*be) - gam) ); psi( fixup ) = nan;
albe = sqrt(al./be);

% Calculate A_k's ... p and m refer to k = +1 and -1. i is sqrt(-1).
A0 =  -chi.*V.*(1 + i*(1 + chi)) ./ ( 2*chi.^2 + 3*chi + 2 ); A0(fixup) = nan;
u0s = albe .* real(A0);   % Symmetric surface wind component
v0s =         imag(A0);
Am = -( (1 + (1+i)*eta)./albe ...
      + (2 + (1+i)*eta) ).*psi.*Ut ...
   ./ ( (2 + 2*i)*(1 + eta.*psi) + 3*psi + 3*i*eta ); 
AmIII = -( (1 + (1+i)*eta)./albe ...
           + (2 + (1+i)*eta) ).*psi.*Ut ...
       ./ ( 2 - 2*i + 3*(psi + eta) + (2 + 2*i)*eta.*psi );
Am(III) = AmIII(III);
Am(fixup) = nan;
ums = albe .* real( Am .* exp(-i*lam) );  % First asymmetric surface component
vms =         imag( Am .* exp(-i*lam) );
Ap = -( (1 + (1+i)*psi)./albe ...
      - (2 + (1+i)*psi) ).*eta.*Ut ...
   ./ ( (2 + 2*i)*(1 + eta.*psi) + 3*eta + 3*i*psi ); 
ApIII = -(  (1 + (1-i)*psi)./albe ...
           - (2 + (1-i)*psi) ).*eta.*Ut ...
       ./ ( 2 + 2*i + 3*(eta + psi) + (2 - 2*i)*eta.*psi );
Ap(III) = ApIII(III);
Ap(fixup) = nan;
ups = albe .* real( Ap .* exp(i*lam) ); % Second asymmetric surface component
vps =         imag( Ap .* exp(i*lam) );

% Total surface wind in (moving coordinate system)
us =     u0s + ups + ums;
vs = V + v0s + vps + vms;

% Total surface wind in stationary coordinate system (f for fixed)
usf = us + Ut*cos(lam);
vsf = vs - Ut*sin(lam);
Uf =   + Ut*cos(lam);
Vf = V - Ut*sin(lam);

%figure(1); clf;
%subplot(2,2,1); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%[cs,h] = contourf(xkm,ykm,us,[-20:2:10],'k-'); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm,us,[-20:2:10],'k-'); clabel(cs,h,[-20:4:10]);
%plot(cx,cy,'w');
%xlabel('Storm-relative surface u');
%title(['V_m=' num2str(vm) ' r_m=' num2str(rm*1e-3) ' b=' num2str(b) ...
%      ' lat=' num2str(lat) ' K=' num2str(K) ' C=' num2str(C) ]);
%subplot(2,2,2); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%[cs,h] = contourf(xkm,ykm,vs,[0:5:50],'k-'); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm,vs,[0:5:50],'k-'); clabel(cs,h,[0:10:50]);
%plot(cx,cy,'w');
%xlabel('Storm-relative surface v');
%subplot(2,2,3); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%[cs,h] = contourf(xkm,ykm,usf,[-20:2:10],'k-'); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm,usf,[-20:2:10],'k-'); clabel(cs,h,[-20:4:10]);
%plot(cx,cy,'w');
%xlabel('Earth-relative surface u');
%subplot(2,2,4); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%[cs,h] = contourf(xkm,ykm,vsf,[0:5:50],'k-'); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm,vsf,[0:5:50],'k-'); clabel(cs,h,[0:10:50]);
%plot(cx,cy,'w');
%xlabel('Earth-relative surface v');

swf1 = sqrt( usf.^2 + vsf.^2 ) ;
% Wind speed (norm) for a moving coordinate system
swf2 = sqrt( us.^2 + vs.^2 ) ;

figure(1); clf;
subplot(1,1,1); cla; hold on; set(gca,'DataAspectRatio',[2 2 2]); xlim([-150,150]); ylim([-150,150]);
[cs,h] = contourf(xkm,ykm, swf1); set(h,'LineStyle','none');
[cs,h] = contour(xkm,ykm, swf1); clabel(cs,h);
plot(cx,cy,'w');
xlabel(' Surface Earth Relative Wind Speed in m/s ');
%
print  andrew_at_ldfl_150_case2
%
%figure(2); clf;
%subplot(1,1,1); cla; hold on; set(gca,'DataAspectRatio',[2 2 2]); xlim([-150,150]); ylim([-150,150]);
%[cs,h] = contourf(xkm,ykm, swf2); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm, swf2); clabel(cs,h);
%plot(cx,cy,'w');
%xlabel(' Surface Wind Speed in storm relative coordinate (m/s) ');
%
print  andrew_at_ldfl_50_case2
%
% write data to a file

fid = fopen('storm_xkm.txt','w');
fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',xkm);
fclose(fid)
fid = fopen('storm_ykm.txt','w');
fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',ykm);
fclose(fid)
fid = fopen('storm_us.txt','w');
fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',us);
fclose(fid)
fid = fopen('storm_vs.txt','w');
fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',vs);
fclose(fid)
%us(1,3)
%us(1,11)
%us(11,1)
%us(3,1)
%swf2(1,1)
%swf2(2,3)
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% +++++++++++++++++++++3 d part now ++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% Coefficients p_k for k = 0, 1, -1.
p0 = -(1 + i) * (al .* be).^(1/4);
pp = -(1 + i) * sqrt( sqrt(al .* be) + gam );
pm = -(1 + i) * sqrt( abs( sqrt(al .* be) - gam ) );
pm(III) = -i * pm(III);

% Set up a 3-d grid and calculate winds and horizontal divergence,
% then vertical velocity.
nz = 5; z = [1 : nz]*10;

[ny,nx] = size(x);
nz = max(size(z));
u = zeros(nz,ny,nx);
v = zeros(nz,ny,nx);
ux = zeros(nz,ny,nx);  % x for cartesian wind
vx = zeros(nz,ny,nx);
w = zeros(nz,ny,nx);
hdiv = zeros(nz,ny,nx);
uxa = zeros(nz,ny,nx);  % Asymmetric parts have trailing a
vxa = zeros(nz,ny,nx);
wa = zeros(nz,ny,nx);
hdiva = zeros(nz,ny,nx);

for (kk = 1:nz),
   % These w's are really omega's!
   wm = Am .* exp(z(kk)*pm - i*lam); wm(fixup) = nan;
   w0 = A0 .* exp(z(kk)*p0        ); w0(fixup) = nan;
   wp = Ap .* exp(z(kk)*pp + i*lam); wp(fixup) = nan;

   u(kk,:,:) =     albe.*real(wm + w0 + wp);
   v(kk,:,:) = V +       imag(wm + w0 + wp);
   ua(kk,:,:) = albe.*real(wm + wp);
   va(kk,:,:) =       imag(wm + wp);

   ux(kk,:,:) = ( squeeze(u(kk,:,:)).*x - squeeze(v(kk,:,:)).*y ) ./ r;
   vx(kk,:,:) = ( squeeze(v(kk,:,:)).*x + squeeze(u(kk,:,:)).*y ) ./ r;
   hdiv(kk,:,:) = div(squeeze(ux(kk,:,:)), squeeze(vx(kk,:,:)), ds);
   uxa(kk,:,:) = ( squeeze(ua(kk,:,:)).*x - squeeze(va(kk,:,:)).*y ) ./ r;
   vxa(kk,:,:) = ( squeeze(va(kk,:,:)).*x + squeeze(ua(kk,:,:)).*y ) ./ r;
   hdiva(kk,:,:) = div(squeeze(uxa(kk,:,:)), squeeze(vxa(kk,:,:)), ds);
end;

% Integrate the hydrostatic equation - these w's are vertical velocity.
w (1,:,:) = -hdiv (1,:,:) * 10.0;
wa(1,:,:) = -hdiva(1,:,:) * 10.0;
for (kk = 2:nz),
  w (kk,:,:) = w (kk-1,:,:) - hdiv (kk,:,:) * (z(kk) - z(kk-1));
  wa(kk,:,:) = wa(kk-1,:,:) - hdiva(kk,:,:) * (z(kk) - z(kk-1));
end;

ucon = -20:2:10;
vcon = 0:5:40;
wcon = -0.2:0.02:0.5;

%

figure(3); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(u(:,lp+1,:)),ucon); hold on;
[cs,h] = contour(xx,z,squeeze(u(:,lp+1,:)),[0 0],'k'); set(h,'LineWidth',2);
xlabel('Andrew Storm-relative u WE section');
print andrew_u_WE_section
%
figure(4); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(u(:,:,lp+1)),ucon); hold on; 
[cs,h] = contour(xx,z,squeeze(u(:,:,lp+1)),[0 0],'k'); set(h,'LineWidth',2);
xlabel('Andrew Storm-relative u SN section');
print andrew_u_SN_section
%
figure(5); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(v(:,lp+1,:)),vcon); hold on;
[cs,h] = contour(xx,z,squeeze(v(:,lp+1,:)),[0 0],'k'); set(h,'LineWidth',2);
xlabel('Storm-relative v WE section');
print andrew_v_WE_section
%
figure(6); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(v(:,:,lp+1)),vcon); hold on;
[cs,h] = contour(xx,z,squeeze(v(:,:,lp+1)),[0 0],'k'); set(h,'LineWidth',2);
xlabel('Storm-relative v SN section');
print andrew_v_SN_section
%
norm=sqrt( u.^2 + v.^2 ) ;
%
figure(7); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(norm(:,:,lp+1)),vcon); hold on; xlim([-250,250]); ylim([10,500]);
%[cs,h] = contour(xx,z,squeeze(norm(:,:,lp+1)),[0 0],'k'); set(h,'LineWidth',2);
[cs,h] = contour(xx,z,squeeze(norm(:,:,lp+1))); set(h,'LineWidth',2);
[cs,h] = contour(xx,z,squeeze(norm(:,:,lp+1)));clabel(cs);
xlabel('Wind Speed for Andrew at landfall: North-South section (m/s)');
print andrew_profile_SN_section
%
%
figure(8); clf;
subplot(1,1,1); cla; 
contourf(xx,z,squeeze(norm(:,lp+1,:)),vcon); hold on;xlim([-250,250]); ylim([10,500]);
%[cs,h] = contour(xx,z,squeeze(norm(:,lp+1,:)),[0 0],'k'); set(h,'LineWidth',2);
[cs,h] = contour(xx,z,squeeze(norm(:,lp+1,:))); set(h,'LineWidth',2);
[cs,h] = contour(xx,z,squeeze(norm(:,lp+1,:)));clabel(cs);
xlabel('Wind Speed for Andrew at landfall: East-West section (m/s)');
print andrew_profile_WE_section

fid = fopen('storm_prfl.txt','w');
fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',norm);
fclose(fid)
%norm(1,3)
%norm(1,11)
%
%
%

%[cs,h] = contourf(xkm,ykm, swf2); set(h,'LineStyle','none');
%[cs,h] = contour(xkm,ykm, swf2); clabel(cs,h);
%xlabel('Storm-relative asymetric w WE section');
%subplot(2,2,4); cla; 
%contourf(xx,z,squeeze(wa(:,:,lp+1)),wcon,'k'); hold on;
%[cs,h] = contour(xx,z,squeeze(wa(:,:,lp+1)),[0 0],'k'); set(h,'LineWidth',2);
%xlabel('Storm-relative asymetric w SN section');

[vmax,imax] = max(v,[],1); vmax = shiftdim(vmax); imax = shiftdim(imax);
zmax = z(imax);
fac1 = vmax./V;
fac2 = ( vmax - Ut*sin(lam) ) ./ ( V - Ut*sin(lam) );

%figure(5); clf;
%subplot(2,2,1); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%contourf(xx,yy,fac1); plot(cx,cy,'w');
%title(['V_m=' num2str(vm) ' r_m=' num2str(rm*1e-3) ' b=' num2str(b) ...
%      ' lat=' num2str(lat) ' K=' num2str(K) ' C=' num2str(C) ...
%      ' U_t=' num2str(Ut) ]);
%xlabel('Storm-relative jet factor');
%subplot(2,2,2); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%contourf(xx,yy,zmax); plot(cx,cy,'w');
%xlabel('Jet height');
%subplot(2,2,3); cla; hold on; set(gca,'DataAspectRatio',[1 1 1]); xlim([-150,150]); ylim([-150,150]);
%contourf(xx,yy,fac2); plot(cx,cy,'w');
%xlabel('Earth-relative jet factor');

%
