% 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;
interp_ok = 0;
folder = dir('20*');

for i_folder = 1:length(folder)
    if (folder(i_folder).isdir)
        cd(folder(i_folder).name);
        folder(i_folder).name
        files = dir('science_*.mat');
        
        for i = 1:length(files)
            disp(['file ' num2str(i)]);
            % Get the variables from the file we care about...
            try load(files(i).name, 'depth'); catch; end;
            try load(files(i).name, 'volume_scattering_470_nm');catch; end;
            try load(files(i).name, 'volume_scattering_650_nm');catch; end;
            try load(files(i).name, 'mass_concentration_of_chlorophyll_in_sea_water');catch; end;
            try load(files(i).name, 'sea_water_temperature');catch; end;
            try load(files(i).name, 'sea_water_salinity');catch; end;
            try load(files(i).name, 'sea_water_density');catch; end;
            try load(files(i).name, 'mole_concentration_of_nitrate_in_sea_water');catch; end;
            try load(files(i).name, 'mole_concentration_of_oxygen_in_sea_water');catch; end;
            try load(files(i).name, 'downwelling_photosynthetic_photon_flux_in_sea_water');catch; end;
            try load(files(i).name, 'latitude');catch; end;
            try load(files(i).name, 'longitude');catch; end;
            try
                load(files(i).name, 'DVL_micro');
                alt_value = DVL_micro.height_above_sea_floor.value;
                alt_time = DVL_micro.height_above_sea_floor.time;
            catch
                try
                    load(files(i).name, 'height_above_sea_floor');
                    alt_value = height_above_sea_floor.value;
                    alt_time = height_above_sea_floor.time;
                catch end;
            end;
            try load(files(1).name, 'platform_pitch_angle');catch; end;
            try load(files(1).name, 'platform_orientation');catch; end;
            
            % 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('depth')
                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;
            end;
            
            if exist('volume_scattering_470_nm')
                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;
            end;
            
            if exist('sea_water_temperature')
                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;
            end;
            
            if exist('mole_concentration_of_nitrate_in_sea_water')
                if ~exist('nit')
                    if length(mole_concentration_of_nitrate_in_sea_water.value) > 1
                        nit = mole_concentration_of_nitrate_in_sea_water.value;
                        tnit = mole_concentration_of_nitrate_in_sea_water.time;
                    end;
                else
                    if length(mole_concentration_of_nitrate_in_sea_water.value) > 1
                        nit = [nit mole_concentration_of_nitrate_in_sea_water.value];
                        tnit = [tnit mole_concentration_of_nitrate_in_sea_water.time];
                    end;
                end;
            end;
            
            if exist('mole_concentration_of_oxygen_in_sea_water')
                if ~exist('o2')
                    if length(mole_concentration_of_oxygen_in_sea_water.value) > 1
                        o2 = mole_concentration_of_oxygen_in_sea_water.value;
                        to2 = mole_concentration_of_oxygen_in_sea_water.time;
                    end;
                else
                    if length(mole_concentration_of_oxygen_in_sea_water.value) > 1
                        o2 = [o2 mole_concentration_of_oxygen_in_sea_water.value];
                        to2 = [to2 mole_concentration_of_oxygen_in_sea_water.time];
                    end;
                end;
            end;
            
            if exist('downwelling_photosynthetic_photon_flux_in_sea_water')
                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;
            end;
            
            if exist('latitude')
                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;
            end;
            
            if exist('alt_value')
                if ~exist('alt')
                    if length(alt_value) > 1
                        alt = alt_value;
                        talt = alt_time;
                    end
                else
                    if length(alt_value) > 1
                        alt = [alt alt_value];
                        talt = [talt alt_time];
                    end;
                end;
            end;
            
            if exist('platform_orientation')
                if ~exist('head')
                    if length(platform_orientation.value) > 1
                        head = platform_orientation.value;
                        thead = platform_orientation.time;
                    end
                else
                    if length(platform_orientation.value) > 1
                        head = [head platform_orientation.value];
                        thead = [thead platform_orientation.time];
                    end;
                end;
            end;
            
            if exist('platform_pitch_angle')
                if ~exist('pitch')
                    if length(platform_pitch_angle.value) > 1
                        pitch = platform_pitch_angle.value;
                        tpitch = platform_pitch_angle.time;
                    end
                else
                    if length(platform_pitch_angle.value) > 1
                        pitch = [pitch platform_pitch_angle.value];
                        tpitch = [tpitch platform_pitch_angle.time];
                    end;
                end;
            end;
            
            cd('..');
        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)

if interp_ok
    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);
end;

try; clear('depth'); catch end;
try; clear('volume_scattering_470_nm'); catch end;
try; clear('volume_scattering_650_nm'); catch end;
try; clear('mass_concentration_of_chlorophyll_in_sea_water'); catch end;
try; clear('sea_water_temperature'); catch end;
try; clear('sea_water_salinity'); catch end;
try; clear('sea_water_density'); catch end;
try; clear('mole_concentration_of_nitrate_in_sea_water'); catch end;
try; clear('mole_concentration_of_oxygen_in_sea_water'); catch end;
try; clear('downwelling_photosynthetic_photon_flux_in_sea_water'); catch end;
try; clear('latitude'); catch end;
try; clear('longitude'); catch end;
try; clear('DVL_micro'); catch end;
try; clear('height_above_sea_floor'); catch end;
try; clear('platform_pitch_angle'); catch end;
try; clear('platform_orientation'); catch end;