From b95d41debc06f51da56a2482ab52e979c7e4aa82 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Sat, 27 Jan 2024 09:38:56 -0800
Subject: [PATCH 1/8] Fix comment

---
 Code/diss2combo.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Code/diss2combo.m b/Code/diss2combo.m
index 963a57a..3925e8f 100644
--- a/Code/diss2combo.m
+++ b/Code/diss2combo.m
@@ -12,4 +12,4 @@ function diss2combo(binned, pars)
 
 [a, fnCombo] = save2combo(binned, pars, pars.diss_combo_root);
 save2NetCDF(a, fnCombo, pars);
-end % profiles2combo
\ No newline at end of file
+end % diss2combo
\ No newline at end of file

From 43834cbb801122bf714d8a64bb31ae1b63710604 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Sat, 27 Jan 2024 09:39:21 -0800
Subject: [PATCH 2/8] Add in hooks for chi processing

---
 Code/chi2binned.m      | 104 +++++++++++++++++++++++++++++++++++++++++
 Code/chi2combo.m       |  15 ++++++
 Code/get_info.m        |   6 +++
 Code/process_P_files.m |  61 +++++++++++++++---------
 Code/profile2chi.m     |  73 +++++++++++++++++++++++++++++
 Code/update_paths.m    |  12 +++++
 6 files changed, 250 insertions(+), 21 deletions(-)
 create mode 100644 Code/chi2binned.m
 create mode 100644 Code/chi2combo.m
 create mode 100644 Code/profile2chi.m

diff --git a/Code/chi2binned.m b/Code/chi2binned.m
new file mode 100644
index 0000000..dd753dd
--- /dev/null
+++ b/Code/chi2binned.m
@@ -0,0 +1,104 @@
+% Bin chi into depth bins
+%
+% July-2023, Pat Welch, pat@mousebrains.com
+
+function [row, retval] = chi2binned(row, a, pars)
+arguments (Input)
+    row (1,:) table % row to work on
+    a struct % Output of profile2chi, struct or empty
+    pars struct % Parameters, defaults from get_info
+end % arguments Input
+arguments (Output)
+    row table % row worked on
+    retval (2,1) cell % {filename or missing, empty or binned profile}
+end % arguments Output
+
+retval = {missing, []};
+
+if ~row.qProfileOkay
+    return;
+end % if ~row.qProfileOkay
+
+warning("Chi binning not implemented");
+
+% fnChi = row.fnChi;
+% fnBin = fullfile(pars.chi_binned_root, append(row.name, ".mat"));
+% row.fnChiBin = fnBin;
+% 
+% if isnewer(fnBin, fnChi)
+%     retval = {fnBin, []}; % We want this for combining
+%     fprintf("%s: %s is newer than %s\n", row.name, row.fnBin, row.fnChi);
+%     return;
+% end % if isnewer
+% 
+% if isempty(a)
+%     fprintf("Loading %s\n", row.fnChi);
+%     a = load(row.fnChi);
+% end % if isempty
+% 
+% %% Bin the data into depth bins
+% 
+% pInfo = a.info;
+% profiles = a.profiles;
+% 
+% fprintf("%s: Binning %d profiles\n", row.name, numel(profiles));
+% 
+% if pars.profile_direction == "time" % Bin in time
+%     binSize = seconds(pars.binChi_width); % Bin stepsize in (sec)
+%     keyName = "t";
+%     binFunc = @bin_by_time;
+%     glueFunc = @glue_lengthwise;
+% else % Bin by depth
+%     binSize = pars.binChi_width; % Bin stepsize (m)
+%     keyName = pars.binChi_variable; % Variable to bin by
+%     binFunc = @bin_by_real;
+%     glueFunc = @glue_widthwise;
+% end % if profile_direction
+% 
+% casts = cell(numel(profiles),1);
+% for index = 1:numel(profiles)
+%     profile = profiles{index};
+%     nE = size(profile.e, 2);
+%     prof2 = table();
+%     for name = string(profile.Properties.VariableNames)
+%         sz = size(profile.(name),2);
+%         if ~ismatrix(profile.(name)), continue; end
+%         if sz == 1
+%             prof2.(name) = profile.(name);
+%         elseif sz == nE
+%             for j = 1:sz
+%                 prof2.(append(name, "_", string(j))) = profile.(name)(:,j);
+%             end % for j;
+%         end % if sz
+%     end % for name
+% 
+%     casts{index} = binFunc(binSize, keyName, prof2, pars.binChi_method);
+% end % for index
+% 
+% qDrop = cellfun(@isempty, casts); % This shouldn't happend
+% 
+% if any(qDrop)
+%     casts = casts(~qDrop);
+%     pInfo = pInfo(~qDrop,:);
+% end % any qDrop
+% 
+% if isempty(casts)
+%     row.qProfileOkay = false;
+%     fprintf("%s: No usable casts found in %s\n", row.name, row.fnProf);
+%     return;
+% end
+% 
+% tbl = glueFunc("bin", casts);
+% 
+% binned = struct ( ...
+%     "tbl", tbl, ...
+%     "info", pInfo);
+% if isfield(a, "fp07")
+%     binned.fp07 = a.fp07;
+% end % if isfield
+% 
+% my_mk_directory(fnBin);
+% save(fnBin, "-struct", "binned", pars.matlab_file_format);
+% fprintf("%s: Saving %d profiles to %s\n", row.name, size(binned.info,1), fnBin);
+% retval = {fnBin, binned};
+end % chi2binned
diff --git a/Code/chi2combo.m b/Code/chi2combo.m
new file mode 100644
index 0000000..80470b4
--- /dev/null
+++ b/Code/chi2combo.m
@@ -0,0 +1,15 @@
+% Combine the output of chi2binned into a single table, sorted by time
+%
+% This is a refactorization to deal with parallelization.
+%
+% Oct-2023, Pat Welch, pat@mousebrains.com
+
+function chi2combo(binned, pars)
+arguments (Input)
+    binned (:,1) cell
+    pars struct
+end % arguments Input
+
+[a, fnCombo] = save2combo(binned, pars, pars.chi_combo_root);
+save2NetCDF(a, fnCombo, pars);
+end % chi2combo
\ No newline at end of file
diff --git a/Code/get_info.m b/Code/get_info.m
index b283a7a..1e8ba49 100644
--- a/Code/get_info.m
+++ b/Code/get_info.m
@@ -107,6 +107,8 @@
 addParameter(p, "diss_T2_norm", 1, validPositive); % Value to multiple T2_fast temperature probe by to calculate mean for dissipation estimate
 addParameter(p, "diss_warning_fraction", 0.1); % When to warn about difference of e probes > diss_warning_ratio
 addParameter(p, "diss_epsilon_minimum", 1e-13, validPositive); % Dissipation estimates less than this are set to nan, for bad electronics
