clear;

% specify the filenames:
MPR_filename = 'st1_AE-MPR068a.dat';
HEADER_filename = 'st1_AE-HEADER.dat';
data_workspace = 'AE_hypro.mat';



% read in the MPR file, which contains the tilt/depth/timestamp data for each
% cast ...

  % ... open the file and read in the data:
  fid  = fopen(MPR_filename);
  data = textscan(fid, '%f%f%f%f%f%f%f%f%f%f%f%f%f%s%s%s%s%s','headerlines',1,'delimiter',',');

  % ... pull out the press/tilt/time fields:
  press = data{7};
  tiltx = data{8};
  tilty = data{9};
  press_times = char(data{16}); 
  mins=str2num(press_times(:,4:5)); 
  secs=str2num(press_times(:,7:end));
  
  
  % ... convert pressure times into single julian date:
  press_dates=nan(7253,6);
  press_dates(:,1)=2011;
  press_dates(:,2)=7;
  press_dates(:,3)=19;
  press_dates(:,4)=13;
  press_dates(:,5)=mins;
  press_dates(:,6)=secs;
  press_timestamp=datenum(press_dates);
  
  % ... close the file and clean up:
  fclose(fid);
  clear data

  
% read in the HEADER file, which contains the pressure tare value ...

  % ... open the file:
  fid  = fopen(HEADER_filename);
  data = textscan(fid, '%s%s%s%s','delimiter',',');

  % ... pull out the press/tilt/time fields:
  header_data = data{2};
  press_tare = str2num(header_data{19});
  
  % ... close the file:
  fclose(fid);


% read the profile data from workspaces produced by Haidi ...

  % ... load the workspace:
  load(data_workspace);
  
  % ... rename the variables:
  Lu174_times = p174_time;
  Lu174_data = Lu_174;
  
  Ed186_times = p186_time;
  Ed186_data = Ed_186;
  
  Ed187_times = p187_time;
  Ed187_data = Ed_187;
  
  
% interpolate data onto a common time/pressure scale:   

  % ... first, interpolate Ed186/Lu174 pressures from the reference
  % pressures:
  Ed186_press=interp1(press_timestamp, press, Ed186_times);
  Lu174_press=interp1(press_timestamp, press, Lu174_times);
 
  % ... compile the data:
  AE_Ed186_data = [Ed186_times Ed186_press Ed186_data]; 
  AE_Lu174_data = [Lu174_times Lu174_press Lu174_data];
  
  % ... interpolate tilt angle data to the Ed186/Lu174 time and 
  % pressure scales:
  Ed186_tiltx = interp1(press_timestamp, tiltx, Ed186_times);
  Ed186_tilty = interp1(press_timestamp, tilty, Ed186_times);

  Lu174_tiltx = interp1(press_timestamp, tiltx, Lu174_times);
  Lu174_tilty = interp1(press_timestamp, tilty, Lu174_times);
  
  
% kill all the profile values with tilt >5deg:
AE_Ed186_data(abs(Ed186_tiltx)>5 | abs(Ed186_tilty)>5,:)=[];
AE_Lu174_data(abs(Lu174_tiltx)>5 | abs(Lu174_tilty)>5,:)=[];
  
% correct the pressures (taking into account the pressure tare as well
% as the distance from the pressure sensor to each radiometer): 
AE_Ed186_data(:,2) = AE_Ed186_data(:,2) - press_tare - 1.1;
AE_Lu174_data(:,2) = AE_Lu174_data(:,2) - press_tare - 0.8;

% remove upcast data and above-surface data ...

  % ... identify the upcast measurements:
  Ed186_diffs = [0; diff(AE_Ed186_data(:,2))];
  Lu174_diffs = [0; diff(AE_Lu174_data(:,2))];
  
  % ... kill all measurements during upcasts: 
  AE_Ed186_data(Ed186_diffs<0,:)=[];
  AE_Lu174_data(Lu174_diffs<0,:)=[];
  
  % ... kill all measurements with negative depths:
  AE_Ed186_data(AE_Ed186_data(:,2)<0,:)=[];
  AE_Lu174_data(AE_Lu174_data(:,2)<0,:)=[];

  
