% Environment variables for oxygen calibration routines
dbstop if error
p_calib=1;  % run calibration (run_oxygen_cal_ml.m) (after Steps 1-2)
p_plot=0;   % Intermediate plots: dv/dt and ctd o2 vs bottles
p_final=1;  % to printout figures
p_filter=9; % filter oxygen voltage over xx points (for match_sgn and calctd)
p_keep_all_btl=1; %keep info for all bottles even if no winkler (AG, Aug 2015)

cruiseName = 'MOORSPICE';
cname=lower(cruiseName);
%driveName = '/Users/marionsofiaalberty/Dropbox/MoorSPICE-share/';
 driveName = '/Users/ganacho/Dropbox/MoorSPICE-share/';

addpath([driveName 'Routines/Alex/']) ; 
addpath([driveName 'Routines/Robbins/']);
addpath([driveName 'MOORSPICE/Oxygen/Calibration/TMP']);

%%%%%%%%%%%%%%%%%%%%%%
% INIT PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%

% Location of CTD data
%----------------------
rootDir = 'data-processing';
% bottle directory
%-----------------
dirbtl = strcat(driveName,cruiseName,filesep,rootDir,...
    filesep,'CTD',filesep);  % dir of _btl.nc

dirctd = dirbtl; % dir of _ctd.nc

%
% working directory
%
dirwork=driveName; %working directory
subdirwork='Calibration';

% NEED TO CREATE SUBDIRECTORY 'Califiles'

%oxygen bottle (.xls) files (xlsx with all stations)
diro2=[dirwork cruiseName '/Oxygen/Chemistry/'];
%oxygen bottle (.xls) files (xlsx with all stations)
excelExtension = 'xlsx'; %system-dependent (linux/mac/windows)
%excelExtension = 'xls';
nbtl=24; %max number of bottle to read from chemistry

%
% output directories (.cal)
%
dirmerge = strcat(dirwork, cruiseName, filesep, 'Oxygen',...
    filesep,subdirwork,filesep,'Califiles', filesep);

%
% graphics
%
dirjpg=[dirmerge 'PNG' filesep];
direps=[dirmerge 'EPS' filesep];

% LIST OF PROFILES TO USE TO DETERMINE CALIBRATION COEFFICIENTS
% (THOSE WHICH HAVE OXYGEN DATA)
% gprof2get=100*(1:length(gicast2get))'+gicast2get';
% and add station 0
gprof2get=[
101
201
301
401
501
601
701
801
901
1001
1101
1201
1301
1401
1501
1601
1701
1801
1901
2001
2101
2201
2301
2401
2501
2601
2701
2801
2901
3001
3101
3201
3301
3401
3501
3601
3701
3801
3901
4001
4101
4201
4202
4301
4501
4601
4701
4801
4901
5001
5006
5101
5201
5301
5401
5501
5601
5701
];

if 0
    %INDEX TO SUBTRACT TO STATION NUMBER TO GET THE TAB IN EXCEL SHEETS
    %by default there is one file every 10 station so that the tab index is
    gistat=floor(gprof2get/100);
    gishiftsheet=10*floor((gistat-1)/10);
    
    % sheet index corresponding to the station within Excel file
    % but special one for stat 60 which had 2 casts
    % gishiftsheet(gprof2get==6013)=49;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SUBSELECT BY STATIONS IF NECESSARY FOR MERGE_CTDBTL_CHIMIE
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    gstat=1:57; %group of stations to treat
    %SELECT PROFILES FOR THESE STATIONS
    [gprof,gii]=find_profile_fromstat(gprof2get,gstat);
    
    gishiftsheet=gishiftsheet(gii);
else
    p_read_xls_all=1; %read all xls chemistry information at once (eg for Bifurcation)
    p_read_xls_all=0; %read xls chemistry with one file every 10 station (eg for MoorSPICE)
    gprof=gprof2get;
    %gprof=5701; disp('ONLY 5701 !!!'); pause
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ELIMINATE SOME PROFILES IF NECESSARY
%   remove some stations 7/11 problems with sampling; no 20,30,35
% for gprofbad=[701,1101,2001,3001,3501,8101] 
%     gishiftsheet(find(gprof==gprofbad))=[];
%     gprof(find(gprof==gprofbad))=[];
% end

% sensors to use (necessarily for the whole group for O2 and for outputs)
gisensorO=1; 
gisensorS=1;
gisensorT=1;