+%% Chi parameters
+addParameter(p, "chi_enable", false, validLogical); % Calculate chi
 %% Binning parameters for profiles, non-dissipation
 addParameter(p, "bin_method", "mean", validMethod); % Which method to use to combine bins together
 addParameter(p, "bin_width", 1, validPositive); % Bin width in (m)
@@ -115,6 +117,10 @@
 addParameter(p, "binDiss_method", "mean", validMethod); % Which method to use to combine bins together
 addParameter(p, "binDiss_width", 1, validPositive); % Bin width in (m)
 addParameter(p, "binDiss_variable", "depth", validString); % Which variable to bin by, unless direction==time
+%% Binning parameters for chi
+addParameter(p, "binChi_method", "mean", validMethod); % Which method to use to combine bins together
+addParameter(p, "binChi_width", 1, validPositive); % Bin width in (m)
+addParameter(p, "binChi_variable", "depth", validString); % Which variable to bin by, unless direction==time
 %% CTD time binning parameters
 addParameter(p, "ctd_bin_enable", true, validLogical); % Should CTD be binned outside of profiles?
 addParameter(p, "ctd_bin_dt", 0.5, validPositive); % Width in seconds of CTD binning
diff --git a/Code/process_P_files.m b/Code/process_P_files.m
index c04aa5a..c6c6c6b 100644
--- a/Code/process_P_files.m
+++ b/Code/process_P_files.m
@@ -64,9 +64,10 @@
 
     params = parallel.pool.Constant(pars); % Doesn't change from here on
 
-    binnedProfile = cell(size(p_filenames,1),1);
-    binnedDiss = cell(size(binnedProfile));
-    binnedCTD = cell(size(binnedProfile));
+    profileBinned = cell(size(p_filenames,1),1);
+    dissBinned = cell(size(profileBinned));
+    chiBinned = cell(size(profileBinned));
+    ctdBinned = cell(size(profileBinned));
 
     parfor index = 1:size(p_filenames,1)
         st = tic();
@@ -80,7 +81,12 @@
         try
             warning("off", "MATLAB:dispatcher:nameConflict");
             if row.qMatOkay
-                [row, binnedProfile{index}, binnedDiss{index}, binnedCTD{index}] = P_to_binned_profile(row, params.Value);
+                [row, ...
+                    profileBinned{index}, ...
+                    dissBinned{index}, ...
+                    chiBinned{index}, ...
+                    ctdBinned{index} ...
+                    ] = P_to_binned_profile(row, params.Value);
                 qMatOkay(index) = row.qMatOkay;
                 qProfileOkay(index) = row.qProfileOkay;
             else
@@ -109,22 +115,27 @@
     my_mk_directory(qUseDB, pars.debug);
     save(qUseDB, "qUse", pars.matlab_file_format);
 
-    binnedProfile = binnedProfile(~cellfun(@isempty, binnedProfile)); % Prune empty bins
-    binnedDiss = binnedDiss(~cellfun(@isempty, binnedDiss)); % Prune empty bins
-    binnedCTD = binnedCTD(~cellfun(@isempty, binnedCTD)); % Prune empty bins
+    profileBinned = profileBinned(~cellfun(@isempty, profileBinned)); % Prune empty bins
+    dissBinned = dissBinned(~cellfun(@isempty, dissBinned)); % Prune empty bins
+    chiBinned = chiBinned(~cellfun(@isempty, chiBinned)); % Prune empty bins
+    ctdBinned = ctdBinned(~cellfun(@isempty, ctdBinned)); % Prune empty bins
 
-    qProf = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), binnedProfile); % Valid data for profile depth binning
-    qDiss = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), binnedDiss);    % Valid data for diss depth binning
-    qCTD  = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), binnedCTD);     % Valid data for CTD time binning
+    qProf = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), profileBinned); % Valid data for profile depth binning
+    qDiss = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), dissBinned);    % Valid data for diss depth binning
+    qChi  = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), chiBinned);     % Valid data for chi depth binning
+    qCTD  = ~cellfun(@(x) isempty(x{1}) || ismissing(x{1}), ctdBinned);     % Valid data for CTD time binning
 
     if any(qProf)
-        profiles2combo(binnedProfile(qProf), pars);
+        profiles2combo(profileBinned(qProf), pars);
     end % if any qProf
     if any(qDiss)
-        diss2combo(binnedDiss(qDiss), pars);
+        diss2combo(dissBinned(qDiss), pars);
+    end % if any qDiss
+    if any(qChi)
+        chi2combo(chiBinned(qChi), pars);
     end % if any qDiss
     if ~ismissing(pars.CT_T_name) && ~ismissing(pars.CT_C_name) && any(qCTD)
-        ctd2combo(binnedCTD(qCTD), pars);
+        ctd2combo(ctdBinned(qCTD), pars);
     end % if any qCTD
 
     fprintf("\n********* Finished at %s in %.0f seconds **********\n", datetime(), toc(stime));
@@ -139,21 +150,23 @@
 warning(warningState); % Restore the warning status
 end % process_P_files
 
-function [row, binnedProfile, binnedDiss, ctdBinned] = P_to_binned_profile(row, pars)
+function [row, profileBinned, dissBinned, chiBinned, ctdBinned] = P_to_binned_profile(row, pars)
 arguments (Input)
     row (1,:) table
     pars struct
 end % arguments Input
 arguments (Output)
     row (1,:) table
