function [s] = parse_spraysat2mat(satfpath, cal)


% 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 
%
% 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 

% choose filename
satfname = '0029.sat';
missionname = '19402901'; 
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 = '19402901';
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([satfname,' 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)

% *************************************************************************
% *********************** 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;

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



% *************************************************************************
% ******************** 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);

% 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));

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))
    disp(['Removing prof ',num2str(p.profinfo(ipdup,1)),' for ctd']);
    % delete ip
    iremove = false(size(p.rawdata,1));
    
    % loop through if there are multiple duplicate profiles
    for i = 1:length(ipdup)
        irowstart = sum(p.profinfo(1:ipdup-1,2))+1;
        irowend = irowstart+p.profinfo(ipdup,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))
    disp(['Removing prof ',num2str(pH.profinfo(ipHdup,1)),' for pH']);
    % delete ipH
    iremove = false(size(pH.rawdata,1));
    
    % loop through if there are multiple duplicate profiles
    for i = 1:length(ipHdup)
        irowstart = sum(pH.profinfo(1:ipHdup-1,2))+1;
        irowend = irowstart+pH.profinfo(ipHdup,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; % 

% 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);
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
pH_calculated = phcalc_jp(data(:,ivrs), data(:,ipres), data(:,itc), data(:,ipsal), cal.k0, cal.k2, cal.fP);
% add pH to the dat 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 SBE63
s.moxy = NaN(msize);
s.vrs  = NaN(msize);
s.vk   = NaN(msize);
s.Ik   = NaN(msize);
s.Ib   = NaN(msize); 
s.pHin = NaN(msize); % in situ pH

% start loop
for i = 1:msize(2)
    
    % get divenum and npts
    divenum = p.profinfo(i,1);
%     npts = p.profinfo(i,2); % commented 4/21/2019, and redefined below. 
    
    % ***!!!!!! 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),')!!!!']);
    end
    s.sdn_(1,i) = gps.sdn(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
    s.sdn_(2,i) = gps.sdn(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive
    % lat
    s.lat_(1,i) = gps.lat(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
    s.lat_(2,i) = gps.lat(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive
    % lon
    s.lon_(1,i) = gps.lon(igps & gps.mstatus == 1); % check for mission status; 1 = start of dive
    s.lon_(2,i) = gps.lon(igps & gps.mstatus == 2); % check for mission status; 2 = end of dive

 
% ================== COMMENTED OUT FOR BETTER VERSION 4/21/2019 ===========
%    % 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. 
% 
%         
%     s.pres(1:npts,i) = flipud(data(istart:iend,ipres)); % flip so it writes shallow to deep
%     s.tc(1:npts,i)   = flipud(data(istart:iend,itc)); % flip so it writes shallow to deep
%     s.psal(1:npts,i) = flipud(data(istart:iend,ipsal)); % flip so it writes shallow to deep
%     s.opt(1:npts,i)  = flipud(data(istart:iend,iopt)); % flip so it writes shallow to deep
%     s.moxy(1:npts,i) = flipud(data(istart:iend,imoxy)); % flip so it writes shallow to deep
%     s.vrs(1:npts,i)  = flipud(data(istart:iend,ivrs)); % flip so it writes shallow to deep
%     s.vk(1:npts,i)   = flipud(data(istart:iend,ivk));  % flip so it writes shallow to deep
%     s.Ik(1:npts,i)   = flipud(data(istart:iend,iIk)); % flip so it writes shallow to deep
%     s.Ib(1:npts,i)   = flipud(data(istart:iend,iIb)); % flip so it writes shallow to deep
%     s.pHin(1:npts,i) = flipud(data(istart:iend,ipH)); % flip so it writes shallow to deep
% =========================================================================
    % extract matrix data
    % better way of indexing. 4/21/2019 YT
    iprof = data(:,1)==i;
    npts = length(find(iprof));
    
    % Write into 's' structure from data. 
    s.pres(1:npts,i) = flipud(data(iprof,ipres)); % flip so it writes shallow to deep
    s.tc(1:npts,i)   = flipud(data(iprof,itc)); % flip so it writes shallow to deep
    s.psal(1:npts,i) = flipud(data(iprof,ipsal)); % flip so it writes shallow to deep
    s.opt(1:npts,i)  = flipud(data(iprof,iopt)); % flip so it writes shallow to deep
    s.moxy(1:npts,i) = flipud(data(iprof,imoxy)); % flip so it writes shallow to deep
    s.vrs(1:npts,i)  = flipud(data(iprof,ivrs)); % flip so it writes shallow to deep
    s.vk(1:npts,i)   = flipud(data(iprof,ivk));  % flip so it writes shallow to deep
    s.Ik(1:npts,i)   = flipud(data(iprof,iIk)); % flip so it writes shallow to deep
    s.Ib(1:npts,i)   = flipud(data(iprof,iIb)); % flip so it writes shallow to deep
    s.pHin(1:npts,i) = flipud(data(iprof,ipH)); % flip so it writes shallow to deep
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
s.dens = sw_dens(s.psal, s.tc, s.pres) - 1000; % in kg/L
s.doxy = s.moxy.*31.998 ./ ((1000+s.dens)./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);

% 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', 'dens_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; 


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

% save mat file
save([depname,'_sat.mat'], 's', 'c', 'gps', 'pH', 'p', 'data');
% save kml file for spray coordinates

kmlwritepoint(depname, s.lat, s.lon);

return; 
