% InterpSondes_P.m - reads sondes for each flight and bin-averages to fixed
% pressure grid with 20mb bin width and 41 bins from 1000 to 200 mb (bin centers).

% matlab script path
path_prog = '/Users/byronb/Documents/MATLAB/2016_ENRR';
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
rehash;

path_data = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/data2/';
cd(path_data);

p_bin = 190:20:1010; p_bin = p_bin'; % 42 level fixed pressure bin edges 190-1010 mb

for ii=1:22 % iterate through 22 flights
    % read raw profile data
    fname = ['*_RF',sprintf('%02.0f',ii),'.mat'];
    file = dir(fname);
    disp(['Processing ',file.name]);
    snd_raw = load(fullfile(path_data,file.name), 'jdxh', 'zhh','Uhh','rhahh','tahh','lathh','lonhh','Phh','qahh','udirhh');
    N = length(snd_raw.jdxh); % number of drops on this flight
    z_raw = snd_raw.zhh;      % altitude vectors, bottom up
    P_filt = sgolayfilt(snd_raw.Phh,1,25);
    
    % wind speed
    U_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.Uhh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        U_bin(:,kk) = flipud(temp(:,2));
        if kk==1
            p_ref = flipud(temp(:,1)) + 10;  % reference pressure bin centers
        end
    end
    
    % wind direction 'from'
    dir_raw = snd_raw.udirhh;
    dir_raw = dir_raw*pi/180;
    dir_raw = unwrap(dir_raw);
    dir_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(dir_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        dir_bin(:,kk) = flipud(temp(:,2));
        dir_bin(:,kk) = dir_bin(:,kk)*180/pi;
        dir_bin(:,kk) = mod(dir_bin(:,kk)+8*360, 360);
    end
    
    % compute Un/Ue components, where Un is positive S to N and Ue is
    % positive for E to W (right handed convention).
    Un = U_bin .* cos(dir_bin*pi/180);
    Ue = U_bin .* sin(dir_bin*pi/180);
    
    % air temp
    ta_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.tahh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        ta_bin(:,kk) = flipud(temp(:,2));
    end
    
    % rh
    rh_raw = snd_raw.rhahh;
    rh_raw(isnan(rh_raw) & z_raw<12000) = 1;
    rh_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(rh_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        rh_bin(:,kk) = flipud(temp(:,2));
    end
    
    % qair
    qa_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.qahh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        qa_bin(:,kk) = flipud(temp(:,2));
    end
    
    % z
    z_bin = NaN(41,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(z_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        z_bin(:,kk) = flipud(temp(:,2));
    end
    
    % save interpolated profiles to new structure
    sname = [file.name(1:28),'_pinterp.mat'];
    snd_pinterp = struct('jd',snd_raw.jdxh,'lon',snd_raw.lonhh,'lat',snd_raw.lathh,...
        'P',p_ref,'z',z_bin,'U',U_bin,'udir',dir_bin,'Un',Un,'Ue',Ue,'ta',ta_bin,...
        'rh',rh_bin,'qa',qa_bin);
    save(sname,'snd_pinterp');
end



