% process_autosal
%
% process salts from Brown
%
% must be called while matlab is in the directory where the raw autosal
% data is located.  This is typically in the CTD data subdirectory called
% salts.
%
%%%%%%
% USAGE NOTES:
% Update some fields in the *.m file
%   1. output_file: defaults to process_output.txt
%   2. standard_sea_water_batch: Make sure your standard seawater bathc number is listed below.
%      If not, add a new batch number where XXX is and a ration value where YYY appears and remove the comment
%   "%" symbol in the first column of the line.
%   3. i_want_many_plots =1 for plots, =0 for less plots
%   4. correction_method = 5 = (default)
%
%%%%%
%+
% Beginning of update section
%
output_file = 'process_output';

standard_sea_water_batch(130) = 0.99997;  %
standard_sea_water_batch(131) = 0.99986;  %
standard_sea_water_batch(132) = 0.99993;  %
standard_sea_water_batch(133) = 0.99986;  %
standard_sea_water_batch(134) = 0.99989;  %
standard_sea_water_batch(135) = 0.99992;  %
standard_sea_water_batch(136) = 0.99996;  %
standard_sea_water_batch(137) = 0.99995;  %
standard_sea_water_batch(138) = 0.99994;  %
standard_sea_water_batch(139) = 0.99993;  %
standard_sea_water_batch(140) = 0.99991;  %
standard_sea_water_batch(141) = 0.99993;  %
standard_sea_water_batch(142) = 0.99991;  %
standard_sea_water_batch(143) = 0.99989;  %
standard_sea_water_batch(144) = 0.99987;  %

standard_sea_water_batch(147) = 0.99987;  % 34.993

standard_sea_water_batch(149) = 0.99984;  % 34.994
standard_sea_water_batch(150) = 0.99978;  %
standard_sea_water_batch(151) = 0.99997;  % 34.999
standard_sea_water_batch(152) = 0.99981;  %
standard_sea_water_batch(153) = 0.99979;  %
standard_sea_water_batch(154) = 0.99990;  % 34.996
standard_sea_water_batch(155) = 0.99981;  % 34.993
%standard_sea_water_batch(XXX) = YYY; 
%
% change the plots if you want
i_want_many_plots = 0;  % yes=1 (plot out each reading); no=0 (only plot summary and duplicates)
%
% End update section
%-
%%%%%%


%%%%%
%+
% Open existing salts data base if one exists in this directory
%
clear autosal_salts
if exist ('autosal_salts_db.mat','file')
    load autosal_salts_db.mat
end
%-
%%%%%


%%%%%
%+
% Get a list of all *.dat files
% These are assumed to all be autosal output files.
% Files can have *any* filename structure.
% Assumes that all *.dat files have the corresponding *.hrd and *.raw files
% with the same root filename.
%
flist = dir('*.dat');

if length(flist) < 1
    sprintf('Warning: no *.dat files in this directory \n \t %s \nExiting now. \n',pwd)
    return
end
%-
%%%%%


%%%%%
%+
%
% load and process only new *.dat files
%
% First find out which files do not appear in autosal_salts_db.mat already
%   This is kept as an index to the filename that is now (as true=1)
%   When index_new_files(i) = 0, that file is read in.
%
%
index_new_files = zeros(length(flist),1);
if exist('autosal_salts','var')
    for II = 1: length(flist)
    [pathstr,this_name,ext] = fileparts(flist(II).name);
        for JJ = 1: length(autosal_salts)
            index_new_files(II) = max(strcmpi(this_name,autosal_salts(JJ).file),index_new_files(II));
        end
    end
else
    index_new_files = zeros(length(flist),1);
end
%-
%%%%%