-    binnedProfile (2,1) cell % {filename, data}
-    binnedDiss (2,1) cell
-    ctdBinned    (2,1) cell
+    profileBinned (2,1) cell % {filename, data}
+    dissBinned (2,1) cell
+    chiBinned (2,1) cell
+    ctdBinned  (2,1) cell
 end % arguments Input
 
-binnedProfile = {missing, []}; % Filename of binned profile information, binned profile data
-binnedDiss = {missing, []};
+profileBinned = {missing, []}; % Filename of binned profile information, binned profile data
+dissBinned = {missing, []};
 ctdBinned = {missing, []};
+chiBinned = {missing, []};
 gps = [];
 
 [row, mat] = convert2mat(row, pars); % Convert P file to mat via odas_p2mat
@@ -170,11 +183,17 @@
 
 if ~row.qProfileOkay, return; end % Nothing more to do
 
-[row, binnedProfile] = profile2binned(row, profiles, pars); % Bin profiles by depth
+[row, profileBinned] = profile2binned(row, profiles, pars); % Bin profiles by depth
 
 % Calculate the dissipations for each profile
 [row, diss] = profile2diss(row, profiles, pars); % Calculate dissipations
-[row, binnedDiss] = diss2binned(row, diss, pars); % Bin the dissipation
+[row, dissBinned] = diss2binned(row, diss, pars); % Bin the dissipation
+
+% Calculate Chi
+if pars.chi_enable
+    [row, chi] = profile2chi(row, profiles, diss, pars); % Calculate Chi for each profile
+    [row, chiBinned] = chi2binned(row, chi, pars); % Bin the chi estimates
+end % if chi
 end %s process_P_to_binned_profile
 
 function startProcessPool(pars)
diff --git a/Code/profile2chi.m b/Code/profile2chi.m
new file mode 100644
index 0000000..3791a72
--- /dev/null
+++ b/Code/profile2chi.m
@@ -0,0 +1,73 @@
+%
+% For each profile, calculate dissipation estimates
+%
+% Nov-2023, Pat Welch, pat@mousebrains.com
+
+function [row, chiInfo] = profile2chi(row, profileInfo, dissInfo, pars)
+arguments (Input)
+    row (1,:) table % row to work on
+    profileInfo struct % Output of mat2profile
+    dissInfo struct % Output of profile2diss
+    pars struct % Parameters, defaults from get_info
+end % arguments Input
+arguments (Output)
+    row table % row worked on
+    chiInfo  % chi estimates structure or empty
+end % arguments Output
+
+chiInfo = [];
+
+fnChi = fullfile(pars.chi_root, append(row.name, ".mat"));
+row.fnChi = fnChi;
+
+if isnewer(fnChi, row.fnDiss) % fnChi is newer than fnDiss
+    fprintf("%s: %s newer than %s\n", row.name, fnChi, row.fnDiss);
+    return;
+end
+
+if isempty(profileInfo) % we need to load profileInfo
+    fprintf("Loading %s\n", row.fnProf);
+    profileInfo = load(row.fnProf);
+end
+
+if isempty(dissInfo) % we need to load dissInfo
+    fprintf("Loading %s\n", row.fnDiss);
+    dissInfo = load(row.fnDiss);
+end
+
+warning("Chi not implemented");
+% return;
+% 
+% stime = tic();
+% 
+% profiles = profileInfo.profiles;
+% pInfo = profileInfo.pInfo;
+% nProfiles = numel(profiles);
+% 
+% tbl = cell(nProfiles, 1);
+% dInfo = cell(size(tbl));
+% 
+% for index = 1:nProfiles
+%     [dInfo{index}, tbl{index}] = calc_diss_shear(profiles{index}, pInfo(index,:), pars);
+% end % for index
+% 
+% qEmpty = cellfun(@isempty, tbl);
+% if all(qEmpty) % No dissipation estimates
+%     fprintf("%s: No dissipation estimates found", row.name);
+%     return;
+% end
+% 
+% tbl = tbl(~qEmpty);
+% dInfo = dInfo(~qEmpty);
+% 
+% dissInfo = struct();
+% dissInfo.info = vertcat(dInfo{:});
+% dissInfo.profiles = tbl;
+% if isfield(profileInfo, "fp07Lags")
+%     dissInfo.fp07 = profileInfo.fp07Lags;
+% end % if isfield
+% 
+% my_mk_directory(fnChi);
+% save(fnChi, "-struct", "chiInfo", pars.matlab_file_format);
+% fprintf("Took %.2f seconds to create %s\n", toc(stime), fnChi);
+end % profile2diss
\ No newline at end of file
diff --git a/Code/update_paths.m b/Code/update_paths.m
index ea6a828..866e9ce 100644
--- a/Code/update_paths.m
+++ b/Code/update_paths.m
@@ -34,6 +34,10 @@
     | startsWith(names, "diss_"); % Dissipation parameters
 qDissBin = qDiss ...
     | startsWith(names, "binDiss_"); % how to bin dissipation
+qChi = qDiss ...
+    | startsWith(names, "chi_"); % Chi related parameters
+qChiBin = qChi ...
+    | startsWith(names, "binChi_"); % how to bin chi
 qCTD = qP2Mat ...
     | (qProfileGen & (pars.profile_direction == "down")) ...
     | qCT ...
@@ -51,6 +55,9 @@
 [hashDiss,      jsonDiss]      = mk_hash_json(pars, names(qDiss));
 [hashDissBin,   jsonDissBin]   = mk_hash_json(pars, names(qDissBin));
 [hashDissCombo, jsonDissCombo] = mk_hash_json(pars, names(qDissBin | qPfiles | qNetCDF));
