Added WIP publication scripts and hodogram plotting functions

PlanetMag.m

@@ -272,7 +272,7 @@
             end
 
             BD = ICAdecomposition(moonName, parentName, t_h*3600, BvecMoon, magModelDescrip, ...
-                SPHOUT, 1, 1, LIVE_PLOTS, figDir, figXtn);
+                SPHOUT, 1, 1, 1, LIVE_PLOTS, figDir, figXtn);
 
             T_h = 1 ./ BD.f / 3600;
             npeaks = length(T_h);

comparison/CompareSatModels.m

@@ -1,5 +1,5 @@
-function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, figXtn, ...
-    magDataMat, LOAD_PDS_ASCII)
+function CompareSatModels(LIVE_PLOTS, LOAD_PDS_ASCII, yearRange, scName, SEQUENTIAL, coeffPath, ...
+    figDir, figXtn, magDataMat)
 % Compare magnetic field measurements from spacecraft near Saturn against each implemented magnetic
 % field model.
 %
@@ -11,8 +11,13 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
 % ----------
 % LIVE_PLOTS : bool, default=0
 %   Whether to load interactive figure windows for plots (true) or print them to disk (false).
-% scNames : string, 1xS, default=["Voyager 1", "Voyager 2"]
-%   Spacecraft names for which magnetic field data will be compared against implemented models. A
+% LOAD_PDS_ASCII : bool, default=0
+%   Whether to load in .TAB files downloaded from PDS (true) or use pre-converted and compressed
+%   .mat summery of just the important bits (false).
+% yearRange : int, 1xG, default=4:17
+%   Years over which to compare Cassini data. Default includes all measurements.
+% scName : string, 1xS, default=["Cassini"]
+%   Spacecraft name for which magnetic field data will be compared against implemented models. A
 %   directory must exist with each name in the ``MAG`` directory within the top-level P\lanetMag
 %   directory. These directories will be searched for data files with the ``.tab`` extension, and
 %   each of these files will be loaded.
@@ -27,9 +32,6 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
 %   Extension to use for figures, which determines the file type.
 % magDataMat : char, 1xF, default='MAG/scName/ALL_FGM_KRTP_1M'
 %   File path to use for save/reload of PDS data in .mat file.
-% LOAD_PDS_ASCII : bool, default=0
-%   Whether to load in .TAB files downloaded from PDS (true) or use pre-converted and compressed
-%   .mat summery of just the important bits (false).
 
 % Part of the PlanetMag framework for evaluation and study of planetary magnetic fields.
 % Created by Corey J. Cochrane and Marshall J. Styczinski
@@ -38,14 +40,15 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     if ~exist('LIVE_PLOTS', 'var'); LIVE_PLOTS = 0; end
+    if ~exist('LOAD_PDS_ASCII', 'var'); LOAD_PDS_ASCII = 0; end
+    if ~exist('yearRange', 'var'); yearRange = 4:17; end
     if ~exist('scName', 'var'); scName = "Cassini"; end
     if ~exist('SEQUENTIAL', 'var'); SEQUENTIAL = 0; end
     if ~exist('coeffPath', 'var'); coeffPath = 'modelCoeffs'; end
     if ~exist('figDir', 'var'); figDir = 'figures'; end
     if ~exist('figXtn', 'var'); figXtn = 'pdf'; end
     sc = char(scName);
     if ~exist('magDataMat', 'var'); magDataMat = fullfile('MAG', sc, 'ALL_FGM_KRTP_1M'); end
-    if ~exist('LOAD_PDS_ASCII', 'var'); LOAD_PDS_ASCII = 0; end
 
     cspice_kclear;
     parentName = 'Saturn';
@@ -65,7 +68,7 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
         magFormatSpec = '%19s%11f%11f%11f%11f%8f%7f%7f%5f%3d%[^\n\r]';
         disp(['Importing PDS measurements over ' parentName ' orbits.'])
         [t_UTC, BrSC, BthSC, BphiSC] = deal([]);
-        for i=4:17
+        for i=yearRange
             datFile = fullfile('MAG', sc, [num2str(2000+i) '_FGM_KRTP_1M.TAB']);
             disp(['Loading ' sc ' MAG data from ' datFile '.'])
             fileID = fopen(datFile,'r');
@@ -120,7 +123,6 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
     BrSC(exclude) = [];
     BthSC(exclude) = [];
     BphiSC(exclude) = [];
-    rP_Rp(exclude) = [];
 
     npts = length(ets);
     t_h = ets / 3600;
@@ -140,16 +142,10 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
     end
 
     %% Plot trajectory
