function [ ] = printNCLCrossSection( filename, date, freq, times,...
    coords, outputFile )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Developed by Matthew Sienkiewicz
%  Stony Brook University
%  11 June 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example: To project a flight segment onto an X-Z plane and print data
%          to an ASCII file in order to be imported into an NCL script
%          that will plot the data:
%
% printNCLCrossSection('20130621.txt','21-Jun-2013',20,times,...
%                    coords,outputFile);
%
% where times is a string array containing the strings of times used to
% parse the flight file:
%
% times = ['21-Jun-2013 18:53:00';'21-Jun-2013 19:24:30'];
%
% and coords is an array containing the end points for the cross-section:
%
% coords = [41.092366,-70.435517;41.547146,-70.288170];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    point1.lat    = coords(1,1);
    point1.lon    = coords(1,2);
    point2.lat    = coords(2,1);
    point2.lon    = coords(2,2);
    
    startTime     = times(1,:);
    endTime       = times(2,:);

    flight = loadFlight(filename,date);
    flight = avgFlight(flight,freq);
    
    [~,indStart,~] = intersect(datenum(datestr(flight.time,...
        'dd-mmm-yyyy HH:MM:SS')),datenum(startTime,...
        'dd-mmm-yyyy HH:MM:SS'));
    [~,indEnd,~] = intersect(datenum(datestr(flight.time,...
        'dd-mmm-yyyy HH:MM:SS')),datenum(endTime,'dd-mmm-yyyy HH:MM:SS'));

    leg1 = [indStart:1:indEnd]';

    d = haversine(point1.lat,point1.lon,point2.lat,point2.lon);
    d = ceil(d/5000)*5;

    dist = getDistance(flight.lat(leg1),flight.lon(leg1),...
        point1.lat,point1.lon,point2.lat,point2.lon)./1000;

    theta = (flight.tmp + 273.15) .* ...
        ( ( 100000 ./ flight.prs ) .^ ( 287 / 1004 ) );

    avgInt = 2;

    for i=1:1:length(leg1)
        alt(i) =  mean(flight.alt(leg1(i)-avgInt:1:leg1(i)+avgInt));
        thet(i) =  round(mean(theta(leg1(i)-avgInt:1:leg1(i)+avgInt))...
            *10)/10;
        uWnd(i) =  round(mean(flight.u(leg1(i)-avgInt:1:leg1(i)+avgInt))...
            *10)/10;
        vWnd(i) =  round(mean(flight.v(leg1(i)-avgInt:1:leg1(i)+avgInt))...
            *10)/10;
    end

    p = figure;
    a = subplot(2,1,1);
    set(a,'XLim',[0 d],'YLim',[0 800],'FontName',...
        'Times New Roman','FontSize',14);
    set(a,'XTick',[0:5:d],'YTick',[0:100:800]);
    title('Potential Temperature [K]','FontName',...
        'Times New Roman','FontSize',20);
    xlabel('Distance [kilometers]','FontName','Times New Roman',...
        'FontSize',16,'FontWeight','b');
    ylabel('Height [meters]','FontName','Times New Roman',...
        'FontSize',16,'FontWeight','b');
    for m=1:1:length(thet)
        text(dist(m),alt(m),num2str(thet(m)),'FontName',...
            'Times New Roman','FontSize',6);
    end
    box on
    b = subplot(2,1,2);
    set(b,'XLim',[0 d],'YLim',[0 800],'FontName','Times New Roman',...
        'FontSize',14);
    set(b,'XTick',[0:5:d],'YTick',[0:100:800]);
    title('Wind Speed [m s^{-1}]','FontName','Times New Roman',...
        'FontSize',20);
    xlabel('Distance [kilometers]','FontName','Times New Roman',...
        'FontSize',16,'FontWeight','b');
    ylabel('Height [meters]','FontName','Times New Roman',...
        'FontSize',16,'FontWeight','b');

    fileID = fopen(horzcat(outputFile),'w');
    for m=1:1:length(uWnd)
        text(dist(m),alt(m),num2str(round(sqrt((uWnd(m)^2)+(vWnd(m)^2))...
            *10)/10),'FontName','Times New Roman','FontSize',6);
        formatSpec = '%8.3f\t%5d\t%5.1f\t%5.1f\t%6.1f\n';
        fprintf(fileID,formatSpec,round(dist(m)*1000)/1000,...
            round(alt(m)),uWnd(m),vWnd(m),thet(m));
    end
    box on
    fclose(fileID);

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 [ dist ] = getDistance( lats, lons, lat1, lon1, lat2, lon2)

    dist = nan(length(lats),1);
    a = haversine(lat1,lon1,lat2,lon2);
    for i=1:1:length(lats)
        b = haversine(lat1,lon1,lats(i),lons(i));
        c = haversine(lats(i),lons(i),lat2,lon2);
        A = acos((a*a+b*b-c*c)/(2*a*b));
        dist(i) = b * sin(degtorad(90) - A);
    end

end

function [ dist ] = haversine( lat1, lon1, lat2, lon2 )

    Lat1 = degtorad(lat1);
    Lat2 = degtorad(lat2);
    delLat = Lat2 - Lat1;
    Lon1 = degtorad(lon1);
    Lon2 = degtorad(lon2);
    delLon = Lon2 - Lon1;
    R = 6371000.;
    a = (sin(delLat/2)*sin(delLat/2))...
        + (cos(Lat1)*cos(Lat2)*sin(delLon/2)*sin(delLon/2));
    c = 2*atan2(sqrt(a),sqrt(1-a));
    dist = R*c;
    clear Lat1 Lat2 delLat Lon1 Lon2 delLon R a c

end