%%%%%
%+
% Now read the new files below
%
% The routine read_autosal_dat_raw_mb (called below)
% creates the autosal_salts structure that looks like:
% tmpsalts (later to become 'salts')
%    file - one file name with *.dat file containing:
%         - each appears once for # of total reads including standards
%    sample_num         - linear number of sample from one run (or
%                           1000/1001 for standard in/out)
%    sample_id_str      - generated by oceanscientific as a combination of
%                           file name and sample_num.
%    bottle_label_str   - input by operator.  If in the following form: "cast# niskin#"
%                           then program will populate those fields,
%                           otherwise it will take a random guess at them.
%                           If only one field is present, assumes this is a
%                           niskin number and guesses the station number.  
%    sample_type        - =0 for standards, =1 for bottle reads
%    reading_num        - =0 for average of =1 one, =2 two, =3 three reads.
%    bath_t_str         - genrated by OceanScientific = bath temperature
%    uncorr_ratio       - uncorrect ratio
%    uncorr_ratio_stnd_dev - uncorrected ratio standard deviation (of all
%                           ten values for each reads)
%    correction         - 
%    adj_ratio          - 
%    calc_s             - for # of total reads including standards
%    calc_s_stnddev     - for # of total reads including standards
%    date_time          - for # of total reads including standards
%    dat_qc             - always 2.
% and these elements:  one for each of 10 values per single "read"
%    raw_sample_id_str
%    raw_reading_num
%    raw_uncorr_ratio
%    raw_adj_ratio
%    raw_calc_s
%    raw_qc
%    unique_sample_str
% All of these values are what was in the original file computed by the
% autosal software.
%
clear tmpsalts
KK_index_of_new_files = find(index_new_files == 0);
num_new_files = length(KK_index_of_new_files);
if (~isempty(KK_index_of_new_files))
    for II = 1: length(KK_index_of_new_files)
        flist(KK_index_of_new_files(II)).name
        tmpsalts(II) = read_autosal_dat_raw_mb (flist(KK_index_of_new_files(II)).name);
    end
else
    sprintf('No new *.dat files not already in autosal_salts_db.mat in this directory \n \t %s \nExiting now. \n',pwd)
    return
end
%-
%%%%%%


%%%%%%
%+
% Do some plotting for new files only
%
if num_new_files > 0
    if i_want_many_plots
        for II = 1: num_new_files
            plot_one_salt_run_raw_data(tmpsalts(II)); 
        end
    end
end
%-
%%%%%

%%%%%%
%+
%   make_one_btl_value_no_txt2
%   Adds the following fields to theResult
%           txt_station_id: [48x1 double]
%        txt_niskin_bottle: [48x1 double]
%             txt_samp_nbr: [48x1 double]
%            txt_tank_temp: [48x1 double]
%           txt_cond_ratio: [48x1 double]
%                 txt_date: [48x1 double]
%                   txt_qc;
%
if num_new_files > 0
    clear tmp
    for II = 1: num_new_files    
        tmp(II) = make_one_btl_value_no_txt3 (tmpsalts(II));  % adds the averages to the file as *txt* elements
    end
    tmpsalts = tmp;
end
%-
%%%%%%