% add cast numbers to the arrays ...

  % ... append a column of nans to the array:
  AE_Ed186_data = [nan(size(AE_Ed186_data,1),1) AE_Ed186_data];
  AE_Lu174_data = [nan(size(AE_Lu174_data,1),1) AE_Lu174_data];
  
  % ... identify the indices where each cast starts: 
  Ed186_casts_index = [1; find(diff(AE_Ed186_data(:,3))<0)+1];
  Lu174_casts_index = [1; find(diff(AE_Lu174_data(:,3))<0)+1];
  
  % ... write in the cast numbers:
  for n=1:length(Ed186_casts_index)
    
    if n<length(Ed186_casts_index)
      AE_Ed186_data(Ed186_casts_index(n):Ed186_casts_index(n+1),1)=n;
      AE_Lu174_data(Lu174_casts_index(n):Lu174_casts_index(n+1),1)=n;
    else
      AE_Ed186_data(Ed186_casts_index(n):end,1)=n;
      AE_Lu174_data(Lu174_casts_index(n):end,1)=n;
    end
    
  end
  
  
% bin the individual profiles from each cast ...
 
  % create bin vectors:
  bin1_step = 0.20; bins1 = 0.10:bin1_step:20; 
  bin2_step = 1.00; bins2 = 0.50:bin2_step:20;    

  % pre-allocate space for the binned cast arrays:
  AE_Ed186_bin1 = nan(length(bins1),138,12);
  AE_Ed186_bin2 = nan(length(bins2),138,12);
  
  AE_Lu174_bin1 = nan(length(bins1),138,12);
  AE_Lu174_bin2 = nan(length(bins2),138,12);
  
  
  % loop through each cast and bin vector and average the profile data:
  for n = 1:12 
    
    % pull out the data for this cast:
    Ed186_cast = AE_Ed186_data(AE_Ed186_data(:,1)==n,:);
    Lu174_cast = AE_Lu174_data(AE_Lu174_data(:,1)==n,:);
    
    % loop through bins1 and average the data:
    for m = 1:length(bins1)
    
        Ed186_bin1_index = find( ... 
          Ed186_cast(:,3)>=bins1(m)-bin1_step/2 & ... 
          Ed186_cast(:,3)<=bins1(m)+bin1_step/2);
        Ed186_bin1_data = Ed186_cast(Ed186_bin1_index, 3:end);
        Ed186_bin1_data_avg = nanmean(Ed186_bin1_data,1);
        AE_Ed186_bin1(m,:,n) = Ed186_bin1_data_avg;
        
        Lu174_bin1_index = find( ... 
          Lu174_cast(:,3)>=bins1(m)-bin1_step/2 & ... 
          Lu174_cast(:,3)<=bins1(m)+bin1_step/2);
        Lu174_bin1_data = Lu174_cast(Lu174_bin1_index, 3:end);
        Lu174_bin1_data_avg = nanmean(Lu174_bin1_data,1);
        AE_Lu174_bin1(m,:,n) = Lu174_bin1_data_avg;
    
    end
    
    % loop through bins2 and average the data:
    for m = 1:length(bins2)
    
        Ed186_bin2_index = find( ... 
          Ed186_cast(:,3)>=bins2(m)-bin2_step/2 & ... 
          Ed186_cast(:,3)<=bins2(m)+bin2_step/2);
        Ed186_bin2_data = Ed186_cast(Ed186_bin2_index, 3:end);
        Ed186_bin2_data_avg = nanmean(Ed186_bin2_data,1);
        AE_Ed186_bin2(m,:,n) = Ed186_bin2_data_avg;
        
        Lu174_bin2_index = find( ... 
          Lu174_cast(:,3)>=bins2(m)-bin2_step/2 & ... 
          Lu174_cast(:,3)<=bins2(m)+bin2_step/2);
        Lu174_bin2_data = Lu174_cast(Lu174_bin2_index, 3:end);
        Lu174_bin2_data_avg = nanmean(Lu174_bin2_data,1);
        AE_Lu174_bin2(m,:,n) = Lu174_bin2_data_avg;
    
    end
    
  end
   
  