% OPTION TO LIST CHANNEL TO USE IN THE .BTL FILE (UPCAST) FOR EACH STATION 
% FOR MERGE CTD BOTTLE ONLY
gisensorSu=ones(size(gprof));
gisensorTu=gisensorSu;

if exist('p_merge')&p_merge
    return %stops here for parameters necessary to merge_ctdbtl_chimie
end

%%%%%%%%%%%%%%%%%%%%%
% STEP 2 : MATCH_SGN
%%%%%%%%%%%%%%%%%%%%%%

%Equation to use for initial oxygen calculation (see oxyg_eq.m)
p_equation=1; %PMEL
% p_equation=2; %SEABIRD 2 (ABC equation)
% p_equation=3; %SEABIRD 3 (pcor equation)

if gisensorO==1 %SN 0023 : primary circuit *******
    switch p_equation
        % see routine
        % Below is calibration data from O2 sensor 0023, pre-cruise
        % calibration on 09 November 2013
        case {1} % PML Equation (from a mix of the two calibrations below)
            soc=0.4646;         % from eq. 3
            voffset=-0.5200;    % from eq. 3
            E=3.6000e-002;      % from eq. 2 (not exactly pcor)
            tau_20=1.3600;      % from eq. 2 (response time at surface/20dC (seconds))
            tcor=-0.0004;       % from eq.3
            %fixed coefficients
            sbecoefd.d1=1.92634e-004;   % response time vs pressure
            sbecoefd.d2=-4.64803e-002;  % response time vs temperature
            sbecoef0(1)=soc;
            sbecoef0(2)=voffset;
            sbecoef0(3)=tau_20;
            sbecoef0(4)=tcor;
            sbecoef0(5)=E;
        case {2} % ABC equation
            %Soc*(V+Voffset)*(1.0+A*T+B*T2 +C*T3)*OxSol(T,S)*exp(E*P/K)
            soc=0.4646;
            voffset=-0.5200;
            E=3.6000e-002;      % or E for that equation
            tau_20=1.3600;
            %fixed coefficients
            sbecoefd.A=-4.0206e-003;
            sbecoefd.B=1.9156e-004;
            sbecoefd.C=-2.8340e-006;
            sbecoefd.d0=2.5826e000;
            sbecoefd.d1=1.92634e-004;   % response time vs pressure
            sbecoefd.d2=-4.64803e-002;  % response time vs temperature
            sbecoef0(1)=soc;
            sbecoef0(2)=voffset;
            sbecoef0(3)=NaN;
            sbecoef0(4)=NaN;
            sbecoef0(5)=E;
        case 3 % pcor equation
            %(Soc * (V + Voffset)) * exp(Tcor * T) * Oxsat(T,S) * exp(Pcor * P)
            soc=0.4646;
            boc=0.0000;
            voffset=-0.5200;
            tcor=-0.0004;
            pcor=1.350e-04;
            sbecoefd.d1=1.92634e-004;  %response time vs pressure
            sbecoefd.d2=-4.64803e-002; %response time vs temperature
            sbecoef0(1)=soc ;
            sbecoef0(2)=voffset ;
            sbecoef0(3)=tcor ;
            sbecoef0(4)=pcor ;
            sbecoef0(5)=NaN ;
        otherwise
            error('calibration coefficients not set')
    end
    %
    %
