% This script takes all the 'science*.mat' files in a folder, opens them,
% extracts key variables, and concatenates the variables into continuous
% vectors.   Ends up with a much smaller, more managable variable space.
%
% J. G. Bellingham
% Feb 13, 2014

clear;
fix_time = 0;
files = dir('science*.*');

for i = 1:length(files)
    i
    % Get the variables from the file we care about...
    load(files(i).name, 'depth');
    load(files(i).name, 'volume_scattering_470_nm');
    load(files(i).name, 'volume_scattering_650_nm');
    load(files(i).name, 'mass_concentration_of_chlorophyll_in_sea_water');
    load(files(i).name, 'sea_water_temperature');
    load(files(i).name, 'sea_water_salinity');
    load(files(i).name, 'sea_water_density');
    load(files(i).name, 'downwelling_photosynthetic_photon_flux_in_sea_water');
    load(files(i).name, 'latitude');
    load(files(i).name, 'longitude');
    load(files(i).name, 'DVL_micro');
    
    % Concatenate the variables for the entire section...if the variable
    % does not exist yet, then create it.  After that, just tack the new
    % values on.  Make sure there is actually data to be added on,
    % however...
    if ~exist('dep')
        if length(depth.value) > 1
            dep = depth.value;
            tdep = depth.time;
        end
    else
        if length(depth.value) > 1
            dep = [dep depth.value];
            tdep = [tdep depth.time];
        end;
    end;
    
    if ~exist('eco470')
        if length(volume_scattering_470_nm.value) > 1
            eco470 = volume_scattering_470_nm.value;
            teco = volume_scattering_470_nm.time;
            eco650 = interp1(volume_scattering_650_nm.time, volume_scattering_650_nm.value, volume_scattering_470_nm.time);
            ecochl = interp1(mass_concentration_of_chlorophyll_in_sea_water.time, mass_concentration_of_chlorophyll_in_sea_water.value, volume_scattering_470_nm.time);
        end;
    else
        if length(volume_scattering_470_nm.value) > 1
            eco470 = [eco470 volume_scattering_470_nm.value];
            teco = [teco volume_scattering_470_nm.time];
            eco650 = [eco650 volume_scattering_650_nm.value];
            ecochl = [ecochl mass_concentration_of_chlorophyll_in_sea_water.value];
        end;
    end;
    
    if ~exist('swtemp')
        if length(sea_water_temperature.value) > 1
            swtemp = sea_water_temperature.value;
            tsw = sea_water_temperature.time;
            swsalt = interp1(sea_water_salinity.time, sea_water_salinity.value, sea_water_temperature.time);
            swdens = interp1(sea_water_density.time, sea_water_density.value, sea_water_temperature.time);
        end;
    else
        if length(sea_water_temperature.value) > 1
            swtemp = [swtemp sea_water_temperature.value];
            tsw = [tsw sea_water_temperature.time];
            swsalt = [swsalt interp1(sea_water_salinity.time, sea_water_salinity.value, sea_water_temperature.time)];
            swdens = [swdens interp1(sea_water_density.time, sea_water_density.value, sea_water_temperature.time)];
        end;
    end;
    
    if ~exist('par')
        if length(downwelling_photosynthetic_photon_flux_in_sea_water.value) > 1
            par = downwelling_photosynthetic_photon_flux_in_sea_water.value;
            tpar = downwelling_photosynthetic_photon_flux_in_sea_water.time;
        end;
    else
        if length(downwelling_photosynthetic_photon_flux_in_sea_water.value) > 1
            par = [par downwelling_photosynthetic_photon_flux_in_sea_water.value];
            tpar = [tpar downwelling_photosynthetic_photon_flux_in_sea_water.time];
        end;
    end;
    
    if ~exist('lat')
        if length(latitude.value) > 1
            lat = latitude.value;
            tpos = latitude.time;
            long = interp1(longitude.time, longitude.value, latitude.time);
        end;
    else
        if length(latitude.value) > 1
            lat = [lat latitude.value];
            tpos = [tpos latitude.time];
            long = [long interp1(longitude.time, longitude.value, latitude.time)];
        end;
    end;
    
        if ~exist('alt')
        if length(depth.value) > 1
            dep = depth.value;
            tdep = depth.time;
        end
    else
        if length(depth.value) > 1
            dep = [dep depth.value];
            tdep = [tdep depth.time];
        end;
    end;
end;

% If there are any negative time values in the time variables, they have to
% be fixed, or it totally messes up the process of getting all the
% variables in the same time base.

if fix_time
    dt = diff(tdep);
    ref = tdep(length(tdep));
    dt(find(dt <= 0)) = 1/24/3600/100; % Make negative values 0.01 seconds
    new = [0 cumsum(dt)];
    tdep = new - new(length(new)) + ref;
    
    
    dt = diff(teco);
    ref = teco(length(teco));
    dt(find(dt <= 0)) = 1/24/3600/100; % Make negative values 0.01 seconds
    new = [0 cumsum(dt)];
    teco = new - new(length(new)) + ref;
    
    
    dt = diff(tsw);
    ref = tsw(length(tsw));
    dt(find(dt <= 0)) = 1/24/3600/100; % Make negative values 0.01 seconds
    new = [0 cumsum(dt)];
    tsw = new - new(length(new)) + ref;
    
    
    dt = diff(tpar);
    ref = tpar(length(tpar));
    dt(find(dt <= 0)) = 1/24/3600/100; % Make negative values 0.01 seconds
    new = [0 cumsum(dt)];
    tpar = new - new(length(new)) + ref;
    
    
    dt = diff(tpos);
    ref = tpos(length(tpos));
    dt(find(dt < 0)) = 1/24/3600/100; % Make negative values 0.01 seconds
    new = [0 cumsum(dt)];
    tpos = new - new(length(new)) + ref;
end;

% Now interpolate everything to a common time base (after selecting the
% appropriate data section)

tc = teco;  % the common time base...

cdep = interp1(tdep, dep, tc);

ceco470 = interp1(teco, eco470, tc);
ceco650 = interp1(teco, eco650, tc);
cecochl = interp1(teco, ecochl, tc);

cswtemp = interp1(tsw, swtemp, tc);
cswsalt = interp1(tsw, swsalt, tc);
cswdens = interp1(tsw, swdens, tc);

cpar = interp1(tpar, par, tc);

clat = interp1(tpos, lat, tc);
clong = interp1(tpos, long, tc);