+[hashChi,       jsonChi]       = mk_hash_json(pars, names(qChi));
+[hashChiBin,    jsonChiBin]    = mk_hash_json(pars, names(qChiBin));
+[hashChiCombo,  jsonChiCombo]  = mk_hash_json(pars, names(qChiBin | qPfiles | qNetCDF));
 
 pars.output_root = abspath(pars.output_root); % Get rid of ~ or relative paths
 my_mk_directory(pars.output_root); % Make sure the root path exists
@@ -68,6 +75,11 @@
 pars.diss_root        = mkRootDir(pars.output_root, "diss", hashDiss, jsonDiss);
 pars.diss_binned_root = mkRootDir(pars.output_root, "diss_binned", hashDissBin, jsonDissBin);
 pars.diss_combo_root  = mkRootDir(pars.output_root, "diss_combo", hashDissCombo, jsonDissCombo);
+if pars.chi_enable
+    pars.chi_root = mkRootDir(pars.output_root, "chi", hashChi, jsonChi);
+    pars.chi_binned_root = mkRootDir(pars.output_root, "chi_binned", hashChiBin, jsonChiBin);
+    pars.chi_combo_root  = mkRootDir(pars.output_root, "chi_combo", hashChiCombo, jsonChiCombo);
+end % if chi
 
 pars.log_root = fullfile(pars.output_root, "logs"); % Where to write log files to
 pars.database_root = fullfile(pars.output_root, "database"); % Where to store various databases

From 3c8998d9321b8bf63e4cdcd77e7f68e0664015c5 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Wed, 7 Aug 2024 17:24:13 -0700
Subject: [PATCH 3/8] Associate instrument serial number with each profile

---
 Code/CTD.json      | 5 +++++
 Code/Combo.json    | 5 +++++
 Code/convert2mat.m | 8 ++++++++
 Code/ctd2binned.m  | 4 ++--
 Code/mat2profile.m | 3 ++-
 5 files changed, 22 insertions(+), 3 deletions(-)

diff --git a/Code/CTD.json b/Code/CTD.json
index 5f281bb..f8e2409 100644
--- a/Code/CTD.json
+++ b/Code/CTD.json
@@ -200,6 +200,11 @@
 			"long_name": "Turbidity",
 			"standard_name": "sea_water_turbidity",
 			"coverage_content_type": "physicalMeasurement"
+		},
+		"sn": {
+			"units": "",
+			"long_name": "Serial Number",
+			"coverage_content_type": "auxiliaryInformation"
 		}
 	},
 
diff --git a/Code/Combo.json b/Code/Combo.json
index 3099bfa..2ee24f1 100644
--- a/Code/Combo.json
+++ b/Code/Combo.json
@@ -546,6 +546,11 @@
 			"units": "",
 			"long_name": "ODAS Library Version",
 			"coverage_content_type": "auxiliaryInformation"
+		},
+		"sn": {
+			"units": "",
+			"long_name": "Serial Number",
+			"coverage_content_type": "auxiliaryInformation"
 		}
 	},
 
diff --git a/Code/convert2mat.m b/Code/convert2mat.m
index 2fa214e..d7bc497 100644
--- a/Code/convert2mat.m
+++ b/Code/convert2mat.m
@@ -36,6 +36,14 @@
     a = odas_p2mat(char(row.fn), p2args{:}); % extract P file contents
     my_mk_directory(row.fnMat);
     save(row.fnMat, "-struct", "a", pars.matlab_file_format); % save into a mat file
+
+    row.sn = setupstr(a.setupfilestr, 'instrument_info', 'sn'); % Get the serial number in the instrument_info stanza
+    if isempty(row.sn)
+        row.sn = missing;
+    else
+        row.sn = string(row.sn{1});
+    end % if isempty
+
     row.qMatOkay = true;
     fprintf("Took %.2f seconds to convert %s\n", toc(stime), row.name);
 catch ME
diff --git a/Code/ctd2binned.m b/Code/ctd2binned.m
index 90a2779..3a5c0ee 100644
--- a/Code/ctd2binned.m
+++ b/Code/ctd2binned.m
@@ -144,7 +144,7 @@
 fnCTD = fullfile(pars.ctd_root, append(row.name, ".mat"));
 row.fnCTD = fnCTD;
 
-cInfo = row(:,["name", "t0", "tEnd"]);
+cInfo = row(:,["name", "t0", "tEnd", "sn"]);
 cInfo = renamevars(cInfo, "tEnd", "t1"); % For NetCDF time range
 
 binned = struct("tbl", tbl, "info", cInfo);
@@ -216,4 +216,4 @@
     ctd.lon(q) = gps.lon(t) - dLondt .* dt;
     ctd.lat(q) = gps.lat(t) - dLatdt .* dt;
 end
-end % addGPS
\ No newline at end of file
+end % addGPS
diff --git a/Code/mat2profile.m b/Code/mat2profile.m
index d4f03cc..7fa7669 100644
--- a/Code/mat2profile.m
+++ b/Code/mat2profile.m
@@ -157,7 +157,7 @@
 
 profilesInfo = struct();
 profilesInfo.profiles = profiles;
-profilesInfo.row = row(1,["name", "date", "fClock" "t0", "t1", "tEnd"]);
+profilesInfo.row = row(1,["name", "date", "fClock" "t0", "t1", "tEnd", "sn"]);
 profilesInfo.pInfo = profileInfo;
 if qFP07
     profilesInfo.fp07Lags = fp07_lags;
@@ -181,6 +181,7 @@
 tbl = table();
 tbl.name = repmat(row.name, nProfiles, 1);
 tbl.index = (1:nProfiles)';
+tbl.sn = repmat(row.sn, nProfiles,1);
 tbl.t0 = NaT(nProfiles,1);
 tbl.t1 = NaT(nProfiles,1);
 tbl.n_slow = nan(nProfiles,1);

From e4da96f23fb2c932feaebea3c88c208403c20a44 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Sat, 10 Aug 2024 13:02:11 -0700
Subject: [PATCH 4/8] First cut at some Chi code, which I'm not happy with