-    % Exclude points we won't be plotting
     disp(['Getting ' sc ' trajectories in KSO frame for ' num2str(npts) ' pts.'])
     [~, ~, ~, xyzKSO_km, ~] = GetPosSpice(sc, parentName, t_h, 'KSO');
     xyz_Rp = xyzKSO_km / Rp_km;
     x = xyz_Rp(1,:); y = xyz_Rp(2,:); z = xyz_Rp(3,:);
-    axlim_Rp = 10;
-    outer = find(rP_Rp > axlim_Rp);
-    x(outer) = nan;
-    y(outer) = nan;
-    z(outer) = nan;
 
     %% Plot trajectories
     windowName = 'Spacecraft Saturn trajectories';
@@ -166,6 +162,7 @@ function CompareSatModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, fig
 
     % Plot planet, pole, and rings for illustrative purposes
     pbaspect([1 1 1]);
+    axlim_Rp = 10;
     xlim([-axlim_Rp, axlim_Rp]); ylim([-axlim_Rp, axlim_Rp]); zlim([-axlim_Rp, axlim_Rp]);
     grid on;
     title([parentName ' orbit trajectories']);

functions/GetModelOpts.m

@@ -53,8 +53,8 @@
 %
 % MPopt : int
 %   Index number for magnetopause current magnetic field model. 0 selects the default option. If
-%   any value is passed that does not correspond to a valid model index, no magnetopause current
-%   magnetic field will be modeled. Currently implemented model options are:
+%   any value is passed that does not correspond to a valid model index (i.e. -1), no magnetopause
+%   current magnetic field will be modeled. Currently implemented model options are:
 %
 %       - ``Earth``: No magnetopause current magnetic field model is implemented.
 %       - ``Jupiter``:

functions/HodogramSingle.m

