function [ ] = flightTimeSeriesByLeg( filename, date, freq, chop )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Developed by Matthew Sienkiewicz
%  Stony Brook University
%  11 June 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example: To create flight leg plots for the 23-Jun-2013 flight,
% enter the following command on the command line:
%
% flightTimeSeriesByLeg('20130623.txt','23-Jun-2013',20,chop);
%
% where chop is an array containing the strings of times used to parse the
% flight file:
%
% chop = ['23-Jun-2013 20:36:15';'23-Jun-2013 20:55:00';...
%     '23-Jun-2013 20:59:30';'23-Jun-2013 21:21:00';...
%     '23-Jun-2013 21:31:30';'23-Jun-2013 21:40:15';...
%     '23-Jun-2013 21:52:00';'23-Jun-2013 22:00:00';...
%     '23-Jun-2013 22:03:00';'23-Jun-2013 22:09:15';...
%     '23-Jun-2013 22:18:30';'23-Jun-2013 22:23:30';...
%     '23-Jun-2013 22:28:00';'23-Jun-2013 22:31:00';...
%     '23-Jun-2013 23:07:00';'23-Jun-2013 23:16:00';...
%     '23-Jun-2013 23:22:30';'23-Jun-2013 23:24:00';...
%     '23-Jun-2013 23:24:45'];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Import file into structured array
flight = loadFlight(filename,date);

%Average flight data down to 1 Hz
flight = avgFlight(flight,freq);

%Get Indexes of times to parse file
indexes = getIndexes(flight,chop);
disp(indexes)

% Find nice upper/lower limits for wsp/tmp/rh/alt plots
minWSP = floor(min(flight.wsp)/5)*5;
disp(horzcat('minWSP = ',num2str(minWSP)))
maxWSP = ceil(max(flight.wsp)/5)*5;
disp(horzcat('maxWSP = ',num2str(maxWSP)))

minTMP = floor(min(flight.tmp)/2)*2;
disp(horzcat('minTMP = ',num2str(minTMP)))
maxTMP = ceil(max(flight.tmp)/2)*2;
disp(horzcat('maxTMP = ',num2str(maxTMP)))

minRH = floor(min(flight.hum)/0.1)*0.1;
disp(horzcat('minRH = ',num2str(minRH)))
maxRH = ceil(max(flight.hum)/0.1)*0.1;
disp(horzcat('maxRH = ',num2str(maxRH)))

minAlt = 0;
maxAlt = ceil(max(flight.alt)/500)*500;

% For each flight leg, plot multi-panel figure
[m,~] = size(indexes);
for i=1:1:m-1
    sTime = flight.time(indexes(i,1));
    eTime = flight.time(indexes(i+1,1));
    leg = num2str(i);
    if length(leg) < 2
        leg = horzcat('0',leg);
    end

    figure;
    set(gcf,'Units','inches');
    set(gcf,'Position',[ 0    0.4167   16.6667    8.1563]);
    set(gcf,'PaperUnits','inches');
    set(gcf,'Position',[ 0    0.4167   16.6667    8.1563],...
        'PaperPositionMode','auto');

    subplot(3,2,1);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.alt(indexes(i,1):indexes(i+1,1)),'LineWidth',3,'Color','g');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid','on',...
        'YGrid','on');
    ylabel(gca,'Altitute [m]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[minAlt maxAlt],'YTick',[minAlt:500:maxAlt]);
    datetick(gca,'x','keeplimits');
    
    subplot(3,2,2);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.tmp(indexes(i,1):indexes(i+1,1)),'LineWidth',3,'Color','k');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid','on',...
        'YGrid','on');
    ylabel(gca,'Temperature [^{\circ}C]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[minTMP maxTMP],'YTick',[minTMP:2:maxTMP]);
    datetick(gca,'x','keeplimits');
    
    subplot(3,2,3);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.head(indexes(i,1):indexes(i+1,1)),...
        'LineWidth',3,'Color','m');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid','on',...
        'YGrid','on');
    ylabel(gca,'Heading [deg]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[0 360],'YTick',[0:90:360]);
    datetick(gca,'x','keeplimits');
    
    subplot(3,2,4);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.hum(indexes(i,1):indexes(i+1,1)),'LineWidth',3,'Color','k');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid',...
        'on','YGrid','on');
    ylabel(gca,'Relative Humidity','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[minRH maxRH],'YTick',[minRH:0.1:maxRH]);
    datetick(gca,'x','keeplimits');
    
    subplot(3,2,5);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.wsp(indexes(i,1):indexes(i+1,1)),'LineWidth',3,'Color','k');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid','on',...
        'YGrid','on');
    xlabel(gca,'Time [UTC]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    ylabel(gca,'Wind Speed [m s^{-1}]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[minWSP maxWSP],'YTick',[minWSP:5:maxWSP]);
    datetick(gca,'x','keeplimits');
    
    subplot(3,2,6);
    plot(flight.time(indexes(i,1):indexes(i+1,1)),...
        flight.wdr(indexes(i,1):indexes(i+1,1)),'LineWidth',3,'Color','k');
    set(gca,'FontName','Times New Roman','FontSize',12,'XGrid','on',...
        'YGrid','on');
    xlabel(gca,'Time [UTC]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    ylabel(gca,'Wind Direction [deg]','FontName','Times New Roman',...
        'FontSize',14,'FontWeight','b');
    title(horzcat('Flight Leg ',leg,': ',datestr(sTime,'HH:MM:SS'),...
        ' - ',datestr(eTime,'HH:MM:SS')),'FontName','Times New Roman',...
        'FontSize',12,'FontWeight','b');
    set(gca,'XLim',[sTime eTime]);
    set(gca,'YLim',[0 360],'YTick',[0:90:360]);
    datetick(gca,'x','keeplimits');

%     print(gcf,'-dtiff','-r300',horzcat('C:\Users\Matthew\Desktop\',...
%         datestr(flight.date,'yyyymmdd'),'_leg_',leg));
%     close all;
end

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


function [ output ] = avgFlight( flight, freq )

output.date = flight.date;

data = cat(2,flight.time,flight.lat,flight.lon,flight.alt,flight.wsp,...
    flight.wdr,flight.tmp,flight.hum,flight.prs,flight.u,flight.v,...
    flight.w,flight.head,flight.tas,flight.slip,flight.thetav,flight.mix);

hold_data = nan(length(data),17);
k=1;
for i=freq:freq:length(data)
    hold_data(k,:) = mean(data((i-(freq-1)):i,:));
    k = k + 1;
end

data = hold_data(1:(k-1),:);
output.time = data(:,1);
output.lat = data(:,2);
output.lon = data(:,3);
output.alt = data(:,4);
output.wsp = data(:,5);
output.wdr = data(:,6);
output.tmp = data(:,7);
output.hum = data(:,8);
output.prs = data(:,9);
output.u = data(:,10);
output.v = data(:,11);
output.w = data(:,12);
output.head = data(:,13);
output.tas = data(:,14);
output.slip = data(:,15);
output.thetav = data(:,16);
output.mix = data(:,17);

end


function [ indexes] = getIndexes( flight, chop )

[m,~] = size(chop);
indexes = nan(m,1);
for i=1:1:m
    [~,indexes(i,1),~] = intersect(datenum(datestr(flight.time,...
        'dd-mmm-yyyy HH:MM:SS'),'dd-mmm-yyyy HH:MM:SS'),...
        datenum(chop(i,:),'dd-mmm-yyyy HH:MM:SS'));
end

end