%%%%%%
%+
% correct_autosal_drift in matlab structure
% (don't read in and write out new files)
% This uses the txt values already present in autosal_salts structure
%
if num_new_files > 0
    clear tmp
    for II = 1: num_new_files
        % 2.*standard_sea_water_batch(tmpsalts(II).batch_number)
        tmp(II) = correct_autosal_drift2_matlab_structure_no_file3 (2.*standard_sea_water_batch(tmpsalts(II).batch_number), tmpsalts(II), 7, 0, 0);
        close all;
    end
    tmpsalts = tmp;
    clear tmp;
end
%-
%%%%%%


%%%%%
%+
% Now add the new stations to the existing autosal_salts structure
%
if (exist('autosal_salts','var'))
    II_index_existing_salt_files = length(autosal_salts);
    % now append the new salt values to the end of the existing ones
    autosal_salts(end+1:end+num_new_files) = tmpsalts;
else
    autosal_salts = tmpsalts;
end
%-
%%%%%

%%%%%
%+
% Create new output file
% With the corrected newtxt values create a salt file
%
output = make_salinity_file5 (autosal_salts);
[output, duplicate_S] = average_and_remove_duplicates_aoml4 (output, [output_file '_salts.txt']);

%
% Get all calibration data
%
calibration_data = get_calibrations_from_salts2 (autosal_salts);

save autosal_salts_db.mat autosal_salts
save calibations_db.mat calibration_data

%
% Plot calibration data time series
%
min_date = min(calibration_data(1).txt_date);
max_date = max(calibration_data(end).txt_date);
figure;
for II = 1: length(calibration_data)
    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_date(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'og', 'Markerfacecolor', 'g','MarkerSize',10)
    min_date = min(min(calibration_data(II).txt_date( I_start)), min_date);
    max_date = max(max(calibration_data(II).txt_date( I_start)), max_date);
    hold on;

    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==4);
    plot (calibration_data(II).txt_date(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'ok', 'Markerfacecolor', 'k','MarkerSize',10);
    if ~isempty(I_start)
        min_date = min(min(calibration_data(II).txt_date( I_start)), min_date);
        max_date = max(max(calibration_data(II).txt_date( I_start)), max_date);
    end
    
    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_date( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vr', 'Markerfacecolor', 'r','MarkerSize',10)
    min_date = min(min(calibration_data(II).txt_date( I_start)), min_date);
    max_date = max(max(calibration_data(II).txt_date( I_start)), max_date);

    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==4);
    plot (calibration_data(II).txt_date( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vy', 'Markerfacecolor', 'y','MarkerSize',10)
    if ~isempty(I_start)
        min_date = min(min(calibration_data(II).txt_date( I_start)), min_date);
        max_date = max(max(calibration_data(II).txt_date( I_start)), max_date);
    end
    %
    % plot these good values again so that their symbols appear on the top
    % of the plot
    %
    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_date(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'og', 'Markerfacecolor', 'g','MarkerSize',10)
    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_date( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vr', 'Markerfacecolor', 'r','MarkerSize',10)
end
datetick
set(gca,'xlim',[min_date-0.1*(max_date-min_date) max_date+0.1*(max_date-min_date)])
legend('Good Start','Bad Start','Good End','Bad End','Location','best')
set(gca,'Box','off');
axesposition=get(gca,'Position');
H=gca;
ylimits = get(gca,'ylim');
hNewAxes = axes('Position', axesposition, 'color','none','Ylim',[sw_sals(ylimits(1)/100,24) sw_sals(ylimits(2)/100,24)],'Yaxislocation','right','xtick',[],'box','off');
axes(H)
grid
print -depsc plot/calibrations_vs_date.ps


%%%%%
%+
% Get rid of weird large station numbers for axis limits
%
temp = vertcat(calibration_data.txt_station_id);
temp(find(temp>=900)) = [];
min_station = min(temp);
max_station = max(temp);
%
figure;
for II = 1: length(calibration_data)
    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_station_id(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'og', 'Markerfacecolor', 'g','MarkerSize',10)
    hold on;
    min_station = min(min(calibration_data(II).txt_station_id(I_start)), min_station);
    max_station = max(max(calibration_data(II).txt_station_id(I_start)), max_station);

    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==4);
    plot (calibration_data(II).txt_station_id(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'ok', 'Markerfacecolor', 'k','MarkerSize',10);
    if ~isempty(I_start)
        temp = vertcat(calibration_data.txt_station_id);
        temp(find(temp==999)) = [];
        min_station = min(temp);
        max_station = max(temp);
        %min_station = min(min(calibration_data(II).txt_station_id(I_start)), min_station);
        %max_station = max(max(calibration_data(II).txt_station_id(I_start)), max_station);
    end
    
    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_station_id( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vr', 'Markerfacecolor', 'r','MarkerSize',10)
    min_station = min(min(calibration_data(II).txt_station_id(I_start)), min_station);
    max_station = max(max(calibration_data(II).txt_station_id(I_start)), max_station);

    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==4);
    plot (calibration_data(II).txt_station_id( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vy', 'Markerfacecolor', 'y','MarkerSize',10)
    if ~isempty(I_start)
        temp = vertcat(calibration_data.txt_station_id);
        temp(find(temp==999)) = [];
        min_station = min(temp);
        max_station = max(temp);
        %min_station = min(min(calibration_data(II).txt_station_id(I_start)), min_station);
        %max_station = max(max(calibration_data(II).txt_station_id(I_start)), max_station);
    end
    
    %
    % plot these dgood values again so that their symbols appear on the top
    % of the plot
    %
    I_start = find(calibration_data(II).txt_samp_nbr==1001 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_station_id( I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'vr', 'Markerfacecolor', 'r','MarkerSize',10)
    I_start = find(calibration_data(II).txt_samp_nbr==1000 & calibration_data(II).txt_qc==2);
    plot (calibration_data(II).txt_station_id(I_start), 100/2*calibration_data(II).txt_cond_ratio(I_start), 'og', 'Markerfacecolor', 'g','MarkerSize',10)
    set(gca,'xlim',[min_station max_station])
end
xlabel('Station Number')
set(gca,'xlim',[min_station-0.1*(max_station-min_station) max_station+0.1*(max_station-min_station)])
legend('Good Start','Bad Start','Good End','Bad End','Location','best')
set(gca,'Box','off');
H=gca;
axesposition=get(gca,'Position');
ylimits = get(gca,'ylim');
hNewAxes = axes('Position', axesposition, 'color','none','Ylim',[sw_sals(ylimits(1)/100/2,24) sw_sals(ylimits(2)/100/2,24)],'Yaxislocation','right','xtick',[],'box','off');
axes(H)
set(H,'Box','on')
grid
title('Standard Water Values (x100) vs Station (left psu)') 
print -depsc plot/calibrations_vs_station_num.ps