---
 Code/calc_chi.m     |  35 +++++++++++++
 Code/chi2binned.m   |  47 ++++++++---------
 Code/profile2chi.m  | 119 +++++++++++++++++++++++++++++---------------
 Examples/example0.m |   1 +
 4 files changed, 139 insertions(+), 63 deletions(-)
 create mode 100644 Code/calc_chi.m

diff --git a/Code/calc_chi.m b/Code/calc_chi.m
new file mode 100644
index 0000000..3b8c40d
--- /dev/null
+++ b/Code/calc_chi.m
@@ -0,0 +1,35 @@
+%
+% Calculate chi from temperature gradients given an already calculated epsilon
+%
+% Aug-2024, Pat Welch, pat@mousebrains.com
+
+function [dInfo, tbl] = calc_chi(diss, dInfo)
+arguments (Input)
+    diss table % Profile information
+    dInfo (1,:) table % Summary information about the profile
+end % arguments Input
+arguments (Output)
+    dInfo (1,:) table % pInfo with extra fields
+    tbl table % Tabular form of diss struct
+end % arguments Output
+
+kappa_T = 1.4e-7; % thermal diffusivity [m^2/s]
+
+tbl = table();
+tbl.t = diss.t;
+tbl.depth = diss.depth;
+
+names = string(diss.Properties.VariableNames);
+
+for name = names(startsWith(names, "gradT"))
+    gradT = diss.(name);
+    gradT2 = gradT.^2;
+    % gradT2(gradT <= 0) = NaN;
+    index = extractAfter(name, "gradT");
+    tbl.(name) = gradT;
+    tbl.(sprintf("chi_dT%s_mean", index)) = 2 * kappa_T * diss.N2 .* diss.epsilonMean ./ gradT2;
+    for j = 1:size(diss.e,2)
+        tbl.(sprintf("chi_dT%s_e%d", index, j)) = 2 * kappa_T * diss.N2 .* diss.e(:,j) ./ gradT2;
+    end % for j
+end % for name
+end % calc_chi
\ No newline at end of file
diff --git a/Code/chi2binned.m b/Code/chi2binned.m
index dd753dd..fa92623 100644
--- a/Code/chi2binned.m
+++ b/Code/chi2binned.m
@@ -19,30 +19,31 @@
     return;
 end % if ~row.qProfileOkay
 
-warning("Chi binning not implemented");
+fnChi = row.fnChi;
+fnBin = fullfile(pars.chi_binned_root, append(row.name, ".mat"));
+row.fnChiBin = fnBin;
 
-% fnChi = row.fnChi;
-% fnBin = fullfile(pars.chi_binned_root, append(row.name, ".mat"));
-% row.fnChiBin = fnBin;
-% 
-% if isnewer(fnBin, fnChi)
-%     retval = {fnBin, []}; % We want this for combining
-%     fprintf("%s: %s is newer than %s\n", row.name, row.fnBin, row.fnChi);
-%     return;
-% end % if isnewer
-% 
-% if isempty(a)
-%     fprintf("Loading %s\n", row.fnChi);
-%     a = load(row.fnChi);
-% end % if isempty
-% 
-% %% Bin the data into depth bins
-% 
-% pInfo = a.info;
-% profiles = a.profiles;
-% 
-% fprintf("%s: Binning %d profiles\n", row.name, numel(profiles));
-% 
+if isnewer(fnBin, fnChi)
+    retval = {fnBin, []}; % We want this for combining
+    fprintf("%s: %s is newer than %s\n", row.name, row.fnBin, row.fnChi);
+    return;
+end % if isnewer
+
+if isempty(a)
+    fprintf("Loading %s\n", row.fnChi);
+    a = load(row.fnChi);
+end % if isempty
+
+%% Bin the data into depth bins
+
+pInfo = a.info;
+profiles = a.profiles;
+
+fprintf("%s: Binning %d profiles\n", row.name, numel(profiles));
+
+pInfo
+profiles{1}
+error("GotMe");
 % if pars.profile_direction == "time" % Bin in time
 %     binSize = seconds(pars.binChi_width); % Bin stepsize in (sec)
 %     keyName = "t";
diff --git a/Code/profile2chi.m b/Code/profile2chi.m
index 3791a72..1238fe3 100644
--- a/Code/profile2chi.m
+++ b/Code/profile2chi.m
@@ -1,13 +1,13 @@
 %
-% For each profile, calculate dissipation estimates
+% For each profile, calculate chi estimates
 %
 % Nov-2023, Pat Welch, pat@mousebrains.com
 
-function [row, chiInfo] = profile2chi(row, profileInfo, dissInfo, pars)
+function [row, chiInfo] = profile2chi(row, profileInfo, chiInfo, pars)
 arguments (Input)
     row (1,:) table % row to work on
     profileInfo struct % Output of mat2profile
-    dissInfo struct % Output of profile2diss
+    chiInfo struct % Output of profile2diss
     pars struct % Parameters, defaults from get_info
 end % arguments Input
 arguments (Output)
@@ -30,44 +30,83 @@
     profileInfo = load(row.fnProf);
 end
 
-if isempty(dissInfo) % we need to load dissInfo
+if isempty(chiInfo) % we need to load dissInfo
     fprintf("Loading %s\n", row.fnDiss);
-    dissInfo = load(row.fnDiss);
+    chiInfo = load(row.fnDiss);
 end
 
