function [ ] = printNCLFlightTrack( filename, date, freq, outputFile )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Developed by Matthew Sienkiewicz
%  Stony Brook University
%  11 June 2015
%
%  This MATLAB function will read in an ASCII data file
%  from the AIMMS-20 instrument and output another ASCII
%  file containing the 1-Hz latitude, longitude, u, and v
%  values. The output file is intended for input into the
%  NCL script plotFlightTrack.ncl.
%
%  Example:
%  printNCLFlightTrack('20140723.txt','23-Jul-2014',20,'20140723_ncl.txt');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    flight = loadFlight(filename,date);
    flight = avgFlight(flight,freq);

    fileID = fopen(outputFile,'w');
    for i=1:1:length(flight.lat)
        formatSpec = '%9.6f\t%10.6f\t%10.6f\t%10.6f\n';
        fprintf(fileID,formatSpec,flight.lat(i),flight.lon(i),...
            flight.u(i),flight.v(i));
    end
    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