%Program:  stuve.m 
%Purpose:  To calculate a Stuve diagram for a given sounding.
%          Axes are p^(.2854) and T.
%Input:      ff = stuve(p[mbar],T[C],r[g/kg],Tdpt[C])
%                 p = press
%                 T = air temp
%                 r = mix ratio
%                 Tdpt = dew point
%Reference:  Bolton, D (1980), MWR,108:1046-1053 and Gill, 1980.
%          Note:  The correction to Cpa due to water vapor not made
%          when calculating contour lines of thetae.  This correction
%          is also not made in value for kappa which is used in theta
%          contour lines as well as thetae.
%
%Author: yolande serra
%Date:  February 2  1994
%------------------------------------------------------------

 function ff=stuve(p,T,r,Tdpt)

 if nargin<1
   help stuve
   return;
 end

 Lv=2.501e6;
 Cpa=1005.7;
 kappa=0.2854;
 str='';

%Temperature (x) and pressure (y) ranges for contours
 flag=input('Is this a shallow profile, <y> or <n>?  ', 's');
 if strcmp(flag,'n')
   x=[-20:10:30];
%   y=[420:60:1020];
   y=[360:60:1020];
 else
   x=17:34;
   y=840:10:1020;
 end
%Set pressure to pressure^kappa for Stuve diagram
 y = y.^kappa;
 [X,Y] = meshgrid(x,y);

%Get current axis handle and position or create one
 ax=get(gcf,'children');
 if isempty(ax)
   ax=gca;
   set(gca,'userdata','stuve')
   pos=get(gca,'pos');
 else
   if strcmp(get(gca,'userdata'),'stuve')
     pos=get(gca,'pos');
   else
     clf
     ax=gca;
     set(gca,'userdata','stuve')
     pos=get(ax,'pos');
   end
 end
 c=get(gca,'colororder');

%Contours of theta, thetae and rsat
 Z = (X+273.16).*((1000^kappa)./Y);
%If there is more than one plot on the page use fewer contour lines
 if length(ax)>1
   v=[240 260 280 300 320 340 360];
 else
   v=[240 250 260 270 280 290 300 310 320 330 340 350 360];
 end
 [c1,h1]=contour(x,y,Z,v,'-');
 set(h1,'color',c(1,:))
 hx1=clabel(c1);
 set(hx1,'color',c(1,:))
 set(gca,'YDir','reverse','xlim',[min(x) max(x)],'ylim',[min(y) max(y)])
 ytickv=get(gca,'YTick');
 ylimits=get(gca,'YLim');
 xlimits=get(gca,'XLim');
 ytickl=round(ytickv.^(1/kappa));
 lgth=length(ytickl);
 for i=1:lgth
   if ytickl(i) < 1000.0
     str=[str; [' ' num2str(ytickl(i))]];
   else
     str=[str; num2str(ytickl(i))];
   end
 end
 set(gca,'YTickLabel',str)
 hold on
%saturated equivalent potential temperature contours
 term1 = 17.67*X;
 term2 = X  + 243.5;
 Esat = 6.112.*exp(term1./term2);
 Rsat = 0.622*Esat./(Y.^(1/kappa) - Esat);
 Qsat = Rsat./(Rsat + 1);
 term4 = Lv.*Qsat./Cpa;
 ZZ = Z.*exp(term4./(X+273.16));
 if length(ax)>1
   v2=[255 283.4 335.5];
 else
   v2=[255 267.7 283.4 304.4 335.5 386.2]; 
 end
 [c2,h2]=contour(x,y,ZZ,v2,'-.');
 set(h2,'color',c(2,:))
 hx2=clabel(c2);
 set(hx2,'color',c(2,:))
 hold on
%saturated mixing ratio contours
 Rsat=Rsat.*1000.0;
 if length(ax)>1
   v3=[2 6 10 16 30];
 else
   v3=[1 2 4 6 8 10 12 16 20 30 40 50];
 end
 [c3,h3]=contour(x,y,Rsat,v3,':');
 set(h3,'color',c(3,:))
 hx3=clabel(c3);
 set(hx3,'color',c(3,:))

%Plot soundings of T and Tdpt onto contours
 r=r/1000.0;
 kappax=2/7*(1-0.23.*r);
 ydata=NaN*ones(size(p));
 ydata(isnan(p)==0)=(p(isnan(p)==0)).^kappax(isnan(p)==0);
 p1=plot(T,ydata,'m-','linewidth',0.5);
 hold on
 p2=plot(Tdpt,ydata,'c-','linewidth',1);
 hold off
 title('Stuve Diagram')
 ii=find(isnan(ydata)==0&(p>700&p<800));
 ytx=ydata(ii(1));
 xtx=T(ii(1));
 stext(xtx+0.4,ytx,'\14\times\demi\magenta T');
 ii=find(isnan(ydata)==0&(p>600&p<700));
 ytx=ydata(ii(2));
 xtx=Tdpt(ii(2));
 stext(xtx+0.4,ytx,'\14\times\demi\cyan T_d');
 xlabel('Temperature, ^o C'), ylabel('Pressure, mb')

%Make an axis with height values corresponding to pressure values
 xp=ytickv.^(1/kappa);
%..this polyfit came from snd1 on Jan 17 which gives ht(p) accurately
 H=polyval([0.0060 -20.4436 14557.0],xp);
 H=round(H/100)/10;
 lgth=length(H);
 strH=num2str(H(1));
 text(xlimits(2)+0.02*xlimits(2),ytickv(1),strH(1),'units','data')
 for i=2:lgth
   text(xlimits(2)+0.02*xlimits(2),ytickv(i),num2str(H(i)),'units','data')
 end
 text(1.08,.5,'(approx.) Altitude,  km','units','norm','horiz','center','rot',-90)

