%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%%    Plot the Valentia radiosonde data
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      1      2      3      4      5      6      7     8      9     10      11
% -----------------------------------------------------------------------------
%    PRES   HGHT   TEMP   DWPT   RELH   MIXR   DRCT   SKNT   THTA   THTE   THTV
%     hPa     m      C      C      %    g/kg    deg   knot     K      K      K 
% -----------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clf

%%%   Specify the input file
datafile = '03953-2011111000.dat';

%%%   Load in the data into an array called Table.
Table = load (datafile);

%%%   Select the pressure, height, temperature, etc.
Pressure      = Table(:,1);
Height        = Table(:,2);
Temperature   = Table(:,3);
Dewpoint      = Table(:,4);
Windspeed     = Table(:,8); 
Theta         = Table(:,9); 

%%%   Plot Temperature against Height.

% Open figure(1) for plotting the rotated axes.
figure(1)

plot(Temperature(1:25),Height(1:25),'b')
hold on
plot(Dewpoint(1:25),Height(1:25),'r')
title('03953 ascent')
xlabel('T,T_d')
ylabel('z(m)')
grid on
hold off
pause

%%%   Plot Windspeed and Dewpoint.
plot(Windspeed,Height,'m'); hold on
plot(Temperature,Height,'b')
plot(Dewpoint,Height,'r')
grid on
legend('V(kts)','T(^oC)','T_d(^oC)')
title('03953 ascent')
xlabel('V,T,T_d')
ylabel('z(m)')
hold off
pause

%%%   Plot Potential Temperature.
Tzero = 273.15;
plot(Theta-Tzero,Height,'k')
grid on
legend('\theta(^oC)')
title('03953 ascent')
xlabel('\theta')
ylabel('z(m)')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   Plot log(theta) versus Temperature.

% [left, bottom, width, height]
scrsz = get(0,'ScreenSize');
xmid = scrsz(3)/2;
width = scrsz(3)/2;
height = (40/25)*width;
figure('Position',[xmid 1 width height])

TKelvin = Temperature + Tzero;
h = semilogy(TKelvin,Theta,'b');
axis([200 300 270 480]);

grid on
title('03953 ascent')
xlabel('T')
ylabel('log(\theta)')

hold on
Rgas = 287;
cp = 1004;
kappa = Rgas/cp;
T = 200:300;
for P=10:10:100
  THETA = (P/1000)^(-kappa)*T;
  H = semilogy(T,THETA,'m');
end
for P=100:100:1000
  THETA = (P/1000)^(-kappa)*T;
  H = semilogy(T,THETA,'g');
end
%  Check: Plot 500 and 1000 hPa line in black.
for P=500:500:1000
  THETA = (P/1000)^(-kappa)*T;
  H = semilogy(T,THETA,'k');
end
% axis([200 300 250 1000])
pause


AX = axis;
ASP = height/width;
angle = atan2((AX(2)-AX(1)),ASP*(AX(4)-AX(3)))*(180/pi);
view(angle,90)
pause

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   Plot log(theta) versus Temperature on rotated axes.

% Open a new figure for plotting the rotated axes.
figure(3)

%  Rotate the axes.
XX0 = -30 + Tzero;
YY0 = +30 + Tzero;
XX =  TKelvin - XX0;
YY =  Theta - YY0;

rrt2 = 1/sqrt(2);
xx = XX0 + rrt2*( XX + YY);
yy = YY0 + rrt2*(-XX + YY);

hrot = plot(xx,yy,'b');
grid on
hold on
% title('03953 ascent')
% xlabel('T')
% ylabel('log(\theta)')
pause

%  Plot the isotherms
for TT=100:100:600
  TH = 100:10:600;
  XT = XX0 + rrt2*(( TT-XX0) + (TH-YY0));
  YT = YY0 + rrt2*((-TT-XX0) + (TH-YY0));
  htt = semilogy(XT,YT,'g');
  % axis('square')
end
%  Plot the isentropes
for TH=100:100:600
  TT = 100:10:600;
  XT = XX0 + rrt2*(( TT-XX0) + (TH-YY0));
  YT = YY0 + rrt2*((-TT-XX0) + (TH-YY0));
  htt = semilogy(XT,YT,'c');
  % axis('square')
  % axis('off')
end
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rotate the plot 45 degrees clockwise.
% zaxis = [0 0 1];
% rotate(htt,zaxis,-45)
view(45,90)

pause
%return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rotate the plot 45 degrees clockwise.
% zaxis = [0 0 1];
% rotate(htt,zaxis,-45)
AX = axis;
angle = atan2((AX(2)-AX(1)),(AX(4)-AX(3)))*(180/pi);
view(angle,90)
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