@@ -0,0 +1,151 @@
+function [xx, yy, windowName, titleInfo, xInfo, yInfo, fName, xlims, ylims] = HodogramSingle( ...
+    moonName, parentName, magModelDescrip, Bvec, SPHOUT, COMPARE_SEUFERT, COMPARE_PHIO)
+% Organizes data and labels for plotting a hodogram, showing the projection of the tip of the
+% magnetic field vector in a relevant plane, usually the equatorial plane.
+%
+% Must be plotted with a separate function.
+%
+% Parameters
+% ----------
+% moonName : char, 1xC
+%   Name of moon for which to evaluate excitation moments.
+% parentName : char, 1xD
+%   Name of parent planet for the moon for which to evaluate excitation moments.
+% magModelDescrip : char, 1xE
+%   Text description of the magnetic field model that was evalauted for the input time series.
+% Bvec : double, 3xN
+%   Magnetic field vectors at each measurement time.
+% SPHOUT : bool, default=0
+%   Whether to return vectors aligned to spherical coordinate axes (true) or cartesian (false).
+% COMPARE_SEUFERT : bool, default=1
+%   Whether to plot coordinate directions aligned to the axes presented in Seufert et al. (2011)
+%   https://doi.org/10.1016/j.icarus.2011.03.017. Only has an effect when ``SPHOUT`` is true,
+%   because Seufert et al. used spherical coordinates for evaluating the excitation spectra of
+%   Jupiter's moons.
+% COMPARE_PHIO : bool, default=1
+%   Whether to plot components in the same orientation as PhiO axes, as in past work, and label the
+%   axes with the comparisons (true), or just plot IAU frame axes without comparisons (false).
+%
+% Returns
+% -------
+% xx : double, 1xN
+%   x axis data.
+% yy : double, 1xN
+%   y axis data.
+% windowName : char, 1xC
+%   Name to use for the interactive figure window, which is shown when ``LIVE_PLOTS`` is true.
+% titleInfo : char, 1xD
+%   Description of plot contents to print at the top of the figure.
+% xInfo : char, 1xE
+%   Label for x axis of plot.
+% yInfo : char, 1xF
+%   Label for y axis of plot.
+% fName : char, 1xG
+%   File name to use for output figures, sans extension.
+% xlims : double, 1x2
+%   Minimum and maximum values to display on x axis, respectively.
+% ylims : double, 1x2
+%   Minimum and maximum values to display on y axis, respectively.
+
+% Part of the PlanetMag framework for evaluation and study of planetary magnetic fields.
+% Created by Corey J. Cochrane and Marshall J. Styczinski
+% Maintained by Marshall J. Styczinski
+% Contact: corey.j.cochrane@jpl.nasa.gov
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    if ~exist('SPHOUT', 'var'); SPHOUT = 0; end
+    if ~exist('COMPARE_SEUFERT', 'var'); COMPARE_SEUFERT = 1; end
+    if ~exist('COMPARE_PHIO', 'var'); COMPARE_PHIO = 1; end
+    % The following are defined in SetPlotDefaults. Do NOT reset them anywhere else.
+    global nmTxt
+    global bnmTxt
+    global mathTxt
+    global bmathTxt
+
+    windowName = [moonName ' ' magModelDescrip ' hodogram'];
+    if SPHOUT
+        if strcmp(parentName,'Saturn')
+            BcompMax = max(abs(Bvec(1:2,:)), [], 'all');
+            titleInfo = [bnmTxt moonName ' spin-parent plane hodogram, ' magModelDescrip];
+            if COMPARE_PHIO
+                xx = -Bvec(1,:);
+                yy = -Bvec(2,:);
+                xInfo = [mathTxt '-B_r ' nmTxt parentName ' SIII (' mathTxt '\approx B_y ' ...
+                    nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+                yInfo = [mathTxt '-B_\theta ' nmTxt parentName ' SIII (' mathTxt '\approx B_z ' ...
+                    nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+            else
+                xx = Bvec(1,:);
+                yy = Bvec(2,:);
+                xInfo = [mathTxt 'B_r ' nmTxt parentName ' SIII (nT)'];
+                yInfo = [mathTxt 'B_\theta ' nmTxt parentName ' SIII (nT)'];
+            end
+        else
+            BcompMax = max(abs([Bvec(1,:), Bvec(3,:)]));
+            titleInfo = [moonName ' equatorial plane hodogram, ' magModelDescrip];
+            if COMPARE_SEUFERT
+                xx = -Bvec(3,:);
+                yy =  Bvec(1,:);
+                if COMPARE_PHIO
+                    xInfo = [mathTxt '-B_\phi ' nmTxt parentName ' SIII (' mathTxt ...
+                        '\approx -B_x ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+                    yInfo = [mathTxt 'B_r ' nmTxt parentName ' SIII (' mathTxt ...
+                        '\approx -B_y ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+                else
+                    xInfo = [mathTxt '-B_\phi ' nmTxt parentName ' SIII (nT)'];
+                    yInfo = [mathTxt 'B_r ' nmTxt parentName ' SIII (nT)'];
+                end
+            else
+                if COMPARE_PHIO
+                    xx =  Bvec(3,:);
+                    yy = -Bvec(1,:);
+                    xInfo = [mathTxt 'B_\phi ' nmTxt parentName ' SIII (' mathTxt ...
+                        '\approx B_x ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+                    yInfo = [mathTxt '-B_r ' nmTxt parentName ' SIII (' mathTxt ...
+                        '\approx B_y ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
+                else
+                    xx = Bvec(1,:);
+                    yy = Bvec(3,:);
+                    xInfo = [mathTxt 'B_r ' nmTxt parentName ' SIII (nT)'];
+                    yInfo = [mathTxt 'B_\phi ' nmTxt parentName ' SIII (nT)'];
+                end
+            end
+        end
+    else
+        if strcmp(parentName,'Saturn')
+            BcompMax = max(abs([Bvec(1,:), Bvec(3,:)]));
+            titleInfo = [moonName ' spin-parent plane hodogram, ' magModelDescrip];
+            xx = Bvec(1,:);
+            yy = Bvec(3,:);
+            if COMPARE_PHIO
+                xInfo = [mathTxt 'B_x ' nmTxt 'IAU (' mathTxt '\approx B_y ' nmTxt moonName(1) ...
+                    mathTxt '\phi\Omega' nmTxt ', nT)'];
+                yInfo = [mathTxt 'B_z ' nmTxt 'IAU (' mathTxt '\approx B_z ' nmTxt moonName(1) ...
+                    mathTxt '\phi\Omega' nmTxt ', nT)'];
+            else
+                xInfo = [mathTxt 'B_x ' nmTxt 'IAU (nT)'];
+                yInfo = [mathTxt 'B_z ' nmTxt 'IAU (nT)'];
+            end
+        else
+            BcompMax = max(abs(Bvec(1:2,:)), [], 'all');
+            titleInfo = [moonName ' equatorial plane hodogram, ' magModelDescrip];
+            if COMPARE_PHIO
+                xx = Bvec(2,:);
+                yy = Bvec(1,:);
+                xInfo = [mathTxt 'B_y ' nmTxt 'IAU (' mathTxt '\approx -B_x ' nmTxt moonName(1) ...
+                    mathTxt '\phi\Omega' nmTxt ', nT)'];
+                yInfo = [mathTxt 'B_x ' nmTxt 'IAU (' mathTxt '\approx B_y ' nmTxt moonName(1) ...
+                    mathTxt '\phi\Omega' nmTxt ', nT)'];
+            else
+                xx = Bvec(1,:);
+                yy = Bvec(2,:);
+                xInfo = [mathTxt 'B_x ' nmTxt 'IAU (nT)'];
+                yInfo = [mathTxt 'B_y ' nmTxt 'IAU (nT)'];
+            end
+        end
+    end
+    ylims = [-BcompMax, BcompMax] * 1.05;
+    xlims = ylims;
+
+    fName = [moonName 'Hodogram' magModelDescrip];
+end

functions/ICAdecomposition.m

@@ -1,5 +1,5 @@
 function BD = ICAdecomposition(moonName, parentName, ets, Bvec, magModelDescrip, SPHOUT, ...
-    PLOT_DIAGNOSTIC, COMPARE_SEUFERT, LIVE_PLOTS, figDir, figXtn)
+    PLOT_DIAGNOSTIC, COMPARE_SEUFERT, COMPARE_PHIO, LIVE_PLOTS, figDir, figXtn)
 % Decomposes the input magnetic field vector time series into complex excitation moments using
 % :dfn:`Independent Component Analysis (ICA)`.
 %
@@ -59,6 +59,9 @@
 %   https://doi.org/10.1016/j.icarus.2011.03.017. Only has an effect when ``SPHOUT`` is true,
 %   because Seufert et al. used spherical coordinates for evaluating the excitation spectra of
 %   Jupiter's moons.
+% COMPARE_PHIO : bool, default=1
+%   Whether to plot components in the same orientation as PhiO axes, as in past work, and label the
+%   axes with the comparisons (true), or just plot IAU frame axes without comparisons (false).
 % LIVE_PLOTS : bool, default=0
 %   Whether to load interactive figure windows for plots (true) or print them to disk (false).
 % figDir : char, 1xF, default='figures'
@@ -100,8 +103,10 @@
     if ~exist('SPHOUT', 'var'); SPHOUT = 0; end
     if ~exist('PLOT_DIAGNOSTIC', 'var'); PLOT_DIAGNOSTIC = 1; end
     if ~exist('COMPARE_SEUFERT', 'var'); COMPARE_SEUFERT = 1; end
+    if ~exist('COMPARE_PHIO', 'var'); COMPARE_PHIO = 1; end
     if ~exist('LIVE_PLOTS', 'var'); LIVE_PLOTS = 0; end
     if ~exist('figDir', 'var'); figDir = 'figures'; end
+    if ~exist('figXtn', 'var'); figXtn = 'pdf'; end
     % The following are defined in SetPlotDefaults. Do NOT reset them anywhere else.
     global nmTxt
     global bnmTxt
@@ -258,64 +263,11 @@
             figDir, figXtn, LIVE_PLOTS, figNumBase + 2, 'linear', 'linear', xlims, 'auto', ...
             cfmt, 1, 1);
 
-        windowName = [moonName ' ' magModelDescrip ' hodogram'];
-        if SPHOUT
-            if strcmp(parentName,'Saturn')
-                xx = -Bvec(1,:);
-                yy = -Bvec(2,:);
-                BcompMax = max(abs(Bvec(1:2,:)), [], 'all');
-                xInfo = [mathTxt '-B_r ' nmTxt parentName ' SIII (' mathTxt '\approx B_y ' ...
-                    nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                yInfo = [mathTxt '-B_\theta ' nmTxt parentName ' SIII (' mathTxt '\approx B_z ' ...
-                    nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                titleInfo = [bnmTxt moonName ' spin-parent plane hodogram, ' magModelDescrip];
-            else
-                if COMPARE_SEUFERT
-                    xx = -Bvec(3,:);
-                    yy =  Bvec(1,:);
-                    xInfo = [mathTxt '-B_\phi ' nmTxt parentName ' SIII (' mathTxt ...
-                        '\approx -B_x ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                    yInfo = [mathTxt 'B_r ' nmTxt parentName ' SIII (' mathTxt ...
-                        '\approx -B_y ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                else
-                    xx =  Bvec(3,:);
-                    yy = -Bvec(1,:);
-                    xInfo = [mathTxt 'B_\phi ' nmTxt parentName ' SIII (' mathTxt ...
-                        '\approx B_x ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                    yInfo = [mathTxt '-B_r ' nmTxt parentName ' SIII (' mathTxt ...
-                        '\approx B_y ' nmTxt moonName(1) mathTxt '\phi\Omega' nmTxt ', nT)'];
-                end
-                BcompMax = max(abs([Bvec(1,:), Bvec(3,:)]));
-                titleInfo = [moonName ' equatorial plane hodogram, ' magModelDescrip];
-            end
-        else
-            if strcmp(parentName,'Saturn')
-                xx = Bvec(1,:);
-                yy = Bvec(3,:);
-                BcompMax = max(abs([Bvec(1,:), Bvec(3,:)]));
-                xInfo = [mathTxt 'B_x ' nmTxt 'IAU (' mathTxt '\approx B_y ' nmTxt moonName(1) ...
-                    mathTxt '\phi\Omega' nmTxt ', nT)'];
-                yInfo = [mathTxt 'B_z ' nmTxt 'IAU (' mathTxt '\approx B_z ' nmTxt moonName(1) ...
-                    mathTxt '\phi\Omega' nmTxt ', nT)'];
-                titleInfo = [moonName ' spin-parent plane hodogram, ' magModelDescrip];
-            else
-                xx = Bvec(2,:);
-                yy = Bvec(1,:);
-                BcompMax = max(abs(Bvec(1:2,:)), [], 'all');
-                xInfo = [mathTxt 'B_y ' nmTxt 'IAU (' mathTxt '\approx -B_x ' nmTxt moonName(1) ...
-                    mathTxt '\phi\Omega' nmTxt ', nT)'];
-                yInfo = [mathTxt 'B_x ' nmTxt 'IAU (' mathTxt '\approx B_y ' nmTxt moonName(1) ...
-                    mathTxt '\phi\Omega' nmTxt ', nT)'];
-                titleInfo = [moonName ' equatorial plane hodogram, ' magModelDescrip];
-            end
-        end
-        ylims = [-BcompMax, BcompMax] * 1.05;
-        xlims = ylims;
-        cfmt = {'b'};
-
-        fName = [moonName 'Hodogram' magModelDescrip];
-        PlotGeneric(xx, yy, [], windowName, titleInfo, xInfo, yInfo, fName, figDir, ...
-            figXtn, LIVE_PLOTS, figNumBase + 3, 'linear', 'linear', xlims, ylims, cfmt, 1, 1);
+        %% Hodograms
+        [xx, yy, windowName, titleInfo, xInfo, yInfo, fName, xlims, ylims] = HodogramSingle( ...
+            moonName, parentName, magModelDescrip, Bvec, SPHOUT, COMPARE_SEUFERT, COMPARE_PHIO);
+        PlotGeneric(xx, yy, [], windowName, titleInfo, xInfo, yInfo, fName, figDir, figXtn, ...
+            LIVE_PLOTS, figNumBase + 3, 'linear', 'linear', xlims, ylims, {'b'}, 1, 1);
     end
 
 end

functions/plotting/PlotBandLsq.m

@@ -130,7 +130,7 @@ function PlotBandLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, sc
     figNumBase = 3000 + 100*opt + 10*MPopt;
     windowName = [char(scName) 'Br, ' orbStr ', ' magModelDescrip];
     yy = [Br; BrSC];
-    yInfo = 'Vector component (nT)';
+    yInfo = [mathTxt 'B_r' nmTxt ' component (nT)'];
     legendStrings = [string(magModelDescrip), string(scDataName)];
     titleInfo = commonTitle;
     fName = [char(scName) parentName 'BrComparison' magModelDescrip];
@@ -139,12 +139,14 @@ function PlotBandLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, sc
 
     windowName = [char(scName) 'Bth, ' orbStr ', ' magModelDescrip];
     yy = [Bth; BthSC];
+    yInfo = [mathTxt 'B_\theta' nmTxt ' component (nT)'];
     fName = [char(scName) parentName 'BthComparison' magModelDescrip];
     PlotGeneric(xx, yy, legendStrings, windowName, titleInfo, xInfo, yInfo, fName, ...
         figDir, figXtn, LIVE_PLOTS, figNumBase + 2);
 
     windowName = [char(scName) 'Bphi, ' orbStr ', ' magModelDescrip];
     yy = [Bphi; BphiSC];
+    yInfo = [mathTxt 'B_\phi' nmTxt ' component (nT)'];
     fName = [char(scName) parentName 'BphiComparison' magModelDescrip];
     PlotGeneric(xx, yy, legendStrings, windowName, titleInfo, xInfo, yInfo, fName, ...
         figDir, figXtn, LIVE_PLOTS, figNumBase + 3);