sname=input('Storm number (eg al092003)','s');
sdate='Post';
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\time.out');
time=load(sname1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\radius.out');
radius=load(sname1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\rgraph.out');
rgraph=load(sname1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\scatt.out');
scatt=load(sname1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\autoday.out');
autoday=load(sname1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\track.out');
track=load(sname1);
r=rgraph;
rb=radius(:,1);
rbi=radius(:,2);
rt=radius(:,3);
time1=time(:,1)./24;
at=size(time1);
aday=autoday(1,1);
sname1=strcat('c:\documents and settings\kerry emanuel\my documents\chip\storms\',sname,'\',sdate,'\title.out');
fid=fopen(sname1);
line=fgetl(fid);
tempp=size(line);
temps=tempp(1,2);
aend=temps-4;
abeg=aend-1;
bend=temps-7;
bbeg=bend-1;
cend=temps-1;
cbeg=temps-2;
ds=line(cbeg:cend);
hour0=str2num(ds);
as=line(abeg:aend);
bs=line(bbeg:bend);
month=str2num(bs);
day0=str2num(as);
time1=time1+hour0./24;
months={'January','February','March','April','May','June','July','August','September','October','November','December'};
v=time(:,2);
for ii=1:at,
   uu=time1(ii)-time1(1);
   if uu <= autoday
      jt=ii;
   end 
end 
tshort=time1(1:jt);
vshort=v(1:jt);
vobs=time(:,3);
timeb=scatt(:,1)./24+hour0./24;
va=scatt(:,2);
vb=scatt(:,3);
vc=scatt(:,4);
at=size(timeb);
for ii=1:at,
   uu=timeb(ii)-timeb(1);
   if uu <= autoday
      jb=ii;
   end 
end 
tbshort=timeb(1:jb);
vbshort=vb(1:jb);
%if vship(ags,1) == 0
h=plot(time1,vobs,time1,v);
set(gca,'fontweight','bold','fontsize',12)
%else   
%   h=plot(time1,v,time1,vship,time1,vobs);
%end
set(h,'linewidth',3)
hold on
g=plot(tshort,vshort,'r');
set(g,'linewidth',3)
hold off
axv=axis;
y1=1*floor((axv(1))/1.);
y2=1*ceil((axv(2))/1.);
%axis([y1,y2,-Inf,Inf]);
jg=get(gca,'xticklabel');
hj=str2num(char(jg));
del=hj(2)-hj(1);
delm=max(del,1);
yvec=[y1:delm:y2];
set(gca,'xtick',yvec);
yvect=transpose(yvec);
u=cellstr(int2str(yvect));
set(gca,'xticklabel',u)
xlabel('Time (days)','fontweight','bold')         
title(line,'fontweight','bold')
%if vship(ags,1) == 0
%   legend('CHIPS',line3,0)
%else   
%   legend('CHIPS',line2,line3,0)
%end   
hold on
gs=size(vc);
ags=gs(:,1);
%if vc(ags,1) == 0
   k=scatter(timeb,va,'filled');
%else   
%   k=scatter(timeb,vc,'g','filled');
%   hold on
%   j=scatter(timeb,va,'r','filled');
%end    
hold off
hold on
g=scatter(timeb,vb,'g','filled');
gg=scatter(tbshort,vbshort,'r','filled');
hold off
jg=get(gca,'xticklabel');
hj=str2num(char(jg))+day0;
asiz=size(hj);
asize=asiz(1,1);
iswitch=0;
jswitch=0;
for iss=1:asize,
   if month==4|month==6|month==9|month==11
      if hj(iss)>30
         hj(iss)=hj(iss)-30;
         iswitch=1;
      end   
   elseif month == 3
         if hj(iss)>28
            hj(iss)=hj(iss)-28;
            iswitch=1;
         end
   else
         if hj(iss)>31
            hj(iss)=hj(iss)-31;
            iswitch=1;
         end
      end
      if hj(iss)< 1
         jswitch=1;
         if month == 3
            hj(iss)=hj(iss)+28;
         elseif month == 5|month==7|month==10|month==12
            hj(iss)=hj(iss)+30;
         else
            hj(iss)=hj(iss)+31;
         end
     end    
end      
hjs=cellstr(int2str(hj));
set(gca,'xticklabel',hjs);
legend('Best Track','CHIPS',0)
%legend('Observed','Model',0)
if jswitch ==0
   if iswitch==0
      ixlab=months(month);
   else
      if month<12
         ixlab=strcat(months(month),' - ',months(month+1));
      else
         ixlab={'December - January'};
      end
   end
else
   if iswitch==0
      if month==1
         ixlab={'December - January'};
      else
         ixlab=strcat(months(month-1),' - ',months(month));
      end
   else
      if month==1
         ixlab={'December - February'};
      elseif month==12
         ixlab={'November - January'};
      else
         ixlab=strcat(months(month-1),' - ',months(month+1));
      end
   end   
end
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Maximum surface wind speed (knots)','fontweight','bold')
%
i=1;
while i == 1
it=input('Plot choices: \n \n Evolution: \n   1) V_max \n   3) Track \n   4) V_pot \n   5) V_shear \n   6) r_max  \n   7) p_min  \n   8) Updraft mass flux  \n   9) Downdraft mass flux \n  10) Storm center ocean mixed layer depth \n  11) Storm center SST anomaly \n  22) Temporary variable \n \n Radial plots:  \n   12) Surface azimuthal velocity \n   13) Tropopause azimuthal velocity \n   14) Radial velocities \n   15) Vertical mass fluxes \n   16) Surface pressure \n   17) Entropy \n   18) Mixed layer depth \n   19) Sea surface temperature anomaly \n \n Hovmuller diagrams: \n   20) Updraft mass flux \n   21) Azimuthal velocity \n \n --------------------- \n   50) Print present figure \n    0) Exit \n ');
if it == 1
v=time(:,2);
for ii=1:at,
   uu=time1(ii)-time1(1);
   if uu <= autoday
      jt=ii;
   end 
end 
tshort=time1(1:jt);
vshort=v(1:jt);
vobs=time(:,3);
timeb=scatt(:,1)./24+hour0./24;
va=scatt(:,2);
vb=scatt(:,3);
vc=scatt(:,4);
at=size(timeb);
for ii=1:at,
   uu=timeb(ii)-timeb(1);
   if uu <= autoday
      jb=ii;
   end 
end 
tbshort=timeb(1:jb);
vbshort=vb(1:jb);
%if vship(ags,1) == 0
h=plot(time1,vobs,time1,v);
set(gca,'fontweight','bold','fontsize',12)
%else   
%   h=plot(time1,v,time1,vship,time1,vobs);
%end
set(h,'linewidth',3)
hold on
g=plot(tshort,vshort,'r');
set(g,'linewidth',3)
hold off
set(gca,'xtick',yvec);
xlabel('Time (days)','fontweight','bold')         
title(line,'fontweight','bold')
%if vship(ags,1) == 0
%   legend('CHIPS',line3,0)
%else   
%   legend('CHIPS',line2,line3,0)
%end   
hold on
gs=size(vc);
ags=gs(:,1);
%if vc(ags,1) == 0
   k=scatter(timeb,va,'filled');
%else   
%   k=scatter(timeb,vc,'g','filled');
%   hold on
%   j=scatter(timeb,va,'r','filled');
%end    
hold off
hold on
g=scatter(timeb,vb,'g','filled');
gg=scatter(tbshort,vbshort,'r','filled');
hold off
jj=gca;
jg=get(jj,'xticklabel');
set(jj,'xticklabel',hjs);
legend('Best Track','CHIPS',0)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Maximum surface wind speed (knots)','fontweight','bold')
%
elseif it == 3
   legend off
   levels=[200:100:1000];
   dellat=5;
   x=track(:,6);
   y=track(:,5);
   xmin=min(x)-dellat;
   xmax=max(x)+dellat;
   ymin=min(y)-dellat;
   ymax=max(y)+dellat;
   if xmax>360
      xmax=xmax-360;
   end
   if xmin<0
      xmin=xmin+360;
   end   
   m_proj('gall-peters','long',[xmin xmax], ...
       'lat',[ymin ymax]);
   a=size(x);
   a1=a(1,1);
   b1=a1-5;
   x1=x(1:b1);
   y1=y(1:b1);

    x2=x(b1:a1);
   y2=y(b1:a1);
   m_coast('patch',[.6 .6 .6],'edgecolor','none');
   %  Use following for very high reolution coastlines:
   %
  % m_gshhs_h('patch',[.7 .7 .7],'edgecolor','none');
   hold on
   m_plot(x1,y1,'b.','MarkerSize',18)
   m_plot(x1,y1,'b','linewidth',2)
   m_plot(x2,y2,'b.','MarkerSize',18)
   m_plot(x2,y2,'b','linewidth',2)
   hold off
   m_grid
   xlabel('Longitude','fontweight','bold')
   ylabel('Latitude','fontweight','bold')
   title(line,'fontweight','bold')
   %m_elev('contour',levels,'edgecolor','k');
elseif it == 4
vpot=time(:,5);
h=plot(time1,vpot);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Potential intensity (kts)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 5
vshear=time(:,14);
h=plot(time1,vshear);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Shear (kts)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 6
rmax=time(:,4);
h=plot(time1,rmax);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Radius of maximum winds (km)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 7
pc=time(:,7);
ptheory=time(:,8);
h=plot(time1,ptheory);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Central surface pressure (mb)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 8
mu=time(:,9);
h=plot(time1,mu);
set(h,'linewidth',3);
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Convective updraft mass flux (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 9
md=time(:,10);
h=plot(time1,md);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Convective downdraft mass flux (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 10
hmc=time(:,11);
hmc0=time(:,13);
h=plot(time1,hmc,'-',time1,hmc0,'--');
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Storm center ocean mixed layer depth (m)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 11
sstc=time(:,12);
h=plot(time1,sstc);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
jj=gca;
set(jj,'xticklabel',hjs);
set(jj,'fontweight','bold','fontsize',12)
xlabel(ixlab,'fontweight','bold','fontsize',12)
ylabel('Storm center SST anomaly (K)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 12
v=radius(:,4);
h=plot(rb,v);
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Surface azimuthal velocity (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 22
temp=time(:,14);
h=plot(time1,temp);
set(h,'linewidth',3)
set(gca,'xtick',yvec);
xlabel('Time (hours)','fontweight','bold')         
ylabel('Temporary variable','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 13
vt=radius(:,5);
h=plot(rb,vt);
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Tropopause azimuthal velocity (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 16
p=radius(:,6);
h=plot(rbi,p);
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Surface pressure (mb)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 14
u0=radius(:,7);
u=radius(:,8);
h=plot(rb,u0,'-',rb,u,'--');
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Radial velocities (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 15
mu=radius(:,9);
we=radius(:,10);
h=plot(rbi,mu,'-',rbi,we,'--');
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Updraft and subsidence velocities (m/s)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 17
the=radius(:,11);
thes=radius(:,12);
them=radius(:,13);
h=plot(rbi,the,'-',rbi,thes,'--',rbi,them,'-.');
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('PBL, saturation and mid-trop theta_e','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 18
hm=radius(:,14);
h=plot(rbi,hm);
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Mixed layer depth (m)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 19
sst=radius(:,15);
h=plot(rbi,sst);
set(h,'linewidth',3)
xlabel('Radius (km)','fontweight','bold')
ylabel('Sea surface temperature anomaly (K)','fontweight','bold')
title(line,'fontweight','bold')
elseif it == 20
   load Graphics\whov.out
b=max(max(whov));
a=min(min(whov));
contourf(r,time1,whov,15)
title(['Cumulus mass flux (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
xlabel('Radius (km)','fontweight','bold')
ylabel ('time (hours)','fontweight','bold')
elseif it == 21
   load Graphics\vhov.out
b=max(max(vhov));
a=min(min(vhov));
contourf(r,time1,vhov,15)
title(['Azimuthal velocity (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
xlabel('Radius (km)','fontweight','bold')
ylabel ('time (hours)','fontweight','bold')
elseif it == 50
print -dpsc temp.ps
elseif it == 0
break
end
end