-warning("Chi not implemented");
-% return;
-% 
-% stime = tic();
-% 
-% profiles = profileInfo.profiles;
-% pInfo = profileInfo.pInfo;
-% nProfiles = numel(profiles);
-% 
-% tbl = cell(nProfiles, 1);
-% dInfo = cell(size(tbl));
-% 
-% for index = 1:nProfiles
-%     [dInfo{index}, tbl{index}] = calc_diss_shear(profiles{index}, pInfo(index,:), pars);
-% end % for index
-% 
-% qEmpty = cellfun(@isempty, tbl);
-% if all(qEmpty) % No dissipation estimates
-%     fprintf("%s: No dissipation estimates found", row.name);
-%     return;
-% end
-% 
-% tbl = tbl(~qEmpty);
-% dInfo = dInfo(~qEmpty);
-% 
-% dissInfo = struct();
-% dissInfo.info = vertcat(dInfo{:});
-% dissInfo.profiles = tbl;
-% if isfield(profileInfo, "fp07Lags")
-%     dissInfo.fp07 = profileInfo.fp07Lags;
-% end % if isfield
-% 
-% my_mk_directory(fnChi);
-% save(fnChi, "-struct", "chiInfo", pars.matlab_file_format);
-% fprintf("Took %.2f seconds to create %s\n", toc(stime), fnChi);
+stime = tic();
+
+profiles = profileInfo.profiles;
+pInfo = profileInfo.pInfo;
+nProfiles = numel(profiles);
+
+dProfiles = chiInfo.profiles;
+dInfo = chiInfo.info;
+
+if size(pInfo,1) ~= size(dInfo,1)
+    warning("Mismatched number of profiles in chi calculation, %d ~= %d", size(pInfo,1), size(dInfo,1));
+    return;
+end % if size
+
+tbl = cell(nProfiles, 1);
+dInfo = cell(size(tbl));
+
+for index = 1:nProfiles
+    profile = profiles{index};
+    fast = profile.fast;
+    slow = profile.slow;
+    diss = dProfiles{index};
+    fast.dDepth = interp1(diss.depth, diss.depth, fast.depth, "nearest", "extrap");
+    fast.grp = findgroups(fast.dDepth);
+    a = rowfun(@(x) x(1), fast, ...
+        InputVariables="dDepth", ...
+        GroupingVariables="grp", ...
+        OutputVariableNames="depth");
+    names = string(fast.Properties.VariableNames);
+    names = names(startsWith(names, "gradT"));
+    for name = names
+        a.(name) = rowfun(@(x) median(x, "omitnan"), fast, ...
+            InputVariables=name, ...
+            GroupingVariables="grp", ...
+            OutputFormat="uniform");
+    end % for
+
+    N2 = table();
+    [N2.N2, N2.pMid] = gsw_Nsquared(slow.SA, slow.theta, slow.P_slow, slow.lat);
+    N2.depth = (slow.depth(1:end-1) + slow.depth(2:end)) / 2;
+    N2.dDepth = interp1(diss.depth, diss.depth, N2.depth, "nearest", "extrap");
+    N2.grp = findgroups(N2.dDepth);
+
+    b = rowfun(@(x) x(1), N2, ...
+        InputVariables="dDepth", ...
+        GroupingVariables="grp", ...
+        OutputVariableNames="depth");
+    b.N2 = rowfun(@(x) median(x, "omitnan"), N2, ...
+            InputVariables="N2", ...
+            GroupingVariables="grp", ...
+            OutputFormat="uniform");
+
+    diss = innerjoin(diss, a, Keys="depth", RightVariables=names);
+    diss = innerjoin(diss, b, Keys="depth", RightVariables="N2");
+
+    [dInfo{index}, tbl{index}] = calc_chi(diss, pInfo(index,:));
+end % for index
+
+qEmpty = cellfun(@isempty, tbl);
+if all(qEmpty) % No dissipation estimates
+    fprintf("%s: No dissipation estimates found", row.name);
+    return;
+end
+
+tbl = tbl(~qEmpty);
+dInfo = dInfo(~qEmpty);
+
+chiInfo = struct();
+chiInfo.info = vertcat(dInfo{:});
+chiInfo.profiles = tbl;
+
+my_mk_directory(fnChi);
+save(fnChi, "-struct", "chiInfo", pars.matlab_file_format);
+fprintf("Took %.2f seconds to create %s\n", toc(stime), fnChi);
 end % profile2diss