elseif gisensorO==2 % Sensor SN 1811 : secondary circuit *******
    switch p_equation
        % see routine
        % Below is calibration data from O2 sensor 1811, pre-cruise
        % calibration on 04 September 2013
        case 1 %From a mix of the two calibration below
            soc=0.4761; 
            voffset=-0.4906;
            E=3.6000e-002;
            tau_20=1.4500; 
            tcor=-0.0019;
            %fixed coefficients
            sbecoefd.d1=1.92634e-004;  % response time vs pressure
            sbecoefd.d2=-4.64803e-002; % response time vs temperature
            sbecoef0(1)=soc ;
            sbecoef0(2)=voffset ;
            sbecoef0(3)=tau_20 ;
            sbecoef0(4)=tcor ;
            sbecoef0(5)=E ;
        case 2 %Pre-cruise calibration of 14-Mar-2012; ABC equation
            %Soc*(V+Voffset)*(1.0+A*T+B*T2 +C*T3)*OxSol(T,S)*exp(E*P/K)
            soc=0.4761; 
            voffset=-0.4906;
            E=3.6000e-002;      % not exactly pcor
            tau_20=1.4500;
            %fixed coefficients
            sbecoefd.A=-3.7641e-003;
            sbecoefd.B=1.9546e-004;
            sbecoefd.C=-2.9427e-006;
            sbecoefd.d0=2.5826e000;
            sbecoefd.d1=1.92634e-004;   % response time vs pressure
            sbecoefd.d2=-4.64803e-002;  % response time vs temperature
            sbecoef0(1)=soc ;
            sbecoef0(2)=voffset ;
            sbecoef0(3)=NaN ;
            sbecoef0(4)=NaN ;
            sbecoef0(5)=E ;
        case 3 %Pre-cruise calibration of; pcor equation
            soc=0.4761; 
            voffset=-0.4906;
            %%%%%% BELOW COEFFICIENTS ARE MISSING!! %%%%%%%%%%%%%%%%%%%%%%%
            boc=NaN;
            tcor=NaN;
            pcor=NaN;
            %%%%%% ABOVE COEFFICIENTS ARE MISSING!! %%%%%%%%%%%%%%%%%%%%%%%
            sbecoefd.d1=1.92634e-004;%response time vs pressure
            sbecoefd.d2=-4.64803e-002;%response time vs temperature
            sbecoef0(1)=soc ;
            sbecoef0(2)=voffset ;
            sbecoef0(3)=tcor ;
            sbecoef0(4)=pcor ;
            sbecoef0(5)=NaN ;
        otherwise
            error('calibration coefficients not set')
    end
end

p_plot_upcastbottles=1; %plot upcast bottle o2 on initial residual graphs
%SET A MAXIMUM DEPTH UNDER WHICH WE DO NOT MATCH SIGMAS
% THIS IS MOSTLY FOR STATION 60 WHERE THERE IS NOT T/S GRADIENT BELOW
% 3400M THAT MAKES THE SIGMA MATCH UNSTABLE
p_gidmax_match=7000*ones(length(gprof),1);
p_gidmax_match((floor(gprof/100)==42))=3400; %station 60 because homogeneous 
      % densities below 3400 that generates large erroneous density adjustments


if exist('p_match') & p_match
    p_match=0;
    return %stops here for parameters
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% STEP 3 : run_oxygen_cal_ml -- CALIBRATION: set p_calib=1 above and run this file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibration : optimization SBE coefficients 
%   based on O2ctd-O2sample @ bottle depths
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% HERE WE USE PMEL EQUATION TO MINIMIZE THE RESIDUALS
%
% HAVE TO ASK SEABIRD FOR THE SPECIFIC CALIBRATION COEFFICIENTS WITH TCORR
%
% C. JOHNSON'S EMAIL OF APRIL 2013:
% Yes, we still model the SBE-43 oxygen sensor calibration using that equation.
% We tried the SBE Hysteresis correction to remove the up-down profile
% differences, but found that in the end we got better results by matching
% the up cast water samples to the down cast oxygen voltages using
% apppropriately referenced potential densities. and employing the equation
% below to minimize the difference between those quantities.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET HERE CONFIGURATION FOR CALIBRATION
% SET GROUP OF STATIONS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear csuf
% % Experiment 1 - Global residual is 0.004 std dev is 1.017; rms is 1.016
% gistat_cal=1:57; csuf='1_57_ex1'; iswitch = 1; wswitch = 0; p_maxd = 1e4;

% % Experiment 2 - Global residual is -0.003 std dev is 1.334; rms is 1.333
% gistat_cal=1:57; csuf='1_57_ex2'; iswitch = 0; wswitch = 0; p_maxd = 3510;

% % Experiment 3 - Global residual is 0.004 std dev is 1.081; rms is 1.059
% gistat_cal=1:3; csuf='1_3_ex3'; iswitch = 0; wswitch = 0; p_maxd = 3510;

% % Experiment 4 - Global residual is -0.003 std dev is 0.956; rms is 0.955
% gistat_cal=4:43; csuf='4_43_ex4'; iswitch = 0; wswitch = 0; p_maxd = 3510;

% % Experiment 5 - Global residual is 0.006 std dev is 1.098; rms is 1.095
% gistat_cal=45:57; csuf='45_57_ex5'; iswitch = 0; wswitch = 0; p_maxd = 3510;

% % Experiment 6 - Global residual is -0.003 std dev is 0.763; rms is 0.762
%  gistat_cal=4:43; csuf='4_43_ex6'; iswitch = 1; wswitch = 0; p_maxd = 3510;

% % Experiment 7, used ABC equation in Step 2
% gistat_cal=4:43; csuf='4_43_ex7'; iswitch = 1; wswitch = 0; p_maxd = 3510;

