function [ ] = plotSpiralSounding( filename, date, startTime, endTime )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Developed by Matthew Sienkiewicz
%  Stony Brook University
%  10 July 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example: To create soundings of Theta-V, Mixing Ratio and Wind Speed from
% aircraft spirals, enter the following command on the command line:
%
% plotSpiralSounding('20130623.txt','23-Jun-2013',startTime,endTime);
%
% where startTime and endTime are datestrings:
%
% startTime = '23-Jun-2013 22:03:00';
%
% endTime = '23-Jun-2013 22:09:30';
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

flight = loadFlight(filename,date);
[~,ind1,~] = intersect(datenum(datestr(flight.time,...
        'dd-mmm-yyyy HH:MM:SS'),'dd-mmm-yyyy HH:MM:SS'),...
        datenum(startTime,'dd-mmm-yyyy HH:MM:SS'));
[~,ind2,~] = intersect(datenum(datestr(flight.time,...
        'dd-mmm-yyyy HH:MM:SS'),'dd-mmm-yyyy HH:MM:SS'),...
        datenum(endTime,'dd-mmm-yyyy HH:MM:SS'));
acft(:,1) = round(flight.alt(ind1:ind2,1)./10).*10;
acft(:,2) = flight.wsp(ind1:ind2);
acft(:,3) = flight.thetav(ind1:ind2);
acft(:,4) = flight.mix(ind1:ind2) .* 1000;
blah = nan(length([0:10:2000]),3);
blah(:,1) = (0:10:2000)';
for i=1:1:length(blah(:,1))
    ind = find(acft(:,1)==blah(i,1));
    if isempty(ind)
        blah(i,2) = NaN;
        blah(i,3) = NaN;
        blah(i,4) = NaN;
    else
        blah(i,2) = mean(acft(ind,2));
        blah(i,3) = mean(acft(ind,3));
        blah(i,4) = mean(acft(ind,4));
    end
end
acft = blah;

figure;
subplot(1,3,1);
plot(acft(:,3),acft(:,1),'Color','k','LineWidth',3);
set(gca,'XLim',[288 300],'YLim',[0 1200],'XGrid','on','YGrid','on','FontSize',14,'FontName','Times New Roman');
xlabel('Virtual Potential Temperature [K]','FontSize',20,'FontName','Times New Roman','FontWeight','b');
ylabel('Altitude [meters]','FontSize',20,'FontName','Times New Roman','FontWeight','b');

subplot(1,3,2);
plot(acft(:,4),acft(:,1),'Color','k','LineWidth',3);
set(gca,'XLim',[0 12],'YLim',[0 1200],'XGrid','on','YGrid','on','FontSize',14,'FontName','Times New Roman');
xlabel('Mixing Ratio [g kg^{-1}]','FontSize',20,'FontName','Times New Roman','FontWeight','b');
ylabel('Altitude [meters]','FontSize',20,'FontName','Times New Roman','FontWeight','b');

subplot(1,3,3);
plot(acft(:,2),acft(:,1),'Color','k','LineWidth',3);
set(gca,'XLim',[0 16],'XTick',[0:2:16],'YLim',[0 1200],'XGrid','on','YGrid','on','FontSize',14,'FontName','Times New Roman');
xlabel('Wind Speed [m s^{-1}]','FontSize',20,'FontName','Times New Roman','FontWeight','b');
ylabel('Altitude [meters]','FontSize',20,'FontName','Times New Roman','FontWeight','b');

end

function [ flight ] = loadFlight( filename, date )

rawFlight = load_flight( filename, date);

%Calculate WDIR
wdr = (180 / pi).*atan2((-1.*rawFlight(:,6)),(-1.*rawFlight(:,5)));
wdr(wdr<0) = 360 + wdr(wdr<0);

%Calculate WSP
wsp = sqrt((rawFlight(:,6).^2)+(rawFlight(:,5).^2));

clear flight
flight.time = rawFlight(:,1);
flight.lat = rawFlight(:,7);
flight.lon = rawFlight(:,8);
flight.alt = rawFlight(:,9)+37.7;
flight.wsp = wsp;
flight.wdr = wdr;
flight.tmp = rawFlight(:,2);
flight.hum = rawFlight(:,3);
flight.prs = rawFlight(:,4);
flight.u = rawFlight(:,6);
flight.v = rawFlight(:,5);
flight.w = rawFlight(:,17);
flight.vvel = rawFlight(:,12);
flight.theta = (flight.tmp + 273.15) .* ...
( (100000 ./ flight.prs) .^ ( 287 / 1004) );
e_s = 6.11 .* exp(  ( ( 2.5e6 .* 18.016 ) ./ (1000 .* 8.3145 ) ) ...
.* ( ( 1 ./ 273.15 ) - ( 1 ./ ( flight.tmp + 273.15 ) ) ) );
w_s = 0.622 .* ( e_s ./ ( ( flight.prs ./ 100 ) - e_s ) );
flight.mix = flight.hum .* w_s;
flight.thetav = flight.theta .* ( 1 + ( 0.61 .* flight.mix ) );
flight.head = rawFlight(:,15);
flight.tas = rawFlight(:,16);
flight.slip = rawFlight(:,18);
flight.roll = rawFlight(:,14);
flight.date = ...
    datenum(datestr(flight.time(1),'dd-mmm-yyyy'),'dd-mmm-yyyy');

end

function [ rawFlight ] = load_flight( filename, date )

flightday = datevec(date);

raw_flight = importdata(filename);

times = raw_flight(:,1);
startTime = times(1);
times = times - startTime;
times(times<0) = times(times<0) + 24.0;
times = times + startTime;
hours = floor(times);
times = (times-hours) .* 60;
minutes = floor(times);
seconds = (times - minutes) .* 60;

times = [zeros(length(times),1), zeros(length(times),1), ...
    zeros(length(times),1), hours, minutes, seconds];
times = datenum(times + repmat(flightday,length(times),1));

rawFlight = cat(2,times,raw_flight(:,2:end));

end