function [s] = parse_spraysat2mat(satfpath, cal, missionname)


% Structures that are created:
% 
% c      -> contains calibration (or transformation) coefficients
% gps    -> info about GPS (obtaiend form G line)
% p      -> ctd profile data
% pH     -> DSD pH data
% s      -> final structure with the data 
% eng    -> engineering data
%
% create structure s with all of the information (gonna kind of be like
% netcdf)
%
% Created by: Yui Takeshita
% MBARI
% April 2019
%
% Modifications
% 
% 4/21/2019: YT
% - When there are different n points for p and pH, adjust the number of
% profiles in p.profinfo, and pH.profinfo. This number is used to 
% 8/7/2019: YT
% - Fixed bug for removing duplicate profiles. If there were multiple
% profiles, it was not working. ipdup was not getting indexed in the for
% loop. 
% 7/7/2020: YT
% - added if statement for 20703401 to deal with parsing issue with pH data
% - added if statement to deal with unpacked ctd data for first 2 profiles
%
% 9/13/2020: yT
% - added parsing information for engineering data

% choose filename
satfname = '0029.sat';
% missionname = '19A02901'; 
% satdirectory = ['\\atlas\carbonlab\SPRAY_PROCESSING\DATA\SAT_FILES\',missionname,'\'];

% *************************************************************************
% ****************** USER INPUTS FOR NOW **********************************
% *************************************************************************

% for 19302901, SBE10

% set pH voltage transformation
c.pH = [-2.6, 0.0001];

% for plotting purposes
% depname = '19A02901';
DSDname = 'SBE10';

% *************************************************************************1


% go line by line and put each line into a cell array

% fid = fopen([satdirectory,satfname]);
fid = fopen(satfpath);
tline = fgetl(fid);
tlines = cell(0,1);
while ischar(tline)
    tlines{end+1,1} = tline;
    tline = fgetl(fid);
end
fclose(fid);

disp([missionname,' data imported to Matlab...']);

% parse out each line using space delimiters
% can't have empty lines, so find empty lines
iempty = cellfun(@isempty, tlines);
tlines = tlines(~iempty); % delete all empty lines

trex = cellfun(@(x) textscan(x, '%s', 'delimiter', ' ', 'MultipleDelimsAsOne', 1), tlines);

% get the first value of the line
hdr = cellfun(@(x) x{1}, trex, 'uniformoutput', 0);

ipH = strcmp(hdr, 'pH'); % pH  
ip  = strcmp(hdr, 'p'); % profile (ctd)
iG  = strcmp(hdr, 'G'); % GPS 
iD  = strcmp(hdr, 'D'); % dive start (has n samps in this dive)
iPH = strcmp(hdr, 'PH'); % pH start line (has n samps in this dive)
iEP = strcmp(hdr, 'EP02'); % Engineering pump line 
iEC = strcmp(hdr, 'EC01'); % Engineering communication (SBD) line
iEF = strcmp(hdr, 'EF02'); % Engineering flight line
iEN = strcmp(hdr, 'EN01'); % Engineering navigation line

% *************************************************************************
% *********************** GPS LINES ***************************************
% *************************************************************************

% extract gps data
gps.divenum = cell2mat(cellfun(@(x) str2double(x{2}), trex(iG), 'UniformOutput', false));
gps.mstatus = cell2mat(cellfun(@(x) str2double(x{3}), trex(iG), 'UniformOutput', false));
gps.DTstring = cellfun(@(x) x(10:26), tlines(iG), 'UniformOutput', false);
gps.sdn     = datenum(gps.DTstring);
gps.fixflag = cell2mat(cellfun(@(x) str2double(x{8}), trex(iG), 'UniformOutput', false));
gps.latd    = cell2mat(cellfun(@(x) str2double(x{9}), trex(iG), 'UniformOutput', false));       % latitude degrees
gps.latm    = cell2mat(cellfun(@(x) str2double(x{10}), trex(iG), 'UniformOutput', false));      % latitude minutes
    gps.lat = gps.latd + gps.latm./60; % calculate decimal latitude
gps.lond    = cell2mat(cellfun(@(x) str2double(x{11}), trex(iG), 'UniformOutput', false));      % longitude degrees
gps.lonm    = cell2mat(cellfun(@(x) str2double(x{12}), trex(iG), 'UniformOutput', false));      % longitude minutes
    gps.lon = gps.lond + gps.lonm./60; % calculate decimal longitude

gps.TtoFix  = cell2mat(cellfun(@(x) str2double(x{13}), trex(iG), 'UniformOutput', false));
gps.nsat    = cell2mat(cellfun(@(x) str2double(x{14}), trex(iG), 'UniformOutput', false));
gps.minSNR  = cell2mat(cellfun(@(x) str2double(x{15}), trex(iG), 'UniformOutput', false));
gps.meanSNR = cell2mat(cellfun(@(x) str2double(x{16}), trex(iG), 'UniformOutput', false));
gps.maxSNR  = cell2mat(cellfun(@(x) str2double(x{17}), trex(iG), 'UniformOutput', false));
gps.HDOP    = cell2mat(cellfun(@(x) str2double(x{9}), trex(iG), 'UniformOutput', false));
% Skip some engineering status symbols: GPS health status, Wing Index
gps.lat     = cell2mat(cellfun(@(x) str2double(x{21}), trex(iG), 'UniformOutput', false));      % latitude decimal
gps.lon     = cell2mat(cellfun(@(x) str2double(x{22}), trex(iG), 'UniformOutput', false));      % longitude decimal

disp('Finished importing gps data');

% ================= Delete triplicate gps profiles ========================
% Added 4/8/2019 
% This assumes that the first two gps points are correct, and just removes
% all of the 3rd, 4th, ... gps fixes. Not the perfect way to do it, but
% will keep on doing it until it stops working.

% find triplicates or greater
U = unique(gps.divenum);
nU = NaN(size(U));
for i = 1:length(U)
    nU(i) = length(find(gps.divenum==U(i)));
end

imultgps = find(nU > 2); % get indicies for where there are more than 2 gps fixes

if(~isempty(imultgps))
    disp(['Multiple gps fixes foun on profile ',num2str(U(imultgps)'),'... removed']);
        
    iremove = false(size(gps.divenum)); % presize boolean array

    for i = 1:length(imultgps)
        ngps = nU(imultgps(i)); % how many fixes there were
        igpsstart = find(gps.divenum==U(imultgps(i)),1,'first');

        iremove(igpsstart+2:igpsstart+2+ngps-3) = true; 
    end

    gpsfields = fieldnames(gps);

    for g = 1:length(gpsfields)
        gps.(gpsfields{g})(iremove) = [];
    end

end

clear gpsfields g i iremove igpsstart ngps imultgps nU U;

% =========================================================================

% *************************************************************************
% *********************** ENG LINES ***************************************
% *************************************************************************
% Added 9/13/2020 YT

% Engineering Communication Line
eng.divenum_ec      = cell2mat(cellfun(@(x) str2double(x{2}), trex(iEC), 'UniformOutput', false));  % divenum
eng.ntries          = cell2mat(cellfun(@(x) str2double(x{3}), trex(iEC), 'UniformOutput', false));  % # tries sending an SBD message during the LAST surfacing
eng.nsent           = cell2mat(cellfun(@(x) str2double(x{4}), trex(iEC), 'UniformOutput', false));  % # of SBD messages sent during the LAST surfacing
eng.sbdi_stat       = cell2mat(cellfun(@(x) str2double(x{5}), trex(iEC), 'UniformOutput', false));  % status of the last attempt to send a message
eng.sbdi_shore_stat = cell2mat(cellfun(@(x) str2double(x{6}), trex(iEC), 'UniformOutput', false));  % status of parsing the last shore command; hex
eng.t_sbd           = cell2mat(cellfun(@(x) str2double(x{7}), trex(iEC), 'UniformOutput', false));  % time (sec) that the modem actively tried to send a message
eng.wing            = cell2mat(cellfun(@(x) str2double(x{8}), trex(iEC), 'UniformOutput', false));  % wing used to send the sbd message

% Engineering Flight
eng.divenum_ef    = cell2mat(cellfun(@(x) str2double(x{2}), trex(iEF), 'UniformOutput', false));    % divenum
eng.navg          = cell2mat(cellfun(@(x) str2double(x{3}), trex(iEF), 'UniformOutput', false));    % # of 8s samples avgeraged for each datum in the sbd message
eng.psurf         = cell2mat(cellfun(@(x) str2double(x{4}), trex(iEF), 'UniformOutput', false));    % surface pressure counts
eng.zmax          = cell2mat(cellfun(@(x) str2double(x{5}), trex(iEF), 'UniformOutput', false));    % max depth (dBar) recorded in the bototm turn
eng.pitch         = cell2mat(cellfun(@(x) str2double(x{6}), trex(iEF), 'UniformOutput', false));    % pitch angle (degrees)
eng.alt_min_depth = cell2mat(cellfun(@(x) str2double(x{7}), trex(iEF), 'UniformOutput', false));    % altimeter minimum depth to the bototm (m). if 80, no bottom was detected. if >80, then not sampled
eng.adp_intens    = cell2mat(cellfun(@(x) str2double(x{8}), trex(iEF), 'UniformOutput', false));    % adp altimeter bototm intensity level (dB); =0 for Tritech
eng.zat_alt       = cell2mat(cellfun(@(x) str2double(x{9}), trex(iEF), 'UniformOutput', false));    % spray depth (dBar) at the time of the altimeter reading
eng.roll_err      = cell2mat(cellfun(@(x) str2double(x{10}), trex(iEF), 'UniformOutput', false));   % integrated roll error (degrees) form the PI loop
eng.exc_stat      = cell2mat(cellfun(@(x) str2double(x{11}), trex(iEF), 'UniformOutput', false));   % exception stats bit setting (Hex); decoder below

% Exc_stat decoding key from spray_SAT_format.doc
% 
% EXC_STATUS   Let b0 = bit[0] = Lsbit, b15 = bit[15] = Msbit.
% 	b0 = Pump Recovery was required.            b8 = CF1 Close File error.
% 	b1 = Drop Weight Activated.                 b9 = High pump amps @ 50 m safety pump.
% 	b2 = Pressure>20 m ‘at the surface.’        b10= Press=0 (broken) ‘at the surface.’
% 	b3 = Depth>1500 m detected.                 b11= No SBD sent in last 48 hrs.
% 	b4 = The Altimeter triggered turn-around.	b12= Cannot overcome the current.
% 	b5 = Backed off the surface.                b13= Spurious Reset Detected.
% 	b6 = Extra pumps req’d at end of ascent. 	b14= Altimeter reading is from an ADP.
% 	b7 = Took >700 s to leave the surface.      b15= Not Used.
% 
%   b12=1 means that although current correction is applied, the current is too strong to maintain course.
%   b14 is used to help interpret the altimeter reading in the ‘E’ line.


% Engineering Navigation Line
eng.divenum_en = cell2mat(cellfun(@(x) str2double(x{2}), trex(iEN), 'UniformOutput', false));       % divenum
eng.dr_x       = cell2mat(cellfun(@(x) str2double(x{3}), trex(iEN), 'UniformOutput', false));       % dead-reckoning x distance traveled during the dive (m)
eng.dr_y       = cell2mat(cellfun(@(x) str2double(x{4}), trex(iEN), 'UniformOutput', false));       % dead-reckoning y distance traveled during the dive (m)
eng.wlat       = cell2mat(cellfun(@(x) str2double(x{5}), trex(iEN), 'UniformOutput', false));       % waypoint latitude (dec deg)
eng.wlon       = cell2mat(cellfun(@(x) str2double(x{6}), trex(iEN), 'UniformOutput', false));       % waypoint longitude (dec deg)
eng.tleave     = cell2mat(cellfun(@(x) str2double(x{7}), trex(iEN), 'UniformOutput', false));       % time at start of dive to leave the surface (seconds)
eng.tend       = cell2mat(cellfun(@(x) str2double(x{8}), trex(iEN), 'UniformOutput', false));       % time at end of dive since p<2 m until GPS is on (seconds)
eng.tslow      = cell2mat(cellfun(@(x) str2double(x{9}), trex(iEN), 'UniformOutput', false));       % time off of surface where |w|<5 cm/s (seconds)
eng.des_head   = cell2mat(cellfun(@(x) str2double(x{9}), trex(iEN), 'UniformOutput', false));       % what we tried to steer 
eng.heading    = cell2mat(cellfun(@(x) str2double(x{9}), trex(iEN), 'UniformOutput', false));       % T = true for heading (v0608, 0610), M = magnetic for heading (< 0608)

% Engineering Pump Line
eng.divenum_ep  = cell2mat(cellfun(@(x) str2double(x{2}), trex(iEP), 'UniformOutput', false));      % divenum
eng.maxdepth    = cell2mat(cellfun(@(x) str2double(x{3}), trex(iEP), 'UniformOutput', false));      % max depth (dbar)
eng.avgbattvolt = cell2mat(cellfun(@(x) str2double(x{4}), trex(iEP), 'UniformOutput', false));      % average batt volt
eng.avgamp      = cell2mat(cellfun(@(x) str2double(x{5}), trex(iEP), 'UniformOutput', false));      % average amp draw during bottom turn
eng.nbadamp     = cell2mat(cellfun(@(x) str2double(x{6}), trex(iEP), 'UniformOutput', false));      % # of times that pump current exceeded the max allowed
eng.tmaxamp     = cell2mat(cellfun(@(x) str2double(x{7}), trex(iEP), 'UniformOutput', false)).*10;  % time when max amp occured during bottom turn; sent as s/10, so convert to seconds
eng.maxamps     = cell2mat(cellfun(@(x) str2double(x{8}), trex(iEP), 'UniformOutput', false));      % max pump current recorded in bottom turn
eng.vacuum      = cell2mat(cellfun(@(x) str2double(x{9}), trex(iEP), 'UniformOutput', false));      % -Hg inside Spray pressure case
eng.tpump1      = cell2mat(cellfun(@(x) str2double(x{10}), trex(iEP), 'UniformOutput', false)).*10; % time pumped until Spray was neutral; sent as s/10, so convert to seconds
eng.tpump2      = cell2mat(cellfun(@(x) str2double(x{11}), trex(iEP), 'UniformOutput', false)).*10; % remaining pump time; includes ascent + surface pump time; sent as s/10, so convert to seconds
    eng.tpumptot = eng.tpump1 + eng.tpump2;
eng.venttime    = cell2mat(cellfun(@(x) str2double(x{12}), trex(iEP), 'UniformOutput', false));     % time air vent ran to purge air from the hydraulics; seconds
eng.airflag     = cell2mat(cellfun(@(x) str2double(x{13}), trex(iEP), 'UniformOutput', false));     % air_flag; 1 if LLD detected air; =0 if detected oil.

% =========================================================================

% *************************************************************************
% ******************** TRANFORMATION COEFFICIENTS *************************
% *************************************************************************

% Get calibration coefficients

% index for pressure, temperautre, salinity, optical, and DO cal
% gain/offset for voltage transformation
%
iCP = strncmpi(hdr, 'CP', 2); % pressure 
iCT = strncmpi(hdr, 'CT', 2); % Temperature
iCS = strncmpi(hdr, 'CS', 2); % salinity
iCO = strncmpi(hdr, 'CO', 2); % Optical 
iCX = strncmpi(hdr, 'CX', 2); % O2 sensor 

% get coefficient list
c.P = cell2mat(cellfun(@(x) str2double(x(2:end))', trex(iCP), 'UniformOutput', false));
c.T = cell2mat(cellfun(@(x) str2double(x(2:end))', trex(iCT), 'UniformOutput', false));
c.S = cell2mat(cellfun(@(x) str2double(x(2:end))', trex(iCS), 'UniformOutput', false));
c.O = cell2mat(cellfun(@(x) str2double(x(2:end))', trex(iCO), 'UniformOutput', false));
c.X = cell2mat(cellfun(@(x) str2double(x(2:end))', trex(iCX), 'UniformOutput', false));

disp('Finished importing calibration data');

% *************************************************************************
% *************************** PROFILE DATA ********************************
% *************************************************************************

% extract the cell array for the corresponding variables
pcells  = trex(ip);
pHcells = trex(ipH);

% check to see if ctd data is complete. length should be 8.
pcell_length = cellfun(@length,pcells);
ifix = find(pcell_length < 8); % find all data lines with fewer than 8 columns
% loop through and add {'-999'} for each missing value
for i = 1:length(ifix)
    ncell2add = 8 - pcell_length(i); 
    [cell2append{1:ncell2add}] = deal('-999'); % make
    pcells{ifix(i)} = cell(vertcat(pcells{ifix(i)},cell2append(:)));
    clear cell2append ncell2add
end

% extract D and PH (profile 
p.profinfo  = cell2mat(cellfun(@(x) str2double(x(2:3))', trex(iD), 'UniformOutput', false));
pH.profinfo = cell2mat(cellfun(@(x) str2double(x(2:3))', trex(iPH), 'UniformOutput', false));

% check to see if prof info for p and pH are the same. If not, display
% error
if(~isequal(p.profinfo, pH.profinfo))
    disp('!!!!!!!!!!!pH and p dive num and or npts are not equal!!!!!!!!!');
end

% extract profile data (raw, pretransformed data)
p.rawdata  = cell2mat(cellfun(@(x) str2double(x(2:end))', pcells, 'UniformOutput', false));
pH.rawdata = cell2mat(cellfun(@(x) str2double(x(2:end))', pHcells, 'UniformOutput', false));

p.rawdata(p.rawdata==-999) = NaN; % make the missing data into NaNs. 


disp('Finished importing profile data');

% ====================== REMOVE DUPLICATE PROFILES ========================
% Added 4/8/2019, since multiple profiles for pH and p were coming in for
% 190402901
%
% find indicies for duplicates for ctd
[~,I] = unique(p.profinfo(:,1),'first');
ipdup = 1:length(p.profinfo(:,1));
    ipdup(I) = [];
if(~isempty(ipdup))
    
    % delete ip
    iremove = false(size(p.rawdata,1),1);
    
    % loop through if there are multiple duplicate profiles
    for i = 1:length(ipdup)
        disp(['Removing prof ',num2str(p.profinfo(ipdup(i),1)),' for ctd']);
        irowstart = sum(p.profinfo(1:ipdup(i)-1,2))+1;
        irowend = irowstart+p.profinfo(ipdup(i),2)-1;
        % make sections of the boolean array 'true' where you want to
        % remove
        iremove(irowstart:irowend) = true;      
    end
%     remove profinfo and rawdata
    p.rawdata(iremove,:) = [];
    p.profinfo(ipdup,:) = [];    
end

clear I i iremove irowstart irowend

% find indicies for duplicates for pH
[~,I] = unique(pH.profinfo(:,1),'first');
ipHdup = 1:length(pH.profinfo(:,1));
    ipHdup(I) = [];
if(~isempty(ipHdup))
    
    % delete ipH
    iremove = false(size(pH.rawdata,1),1);
    
    % loop through if there are multiple duplicate profiles
    for i = 1:length(ipHdup)
        disp(['Removing prof ',num2str(pH.profinfo(ipHdup(i),1)),' for pH']);
        irowstart = sum(pH.profinfo(1:ipHdup(i)-1,2))+1;
        irowend = irowstart+pH.profinfo(ipHdup(i),2)-1;
        % make sections of the boolean array 'true' where you want to
        % remove
        iremove(irowstart:irowend) = true;      
    end
    % remove profinfo and rawdata
    pH.rawdata(iremove,:) = [];
    pH.profinfo(ipHdup,:) = [];    
end

clear I i iremove irowstart irowend


% =========================================================================




% *************************************************************************
% *************************** CONVERT DATA ********************************
% *************************************************************************

% now transform raw values to 
% p (ctd) prof data
idivenum = 1;
ipackind = 2;
ipres    = 3;
itc      = 4;
ipsal    = 5;
iopt     = 6;
imoxy    = 7;
% pH indicies; divenum and packind are the same
ivrs     = 10;
ivrs_std = 11;
ivk      = 12;
ivk_std  = 13;
iIk      = 14;
iIb      = 15;
ipH      = 16; % 
% Parsing for 20703401 pH data was off. Bandaid solution, delete this if
% statment in the future. 7/7/2020 YT
if(strcmp(missionname,'20703401'))
    ivrs     = 11;
    ivrs_std = 12;
    ivk      = 13;
    ivk_std  = 14;
    iIk      = 15;
    iIb      = 10;   
end
% check to make sure all profs have the same number of data points
iprofmissdata = find(p.profinfo(:,2) ~= pH.profinfo(:,2)); 

% pad data to the top (end) of the profile
for i = 1:length(iprofmissdata)
    disp(['Profile ',num2str(iprofmissdata(i)),' adjusted so pH and ctd have same number of points: p = '...
        ,num2str(p.profinfo(iprofmissdata(i),2)),', pH = ',num2str(pH.profinfo(iprofmissdata(i),2))]);
    
    % find which one is smaller
    if(p.profinfo(iprofmissdata(i),2) < pH.profinfo(iprofmissdata(i),2))
        % ctd has fewer data points
        % find index of last incidence of this profile
        iprofend = find(p.rawdata(:,1)==iprofmissdata(i),1,'last');
        % get number of rows we need to insert; 
        ddata = p.profinfo(iprofmissdata(i),2) - pH.profinfo(iprofmissdata(i),2);
        tempmat = NaN(abs(ddata),size(p.rawdata,2));
        tempmat(:,1) = iprofmissdata(i);
        % insert row
        p.rawdata = [p.rawdata(1:iprofend,:); tempmat; p.rawdata(iprofend+1:end,:)];
    else
        % pH has fewer data points
         % find index of last incidence of this profile
        iprofend = find(pH.rawdata(:,1)==iprofmissdata(i),1,'last');
        % get number of rows we need to insert; 
        ddata = pH.profinfo(iprofmissdata(i),2) - p.profinfo(iprofmissdata(i),2);
        tempmat = NaN(abs(ddata),size(pH.rawdata,2));
        tempmat(:,1) = iprofmissdata(i);
        % insert rows
        pH.rawdata = [pH.rawdata(1:iprofend,:); tempmat; pH.rawdata(iprofend+1:end,:)];       
%         return; 
    end
    
    
    
end
% 
% p.profinfo(iprofmissdata,:) = [];
% pH.profinfo(iprofmissdata,:) = [];

% data = horzcat(p.rawdata(p.rawdata(:,1)~=iprofmissdata,:), pH.rawdata(pH.rawdata(:,1)~=iprofmissdata,:));

data = horzcat(p.rawdata, pH.rawdata);


% transform data
data(:,ipres)     = data(:,ipres).*c.P(2) + c.P(1);
data(:,itc)     = data(:,itc).*c.T(2) + c.T(1);
data(:,ipsal)     = data(:,ipsal).*c.S(2) + c.S(1) + cal.psal_offset;
data(:,iopt)      = (data(:,iopt) .*c.O(2) + c.O(1)) .* c.O(4); % Optical sensors require a 2-step linear cal. First transformation to Volts, then to Sci units. 
data(:,imoxy)       = data(:,imoxy)  .*c.X(2) + c.X(1);
data(:,ivrs)      = data(:,ivrs) .*c.pH(2) + c.pH(1);
data(:,ivrs_std)  = data(:,ivrs_std) .*c.pH(2) + c.pH(1);
data(:,ivk)       = data(:,ivk) .*c.pH(2) + c.pH(1);
data(:,ivk_std)   = data(:,ivk_std) .*c.pH(2) + c.pH(1);
data(:,iIk)       = data(:,iIk) .*c.pH(2) + c.pH(1);
    data(:,iIk)   = data(:,iIk) ./ 0.00015; % in nA; 150uV / nA 
data(:,iIb)       = data(:,iIb) .*c.pH(2) + c.pH(1);
    data(:,iIb)   = data(:,iIb) ./ 0.001; % in nA; 1mV / nA
% make pH data NaN for when it was off
% ipHoff = data(:,ivrs) > -0.8; % crude way fo doing it
% ipHoff = data(:,ivrs) > -0.6 | data(:,ivrs) < -1.2; % crude way fo doing it
% 
% data(ipHoff,ivrs) = NaN;
% data(ipHoff,ivrs_std) = NaN;
% data(ipHoff,ivk) = NaN;
% data(ipHoff,ivk_std) = NaN;
% data(ipHoff,iIk) = NaN;
% data(ipHoff,iIb) = NaN;

% calculate pH; get total scale
[~,pH_calculated] = phcalc_jp(data(:,ivrs), data(:,ipres), data(:,itc), data(:,ipsal), cal.k0, cal.k2, cal.fP);
% add pH to the data amatrix
data = horzcat(data,pH_calculated);

disp('Finished converting profile data, calculated pH');

% *************************************************************************
% ******************* PARSE DATA TO FINAL STRUCTURE ***********************
% *************************************************************************


% *********************** GPS INFO AT DEPLOYMENT ************************

ideploy = gps.divenum == 0; % the gps fix upon deployment

if(length(find(ideploy)) > 1)
    disp('!!!!!!! MORE THAN 1 DEPLOYMENT GPS FIX !!!!!!!!!!!!!!!!!!!!!!');
end

s.sdn_start = gps.sdn(ideploy);
s.lat_start = gps.lat(ideploy);
s.lon_start = gps.lon(ideploy);


% ************************* START PROFILE DATA ****************************

% presize matrix. Get maximum npts from the highest dive number and npts
% from the profinfos.
msize = [max(p.profinfo(:,2)), max(p.profinfo(:,1))]; % might have to adjust this to choose max of p and pH, if p ever has fewer points than pH. 

% gps data (matrix); contains start and end. Average this at the end to get
% array (to get one corresponding gps per profile data)
s.sdn_  = NaN(2,msize(2));
s.lat_  = NaN(2,msize(2));
s.lon_  = NaN(2,msize(2));

% profile data (matrix)
s.pres = NaN(msize);
s.tc   = NaN(msize);
s.psal = NaN(msize);
s.opt  = NaN(msize); % optical output from seapoint
s.moxy_ml = NaN(msize);
s.vrs  = NaN(msize);
s.vrs_std = NaN(msize);
s.vk   = NaN(msize);
s.vk_std = NaN(msize);
s.Ik   = NaN(msize);
s.Ib   = NaN(msize); 
s.pHin = NaN(msize); % in situ pH

% start loop
% for i = 1:msize(2)
for i = 1:size(p.profinfo,1)
    
    % get divenum and npts
    divenum = p.profinfo(i,1);
    npts = p.profinfo(i,2);
    
    % ***!!!!!! STILL NEED TO ADD CODE TO SEARCH FOR GOOD AND BAD FIXES !!*
    
    % extract gps data
    igps = gps.divenum == divenum;
    if(length(find(igps)) ~= 2)
        disp(['!!!!!!!!!!! NOT 2 GPS POINTS FOR THIS DIVE (',num2str(divenum),')!!!!']);
    else
        s.sdn_(1,divenum) = gps.sdn(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
        s.sdn_(2,divenum) = gps.sdn(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive
        % lat
        s.lat_(1,divenum) = gps.lat(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
        s.lat_(2,divenum) = gps.lat(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive
        % lon
        s.lon_(1,divenum) = gps.lon(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
        s.lon_(2,divenum) = gps.lon(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive

        % extract matrix data

    %     % get start and end index for this dive
    %     if(i==1)
    %         istart = 1;
    %     else
    %         istart = sum(p.profinfo(1:i-1,2))+1;
    %     end
    %     iend = sum(p.profinfo(1:i,2));
        % modified it so it just searches for the profile numbers. 
        istart = find(data(:,1)==divenum,1,'first');
        iend = find(data(:,1)==divenum,1,'last');
        npts = iend-istart+1; 

        s.pres(1:npts,divenum) = flipud(data(istart:iend,ipres)); % flip so it writes shallow to deep
        s.tc(1:npts,divenum)   = flipud(data(istart:iend,itc)); % flip so it writes shallow to deep
        s.psal(1:npts,divenum) = flipud(data(istart:iend,ipsal)); % flip so it writes shallow to deep
        s.opt(1:npts,divenum)  = flipud(data(istart:iend,iopt)); % flip so it writes shallow to deep
        s.moxy_ml(1:npts,divenum) = flipud(data(istart:iend,imoxy)); % flip so it writes shallow to deep
        s.vrs(1:npts,divenum)  = flipud(data(istart:iend,ivrs)); % flip so it writes shallow to deep
        s.vrs_std(1:npts,divenum)  = flipud(data(istart:iend,ivrs_std)); % flip so it writes shallow to deep
        s.vk(1:npts,divenum)   = flipud(data(istart:iend,ivk));  % flip so it writes shallow to deep
        s.vk_std(1:npts,divenum)   = flipud(data(istart:iend,ivk_std));  % flip so it writes shallow to deep
        s.Ik(1:npts,divenum)   = flipud(data(istart:iend,iIk)); % flip so it writes shallow to deep
        s.Ib(1:npts,divenum)   = flipud(data(istart:iend,iIb)); % flip so it writes shallow to deep
        s.pHin(1:npts,divenum) = flipud(data(istart:iend,ipH)); % flip so it writes shallow to deep

    end
    
    
end

% average gps data to get corresponding match for the profile
s.sdn = nanmean(s.sdn_,1);
s.lat = nanmean(s.lat_,1);
s.lon = nanmean(s.lon_,1);

% convert oxygen to umol/kg
% Apply salt correction. Use CTD temp; check with window temp later
Ts    = log((298.15 - s.tc) ./ (273.15 + s.tc)); % Scaled T
% Descending powers for Matlab polyfit
SolB  = [-8.17083e-3 -1.03410e-2 -7.37614e-3 -6.24523e-3];
C0    = -4.88682e-7;
% salinity correction
Scorr = exp((s.psal).*polyval(SolB,Ts) + C0.*(s.psal.^2));
% calculate moxy_ml with salinity correction
s.moxy_ml = s.moxy_ml.*Scorr; % apply salinity correction
s.moxy = s.moxy_ml .* 44.6596; % uncalibrated oxygen in uM
s.pdens = sw_pden(s.psal, s.tc, s.pres, 0) - 1000; % in kg/L
s.doxy = s.moxy ./ ((1000+s.pdens)./1000) .* cal.O2gain; % umol/kg
s.o2satper = s.doxy ./ calcO2sat(s.tc, s.psal) .* 100; % Saturation %

% calculate depth
s.depth = sw_dpth(s.pres, s.lat);

% apply calibration to chl-a 7/19/2020 YT
s.opt_raw = s.opt; % save original raw output
s.opt = (s.opt + cal.opt_offset).* cal.opt_gain;

% start flag assignment

% 0 = good
% 1 = missing or uninspcted
% 4 = questionable
% 8 = bad

% lat/lon QC
s.position_QC = zeros(size(s.lat));
s.sdn_QC = zeros(size(s.sdn));

% default them all to 0
QFvars = {'pres_QC', 'tc_QC', 'psal_QC', 'opt_QC', 'moxy_QC', 'doxy_QC', 'pdens_QC', 'pHin_QC', 'depth_QC', 'o2satper_QC', 'opt_QC'};

for v = 1:length(QFvars)
    % default them all to 0
    s.(QFvars{v}) = zeros(size(s.pres));
    
    % Nans are missing values. make them QF = 1; % 1:end-3 cuts off _QC
    % from the string. 
    s.(QFvars{v})(isnan(s.(QFvars{v}(1:end-3)))) = 1;
end

% a little extra QC for pH and DOXY
for i = 1:size(s.doxy_QC,2)
    % go through each profile, and make the bottom 5 measurments bad for
    % oxy and 2 for pH. Oxy has a slower response time, I think. 
    ilastO2 = find(isfinite(s.doxy(:,i)),1,'last');
    s.doxy_QC(ilastO2-5:ilastO2,i) = 8;
    s.moxy_QC(ilastO2-5:ilastO2,i) = 8;
    s.o2satper_QC(ilastO2-5:ilastO2,i) = 8;
    % now pH
    ilastpH = find(isfinite(s.pHin(:,i)),1,'last');
    s.pHin_QC(ilastpH-2:ilastpH,i) = 8;
end

% % final QC for pH, make first 5 profiles all QF = 8, due to comm issues.
% trex = s.pHin_QC(:,1:5);
% trex(trex~=1) = 8; 
% s.pHin_QC(:,1:5) = trex; 

% run manual QC mfile if it exists
if(isfile(sprintf('%s/QCrt_%s.m',strippath(satfpath,2),missionname)))
    run(sprintf('%s/QCrt_%s.m',strippath(satfpath,2),missionname));
end



disp('Finished parsing data into s structure');

% add mission name
s.depID = missionname; 

% do canyonB
s = calc_carb_spray(s);

% save mat file
save([missionname,'_sat.mat'], 's', 'c', 'gps', 'pH', 'p', 'data', 'eng');
% copyfile([missionname,'_sat.mat'], 'G:\My Drive\Summer Internships\2020\Anna Hughes MBARI Summer\Glider Data\Current_deps\');
movefile([missionname,'_sat.mat'], strippath(satfpath,2), 'f');
% save kml file for spray coordinates

kmlwritepoint(missionname, s.lat_(2,:), s.lon_(2,:));
% copyfile([missionname,'.kml'], 'G:\My Drive\Summer Internships\2020\Anna Hughes MBARI Summer\Glider Data\Current_deps\');
movefile([missionname,'.kml'], strippath(satfpath,2), 'f');



return; 
