function [cast, dcast] = rbrctd(filnam)
% [cast, dcast] = rbrctd(filnam)
% 
% RBRCTD:  Convert RBR CTD data file to a usable Matlab structure
%  INPUT:  filnam - filename of RBR .mat file generated in Ruskin 
%  (File -> Export dataset -> Matlab MAT)
%  OUTPUT: cast - all recorded data for cast
%          dcast - downcast only (based on maximum pressure reading)
%
% Read Ruskin-generated .mat file containing data from several CTD casts 
% and segment data into a structure by cast.  Input is filename (within 
% single quotes), outputs are structures containing data indexed by cast, 
% one for the entire cast, the other for downcast only.
%
% Currently, no error checking is performed.
%
% This function likely only works for thresholded data collected at 6 Hz.
% Data columns may vary depending on what sensors are installed
%
% D. Nowacki (nowacki@uw.edu) September 2011
% 2014-07-02 fix bug in how cast times are computed.

load(filnam)

evnt = find(isnan(RBR.data(:,1))); % get locations of beginning of casts
% this is not documented by RBR, but casts within a single file are 
% delimited by NaNs.
evnt(length(evnt)+1) = length(RBR.data)+1;
diffs = diff(evnt);
% 61 seems to be the number of samples that are collected when the
% instrument is initially deployed).  This threshold also throws out the 6
% samples that are collected between casts.
evnt1 = evnt(find(diffs > 61)) + 1;
% evnt2 = evnt(find(diffs > 61)+1);
len = diffs(find(diffs > 61));

% set locations for beginning and end of each cast
for n = 1:length(evnt1)
    cast(n).start = evnt1(n);
%     cast(n).end = evnt2(n) - 1;
    cast(n).end = evnt1(n) - 2 + len(n);
end

for n = 1:length(cast)
   % compile data from entire cast
   cast(n).cond = RBR.data(cast(n).start:cast(n).end, 1);
   cast(n).temp = RBR.data(cast(n).start:cast(n).end, 2);
   cast(n).pres = RBR.data(cast(n).start:cast(n).end, 3);
   cast(n).volt1 = RBR.data(cast(n).start:cast(n).end, 4);
   cast(n).volt2 = RBR.data(cast(n).start:cast(n).end, 5);
   cast(n).depth = RBR.data(cast(n).start:cast(n).end, 6);
   cast(n).sal = RBR.data(cast(n).start:cast(n).end, 7);
   cast(n).speccond = RBR.data(cast(n).start:cast(n).end, 8);
   cast(n).densanom = RBR.data(cast(n).start:cast(n).end, 9);
   cast(n).sos = RBR.data(cast(n).start:cast(n).end, 10);
   
   % assign start time
   if isfield(RBR, 'sampletimes')
       % RBR.sampletimes is a new variable output by recent versions of
       % Ruskin
       % set average cast time based on total times contained in cast
       cast(n).time = nanmean(datenum(RBR.sampletimes(cast(n).start:cast(n).end), 'dd/mm/yyyy HH:MM:SS.FFF PM'));
   else
       warning('Time assigned to casts might be incorrect. Re-export .mat file with a recent version of Ruskin and re-process with rbrctd.m.');
       % this method is crude and probably incorrect.
       [~, ap] = min(abs(cell2mat(RBR.events(:,3)) - cast(n).start));
       cast(n).time = datenum(RBR.events(ap,2), 'dd/mm/yyyy HH:MM:SS PM');
   end
      
   [~, maxd] = max(cast(n).depth);
   
   % compile data from downcast only, using maximum depth  
   dcast(n).cond = cast(n).cond(1:maxd);
   dcast(n).temp = cast(n).temp(1:maxd);
   dcast(n).pres = cast(n).pres(1:maxd);
   dcast(n).volt1 = cast(n).volt1(1:maxd);
   dcast(n).volt2 = cast(n).volt2(1:maxd);
   dcast(n).depth = cast(n).depth(1:maxd);
   dcast(n).sal = cast(n).sal(1:maxd);
   dcast(n).speccond = cast(n).speccond(1:maxd);
   dcast(n).densanom = cast(n).densanom(1:maxd);
   dcast(n).sos = cast(n).sos(1:maxd);
   dcast(n).time = cast(n).time;
end

end