% % Experiment 8 - Global residual is -0.003 std dev is 0.759; rms is 0.758
% gistat_cal=4:43; csuf='4_43_ex8'; iswitch = 1; wswitch = 0; p_maxd = 3010;

% % Experiment 9 - Global residual is -0.163 std dev is 0.843; rms is 0.857
% gistat_cal=4:43; csuf='4_43_ex9'; iswitch = 1; wswitch = 1; p_maxd = 3010;

% % Experiment 10 - Global residual is -0.370 std dev is 1.490; rms is 1.533
% gistat_cal=4:43; csuf='4_43_ex10'; iswitch = 0; wswitch = 1; p_maxd = 3010;

% % Experiment 11 - Global residual is -0.041 std dev is 0.772; rms is 0.772
% gistat_cal=4:43; csuf='4_43_ex11'; iswitch = 1; wswitch = 2; p_maxd = 3010;

% % Experiment 12 - Global residual is 0.068 std dev is 0.763; rms is 0.764
% gistat_cal=4:43; csuf='4_43_ex12'; iswitch = 1; wswitch = 3; p_maxd = 3010;

% % Experiment 13, this configuration fails for iswitch = 2
% gistat_cal=4:43; csuf='4_43_ex13'; iswitch = 2; wswitch = 0; p_maxd = 3010;

% % Experiment 14 - Global residual is -0.100 std dev is 0.897; rms is 0.901
% gistat_cal=1:57; csuf='1_57_ex14'; iswitch = 1; wswitch = 2; p_maxd = 2000;

% % Experiment 15 - Global residual is -0.089 std dev is 0.892; rms is 0.895
% gistat_cal=1:57; csuf='1_57_ex15'; iswitch = 1; wswitch = 2; p_maxd = 2500;

% % Experiment 16 - Global residual is -0.082 std dev is 0.899; rms is 0.902
% gistat_cal=1:57; csuf='1_57_ex16'; iswitch = 1; wswitch = 2; p_maxd = 3000;

% % Experiment 17 - Global residual is 0.000 std dev is 0.940; rms is 0.939
% gistat_cal=1:57; csuf='1_57_ex17'; iswitch = 1; wswitch = 0; p_maxd = 3510;

% % Experiment 18 - Global residual is 0.004 std dev is 0.902; rms is 0.899
% gistat_cal=45:57; csuf='45_57_ex18'; iswitch = 1; wswitch = 0; p_maxd = 3510;


% Residuals ranges for plotting
gdeprange=[0 500 1500 2500 3510]; 