% bin the combined measurements from all casts ...
 
  % ... create bin vectors:
  bin1_step = 0.20; bins1 = 0.10:bin1_step:20; 
  bin2_step = 1.00; bins2 = 0.50:bin2_step:20;    

  % ... pre-allocate space for the binned cast arrays:
  AE_Ed186_all_bin1 = nan(length(bins1),138);
  AE_Ed186_all_bin2 = nan(length(bins2),138);
  
  AE_Lu174_all_bin1 = nan(length(bins1),138);
  AE_Lu174_all_bin2 = nan(length(bins2),138);
  
  % ... re-order the complete data arrays by increasing depth:
  AE_Ed186_data = sortrows(AE_Ed186_data,3);
  AE_Lu174_data = sortrows(AE_Lu174_data,3);
  
% ... loop through bin and average the profile data:
  for m = 1:length(bins1)

      Ed186_bin1_index = find( ... 
        AE_Ed186_data(:,3)>=bins1(m)-bin1_step/2 & ... 
        AE_Ed186_data(:,3)<=bins1(m)+bin1_step/2);
      Ed186_bin1_data = AE_Ed186_data(Ed186_bin1_index, 3:end);
      Ed186_bin1_data_avg = nanmean(Ed186_bin1_data,1);
      AE_Ed186_all_bin1(m,:) = Ed186_bin1_data_avg;

      Lu174_bin1_index = find( ... 
        AE_Lu174_data(:,3)>=bins1(m)-bin1_step/2 & ... 
        AE_Lu174_data(:,3)<=bins1(m)+bin1_step/2);
      Lu174_bin1_data = AE_Lu174_data(Lu174_bin1_index, 3:end);
      Lu174_bin1_data_avg = nanmean(Lu174_bin1_data,1);
      AE_Lu174_all_bin1(m,:) = Lu174_bin1_data_avg;

  end

  % loop through bins2 and average the data:
  for m = 1:length(bins2)

      Ed186_bin2_index = find( ... 
        AE_Ed186_data(:,3)>=bins2(m)-bin2_step/2 & ... 
        AE_Ed186_data(:,3)<=bins2(m)+bin2_step/2);
      Ed186_bin2_data = AE_Ed186_data(Ed186_bin2_index, 3:end);
      Ed186_bin2_data_avg = nanmean(Ed186_bin2_data,1);
      AE_Ed186_all_bin2(m,:) = Ed186_bin2_data_avg;

      Lu174_bin2_index = find( ... 
        AE_Lu174_data(:,3)>=bins2(m)-bin2_step/2 & ... 
        AE_Lu174_data(:,3)<=bins2(m)+bin2_step/2);
      Lu174_bin2_data = AE_Lu174_data(Lu174_bin2_index, 3:end);
      Lu174_bin2_data_avg = nanmean(Lu174_bin2_data,1);
      AE_Lu174_all_bin2(m,:) = Lu174_bin2_data_avg;

  end

  
  
% interpolate values for binned profiles at the three selected wavelengths

% repeat for casts AE, AF, AG

% save output ...
%     1. array containing processed, unbinned data (all wavelengths)
%     2. array containing individual binned profiles (all wavelengths)
%     3. array containing combined binned profile (all wavelengths)
%     4. array containing processed, unbinned data (3 wavelengths)
%     5. array containing individual binned profiles (3 wavelengths)
%     6. array containing combined binned profile (3 wavelengths)


  % ... interpolate Es (Ed187) measurements to the 
%   Es_data = interp1(Ed186_times, Ed186_press, Lu174_times); 


  
  