% This script will make a Stuve diagram using data from a launch, plotting
% temperature, dewpoint, and altiude related to each other. This is a
% standard meteorological chart used for weather forecasting

% THIS SCRIPT ONLY LOOKS AT ASCENT
% Descent relative humidity readings are often inaccurate

% TO MAKE THIS SCRIPT WORK
% Your data must be formated in column vectors like this:
% dwptC || hghtM || presHPA || relhP || tempC 

% THE ONLY THING YOU NEED TO CHANGE:
% The title of the plot - the very last line of the code

% You will have to save the image on your own, based on where you want it
% to live and what you want to call it

% Written by: Tyler Boyle

figure

burst = find(presHPA == min(presHPA));

% To label the graph with altitude for each pressure level
% Because readings from these flights aren't often at the exact integer, we
% calculate out the altitude closest to the integer
min900 = abs(900 - presHPA);
min900index = find(min900 == min(min900));

min800 = abs(800 - presHPA);
min800index = find(min800 == min(min800));

min700 = abs(700 - presHPA);
min700index = find(min700 == min(min700));

min600 = abs(600 - presHPA);
min600index = find(min600 == min(min600));

min500 = abs(500 - presHPA);
min500index = find(min500 == min(min500));

min400 = abs(400 - presHPA);
min400index = find(min400 == min(min400));

min300 = abs(300 - presHPA);
min300index = find(min300 == min(min300));

min200 = abs(200 - presHPA);
min200index = find(min200 == min(min200));

min100 = abs(100 - presHPA);
min100index = find(min100 == min(min100));

% Create a vector that store the altitude values
axisAlt = [hghtM(min100index,1); hghtM(min200index,1); hghtM(min300index,1);...
    hghtM(min400index,1); hghtM(min500index,1); hghtM(min600index,1); ...
    hghtM(min700index,1); hghtM(min800index,1); hghtM(min900index,1); ...
    hghtM(5,1)];
axisPres = (100:100:1000);

hold on


% Here we are plotting this on 2 y axis. It's just easier to make sure the
% tick marks line up like this

% Left axis is "plotting" the temperature
yyaxis left
semilogy(tempC(1:burst,:),presHPA(1:burst,:),'linewidth',2)
ylim([100 1000])
set(gca,'ydir','reverse')
ylabel('Pressure (mb)')

% Right axis is "plotting" the dewpoint
yyaxis right
semilogy(dwptC(1:burst,:),presHPA(1:burst,:),'linewidth',2)
ylim([100 1000])
yticks(100:100:1000)
yticklabels(axisAlt(:,1))
set(gca,'ydir','reverse')
ylabel('Altitude (m)')

grid on
xticks(-120:10:40)
xlabel('Temperature (°C)')
legend('Temp (°C)','Dwpt (°C)')

set(gca,'Fontsize', 15)

% Change the title based on what you want your chart to read
title('NS62 Stuve Diagram','Fontweight','bold')