if p_calib
    %USE INITIAL COEFFICIENT FROM NEW PRE-CALIBRATION SHEET THAT GO WITH
    %THE TCOR EQUATION (parameters are set in this file, above)
    
    [calmatrixv, btl0]=run_oxygen_cal_ml(sbecoef0,sbecoefd,gistat_cal,iswitch,wswitch,...
        dirjpg,direps,dirmerge,cname,p_final,csuf,p_maxd,gdeprange,p_equation,p_plot);
    
    final=calmatrixv;
    
    eval(sprintf('save %sox_cals_%s_is%d_ws%d_%s.mat final sbecoefd btl0',dirmerge,cname,iswitch,wswitch,csuf))
    
    %Save samples and outliers at bottle levels
    %First file with pressures sigma-matched to downcast
    fid=fopen(sprintf('%ssample_with_qc_%s_is%d_ws%d_%s_dpres.txt',dirmerge,cname,iswitch,wswitch,csuf),'w');
    zoutliers=[btl0.stat,btl0.cast,btl0.btlid, btl0.dpres,btl0.dtemp68,btl0.dsali,btl0.wnklo2um,btl0.wnklo2_flag,...
        btl0.bsali, btl0.oxcal];
    fprintf(fid,'stat  cast bottleID dpres   utemp68  uctdsali o2sample (umol) o2flag (3=questionnable) salisample dctdo2\n');
    fprintf(fid,'%4d %4d %4d %10.1f %10.4f %10.4f %10.4f %4d %10.4f %10.4f\n',zoutliers');
    fclose(fid);

    %Second file with pressures ab tottle level (upres)
    fid=fopen(sprintf('%ssample_with_qc_%s_is%d_ws%d_%s_upres.txt',dirmerge,cname,iswitch,wswitch,csuf),'w');
    zoutliers=[btl0.stat,btl0.cast,btl0.btlid, btl0.upres,btl0.dtemp68,btl0.dsali,btl0.wnklo2um,btl0.wnklo2_flag,...
        btl0.bsali, btl0.oxcal];
    fprintf(fid,'stat  cast bottleID upres   temp68  uctdsali o2sample (umol) o2flag (3=questionnable) salisample dctdo2\n');
    fprintf(fid,'%4d %4d %4d %10.1f %10.4f %10.4f %10.4f %4d %10.4f %10.4f\n',zoutliers');
    fclose(fid);
%
% End calibration parameters/run
%
    return
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 4 : calctd parameters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calibration of the entire CTD profiles from new coefficients
% calctd (with cleanox inside)
%
% !!!!!!!!!!!
% DELETE ALL CALIBRATION FILES THAT ARE UNNECESSARY TO AVOID CONFUSION
% (Califiles/ox_cals_***.mat)
%
clear iswitch wswitch %to avoid confusion: re-set them below along with stations /casts to treat

p_saveoutput=1; % generally 1 except if we do further exploratory experiment after a "final" dataset
p_onlyonecast=0;   % calibrate all casts within range if 1
% statgroup='1_3';
% statgroup='4_13';
% statgroup='14_30';
% statgroup='31_41';
% statgroup='42_44';
statgroup='45_57';


if p_onlyonecast  % new routine to calibrate all profiles
    p_cleanox=1; % do not use cleanox if only one cast (useless)
    p_cleansalt=0; % neither cleansalt (comes afterwards)
    p_cleantemp=0;
    %gistat_nosample{0+1}=1; %CAREFUL: THE INDEX IS STAT+1 TO HAVE STATION ZERO (Use station 1 to get coefficients for station 0)
    
    % Stations we will calibrate
    gistat2cal=1:3; iswitch=0; wswitch=0;
%     gistat2cal=4:43; iswitch=1; wswitch=0;
%     gistat2cal=[45:57]; iswitch=0; wswitch=0;
    
    % Find correct file name for corresponding experiment
    csuf='1_3_ex3';
%     csuf='4_43_ex6';
%     csuf='45_57_ex5';
    
    %NOTE THAT CALIBRATION COEFFICIENTS ARE TAKEN STATION BY STATION (NO
    statgroup='1_3';
%     statgroup='4_43' ;
%     statgroup='45_57' ;
    
    %gprof2cal=1; %for station zero
    
else %all profiles/casts
    % NEEDS WORK
    p_cleanox=0; %do not use cleanox the first time as we verify the theta-S relations
    p_cleansalt=0; %clean salinity after finished cleaning oxygen (uses the o2 ctd output files/station)
    p_cleantemp=1; %clean temperature: after p_cleansalt iteration
    switch statgroup
        case '1_3'
            statbeg=1; castbeg=1;
            statend=3; castend=1;
            iswitch=0; %slope
            wswitch=0; %weight on deep measurements
            csuf='1_3_ex3';
        case '4_13'
            statbeg=4; castbeg=1;
            statend=13; castend=1;
            iswitch=1; %slope
            wswitch=0; %weight on deep measurements
            csuf='4_43_ex6';
        case '14_30'
            statbeg=14; castbeg=1;
            statend=30; castend=1;
            iswitch=1; %slope
            wswitch=0; %weight on deep measurements
            csuf='4_43_ex6';
        case '31_41'
            statbeg=31; castbeg=1;
            statend=41; castend=2;
            iswitch=1; %slope
            wswitch=0; %weight on deep measurements
            csuf='4_43_ex6';
        case '42_44'
            statbeg=42; castbeg=1;
            statend=44; castend=1;
            iswitch=1; %slope
            wswitch=0; %weight on deep measurements
            csuf='4_43_ex6';
            gistat_nosample{44+1}=43; %CAREFUL: THE INDEX IS STAT+1 TO HAVE STATION ZERO
        case '45_57'
            statbeg=45; castbeg=1;
            statend=57; castend=4;
            iswitch=0; %slope
            wswitch=0; %weight on deep measurements
            csuf='45_57_ex5';
    end
end
%   p_auto is for automated interpolation for calctd/cleanox (0 or 1 for each station ; use if many small gaps)
p_auto=0*ones(200,1); %automated interpolation (index=profile)
%set from which calibrated station do we get the coefficients when have not
%sample/calibration
%gistat_nosample{0}=1; %impossible...

p_debug=0; %this plots ctd.doxv, sdvdt, dvdt and filtered for verification during calibration step