\ No newline at end of file
diff --git a/Examples/example0.m b/Examples/example0.m
index 99511bf..31b58c4 100644
--- a/Examples/example0.m
+++ b/Examples/example0.m
@@ -17,6 +17,7 @@
 
     pars = process_P_files( ...
         "debug", true, ...
+        "chi_enable", true, ... % Calculate chi
         "p_file_root", p_file_root, ... % Where the input .P files are located
         "p_file_pattern", "SN*/*", ...  % Glob pattern appended to p_file_root to locate P files
         "output_root", output_root, ... % Where to write output to

From 4252b33f07dd550404c8de71a059ca53269a3d6e Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Sat, 10 Aug 2024 13:34:16 -0700
Subject: [PATCH 5/8] Fix typo

---
 Code/profile2chi.m | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Code/profile2chi.m b/Code/profile2chi.m
index 48be7ab..97ff8c8 100644
--- a/Code/profile2chi.m
+++ b/Code/profile2chi.m
@@ -4,11 +4,11 @@
 %
 % Nov-2023, Pat Welch, pat@mousebrains.com
 
-function [row, chiInfo] = profile2chi(row, profileInfo, chiInfo, pars)
+function [row, chiInfo] = profile2chi(row, profileInfo, dissInfo, pars)
 arguments (Input)
     row (1,:) table % row to work on
     profileInfo struct % Output of mat2profile
-    chiInfo struct % Output of profile2diss
+    dissInfo struct % Dissipation for this this profile
     pars struct % Parameters, defaults from get_info
 end % arguments Input
 arguments (Output)
@@ -42,8 +42,8 @@
 pInfo = profileInfo.pInfo;
 nProfiles = numel(profiles);
 
-dProfiles = chiInfo.profiles;
-dInfo = chiInfo.info;
+dProfiles = dissInfo.profiles;
+dInfo = dissInfo.info;
 
 if size(pInfo,1) ~= size(dInfo,1)
     warning("Mismatched number of profiles in chi calculation, %d ~= %d", size(pInfo,1), size(dInfo,1));

From 19e498ae6f490f56393ec32a9802bb6530dd5051 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Sat, 10 Aug 2024 13:34:32 -0700
Subject: [PATCH 6/8] make chi binning work

---
 Code/chi2binned.m | 109 ++++++++++++++++++++--------------------------
 1 file changed, 47 insertions(+), 62 deletions(-)

diff --git a/Code/chi2binned.m b/Code/chi2binned.m
index 658873e..58fec1a 100644
--- a/Code/chi2binned.m
+++ b/Code/chi2binned.m
@@ -41,66 +41,51 @@
 
 fprintf("%s: Binning %d profiles\n", row.name, numel(profiles));
 
-pInfo
-profiles{1}
-error("GotMe");
-
-% if pars.profile_direction == "time" % Bin in time
-%     binSize = seconds(pars.binChi_width); % Bin stepsize in (sec)
-%     keyName = "t";
-%     binFunc = @bin_by_time;
-%     glueFunc = @glue_lengthwise;
-% else % Bin by depth
-%     binSize = pars.binChi_width; % Bin stepsize (m)
-%     keyName = pars.binChi_variable; % Variable to bin by
-%     binFunc = @bin_by_real;
-%     glueFunc = @glue_widthwise;
-% end % if profile_direction
-% 
-% casts = cell(numel(profiles),1);
-% for index = 1:numel(profiles)
-%     profile = profiles{index};
-%     nE = size(profile.e, 2);
-%     prof2 = table();
-%     for name = string(profile.Properties.VariableNames)
-%         sz = size(profile.(name),2);
-%         if ~ismatrix(profile.(name)), continue; end
-%         if sz == 1
-%             prof2.(name) = profile.(name);
-%         elseif sz == nE
-%             for j = 1:sz
-%                 prof2.(append(name, "_", string(j))) = profile.(name)(:,j);
-%             end % for j;
-%         end % if sz
-%     end % for name
-% 
-%     casts{index} = binFunc(binSize, keyName, prof2, pars.binChi_method);
-% end % for index
-% 
-% qDrop = cellfun(@isempty, casts); % This shouldn't happend
-% 
-% if any(qDrop)
-%     casts = casts(~qDrop);
-%     pInfo = pInfo(~qDrop,:);
-% end % any qDrop
-% 
-% if isempty(casts)
-%     row.qProfileOkay = false;
-%     fprintf("%s: No usable casts found in %s\n", row.name, row.fnProf);
-%     return;
-% end
-% 
-% tbl = glueFunc("bin", casts);
-% 
-% binned = struct ( ...
-%     "tbl", tbl, ...
-%     "info", pInfo);
-% if isfield(a, "fp07")
-%     binned.fp07 = a.fp07;
-% end % if isfield
-% 
-% my_mk_directory(fnBin);
-% save(fnBin, "-struct", "binned", pars.matlab_file_format);
-% fprintf("%s: Saving %d profiles to %s\n", row.name, size(binned.info,1), fnBin);
-% retval = {fnBin, binned};
+if pars.profile_direction == "time" % Bin in time
+    binSize = seconds(pars.binChi_width); % Bin stepsize in (sec)
+    keyName = "t";
+    binFunc = @bin_by_time;
+    glueFunc = @glue_lengthwise;
+else % Bin by depth
+    binSize = pars.binChi_width; % Bin stepsize (m)
+    keyName = pars.binChi_variable; % Variable to bin by
+    binFunc = @bin_by_real;
+    glueFunc = @glue_widthwise;
+end % if profile_direction
+
+casts = cell(numel(profiles),1);
+for index = 1:numel(profiles)
+    profile = profiles{index};
+    prof2 = table();
+    for name = string(profile.Properties.VariableNames)
+        if ~ismatrix(profile.(name)), continue; end
+        prof2.(name) = profile.(name);
+    end % for name
+
+    casts{index} = binFunc(binSize, keyName, prof2, pars.binChi_method);
+end % for index
+
+qDrop = cellfun(@isempty, casts); % This shouldn't happend
+
+if any(qDrop)
+    casts = casts(~qDrop);
+    pInfo = pInfo(~qDrop,:);
+end % any qDrop
+
+if isempty(casts)
+    row.qProfileOkay = false;
+    fprintf("%s: No usable casts found in %s\n", row.name, row.fnProf);
+    return;
+end
+
+tbl = glueFunc("bin", casts);
+
+binned = struct ( ...
+    "tbl", tbl, ...
+    "info", pInfo);
+
+my_mk_directory(fnBin);
+save(fnBin, "-struct", "binned", pars.matlab_file_format);
+fprintf("%s: Saving %d profiles to %s\n", row.name, size(binned.info,1), fnBin);
+retval = {fnBin, binned};
 end % chi2binned

From ae73c3aaf09204df4cb0538906b18635551f057e Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Mon, 2 Sep 2024 09:42:45 -0700
Subject: [PATCH 7/8] Update code to reference version

---
 Code/osgl_get_netCDF.m | 81 +++++++++++++++++++++++++++++++-----------
 1 file changed, 60 insertions(+), 21 deletions(-)

diff --git a/Code/osgl_get_netCDF.m b/Code/osgl_get_netCDF.m
index 04c2811..c58c85d 100644
--- a/Code/osgl_get_netCDF.m
+++ b/Code/osgl_get_netCDF.m
@@ -15,27 +15,21 @@
 %                      for speed reasons
 % Dec-2022, Pat Welch, pat@mousebrains.com, add in compliance with CF-1.8
 %                      metadata standards
-%
+% Jan-2024, Pat Welch, pat@mosuebrains.com, straighten out handling of variable names
 
-function a = osgl_get_netCDF(fn, varargin)
+function a = osgl_get_netCDF(fn, names)
 arguments (Input)
     fn string {mustBeFile} % Input filename
 end % arguments Input
 arguments (Input,Repeating)
-    varargin % Optional arguments
+    names (:,1) string % Optional variable names
 end % arguments Repeating
 arguments (Output)
     a struct % Output structure of data read from fn
 end % arguments output
 
-names = cell(numel(varargin),1); % Pre-allocate
-
-for i = 1:numel(varargin)
-    item = string(varargin{i});
-    names{i} = item(:);
-end % for i
-
-names = string(names); % Convert to an array of strings from a cell array
+names = cellfun(@(x) x(:), names, "UniformOutput", false); % Make sure any arrays are in the same orientation
+names = string(vertcat(names{:}));
 
 ncid = netcdf.open(fn); % Open the netCDF file
 cleanUpObj = onCleanup(@() netcdf.close(ncid)); % In case of errors/warnings/...
@@ -43,6 +37,7 @@
 if isempty(names)
     [names, varids] = loadNames(ncid);
 else
+    names = unique(names); % No duplicates
     [names, varids] = getVarIDs(ncid, names);
 end % if
 
@@ -70,6 +65,15 @@
 end % osgl_get_netCDF
 
 function data = modifyCF(data, ncid, varid)
+arguments (Input)
+    data % Input data column from NetCDF file
+    ncid double % NetCDF id from netcdf.open
+    varid double % Variable id within ncid
+end % arguments Output
+arguments (Output)
+    data % Possibly modified version of input data
+end % arguments Output
+
 qOkay = true(size(data)); % Initially assume all the data is good
 
 try
@@ -137,19 +141,32 @@
 end % modifyCF
 
 function data = ctConvert(data, ncid, varid, dtUnits, timeStr, qOkay)
+arguments (Input)
+    data         % Input from NetCDF file various formats
+    ncid double  % NetCDF id from netcdf.open
+    varid double % Variable id within ncid
+    dtUnits string % Time units from units attribute (The part before "since")
+    timeStr string % Time reference from units attribute (The part after "since")
+    qOkay logical % % Which data entries are good
+end % arguments Input
+arguments (Output)
+    data % possibly datetime, if converted
+end % arguments Output
+
 try
-    calendar = netcdf.getAtt(ncid, varid, "calendar");
-    goodCalendars = ["", "standard", "gregorian", "proleptic_gregorian"];
-    if ~ismember(calendar, goodCalendars)
-        warning("Unsupported calendar, %s, for %s", ...
-            calendar, netcdf.inqVar(ncid, varid));
-        return
-    end % if
+    calendar = netcdf.getAtt(ncid, varid, "calendar"); % Should be there, but somem people don't have it
 catch ME
     getReport(ME)
-    % Do nothing
+    calendar = /Users/pat/Desktop/Chi/perturb/Code/osgl_get_netCDF.m"; % Default
 end
 
+goodCalendars = ["", "standard", "gregorian", "proleptic_gregorian"];
+if ~ismember(calendar, goodCalendars)
+    warning("Unsupported calendar, %s, for %s", ...
+	    calendar, netcdf.inqVar(ncid, varid));
+    return
+end % if
+
 switch lower(dtUnits) % The part before since
     case {"years", "year", "y"}
         dt = data(qOkay) * 365.242198781;
@@ -175,11 +192,16 @@
         return
 end % switch
 
-fmt = "^\s*(\d{4})-(\d{1,2})-(\d{1,2})(\s+\d{1,2}:\d{1,2}:\d{1,2}([.]\d*|)(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)|)\s*$";
-tfmt = "^\s+(\d{1,2}):(\d{1,2}):(\d{1,2}([.]?\d*|))(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)$";
+fmt = strjoin([
+    "^(\d{4})-(\d{1,2})-(\d{1,2})", % yyyy-mm-dd
+    "(|(\s+|[Tt])\d{1,2}:\d{1,2}:\d{1,2}(|[.]\d*))", % (\s|T)HH:MM:SS.ssss
+    "(|\s+([A-Za-z/_]+|[+-]?\d{1,2}(|[:]?\d{1,2})))", % (UTC|[+-]00:00)
+    "\s*$"], "");
+tfmt = "^[Tt\s]+(\d{1,2}):(\d{1,2}):(\d{1,2}([.]?\d*|))(\s+([A-Za-z/_]+|[+-]?\d{1,2}([:]?\d{1,2}|))|)$";
 
 tokens = regexp(timeStr, fmt, "tokens", "once");
 if isempty(tokens)
+    fmt
     warning("Unable to parse time string, %s, for %s", ...
         timeStr, netcdf.inqVar(ncid, varid));
     return
@@ -228,6 +250,14 @@
 end % ctConvert
 
 function [names, varids] = loadNames(ncid)
+arguments (Input)
+    ncid double % NetCDF file id from netcdf.open
+end % arguments Input
+arguments (Output)
+    names (:,1) string  % Variable names in NetCDF file
+    varids (:,1) double % Variable ids in ncid
+end % arguments Output
+
 [~, n] = netcdf.inq(ncid); % Get the number of variables
 names = cell(n,1);
 varids = zeros(size(names)) - 1;
@@ -246,6 +276,15 @@
 end % loadNames
 
 function [names, varids] = getVarIDs(ncid, names)
+arguments (Input)
+    ncid double % NetCDF file id from netcdf.open
+    names (:,1) string  % Variable names to get from the NetCDF file
+end % arguments Input
+arguments (Output)
+    names (:,1) string  % Variable names to get from the NetCDF file
+    varids (:,1) double % Variable ids in ncid
+end % arguments Output
+
 varids = zeros(size(names));
 for i = 1:numel(names)
     varids(i) = netcdf.inqVarID(ncid, names(i));

From 9c2374f9663e6db249aae3a097206f4cd592a2d9 Mon Sep 17 00:00:00 2001
From: Pat Welch <mousebrains@users.noreply.github.com>
Date: Mon, 2 Sep 2024 09:46:37 -0700
Subject: [PATCH 8/8] Fix typo

---
 Code/osgl_get_netCDF.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Code/osgl_get_netCDF.m b/Code/osgl_get_netCDF.m
index c58c85d..30bec60 100644
--- a/Code/osgl_get_netCDF.m
+++ b/Code/osgl_get_netCDF.m
@@ -157,7 +157,7 @@
     calendar = netcdf.getAtt(ncid, varid, "calendar"); % Should be there, but somem people don't have it
 catch ME
     getReport(ME)
-    calendar = /Users/pat/Desktop/Chi/perturb/Code/osgl_get_netCDF.m"; % Default
+    calendar = "proleptic_gregorian"; % Default
 end
 
 goodCalendars = ["", "standard", "gregorian", "proleptic_gregorian"];