Skip to content

Commit

Permalink
Beginning implementation of generic plotting function
Browse files Browse the repository at this point in the history
  • Loading branch information
@itsmoosh
itsmoosh committed Nov 25, 2023
1 parent df3e768 commit 4269ae9
Show file tree
Hide file tree
Showing 14 changed files with 440 additions and 129 deletions.
88 changes: 60 additions & 28 deletions PlanetMag.m
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
function [T_h, B0vec, B1vec1, B1vec2, B1vec3, outFname, header] = PlanetMag(moonName, era, ...
coordType, CALC_NEW, ALL_MODELS, DO_FFT, DO_MPAUSE, LIVE_PLOTS, specificModel, ...
specificMPmodel, outData, coeffPath, fPatternFT, nptsApprox, magPhase_deg)
coordType, CALC_NEW, ALL_MODELS, DO_FFT, DO_MPAUSE, LIVE_PLOTS, nptsApprox, specificModel, ...
specificMPmodel, outData, figDir, coeffPath, fPatternFT, fPatternTseries, figXtn, magPhase_deg)
% Evaluates planetary magnetic field for a time series at the specified moon location and inverts
% for the complex amplitudes of oscillation in that moon's frame.
%
Expand Down Expand Up @@ -37,6 +37,10 @@
% Whether to include a magnetopause screening current model.
% LIVE_PLOTS : bool, default=0
% Whether to load interactive figure windows for plots (true) or print them to disk (false).
% nptsApprox : int, default=12*365.25*3*24
% Desired number of points to use in time series for inversion. A whole number of the period of
% interest (typically synodic period, as it is the strongest oscillation) will ultimately be
% selected, which is why this number is approximate.
% specificModel : int, default=0
% Index number for the magnetospheric model to run. Options depend on the body, and setting to
% 0 selects the default model. See GetModelOpts for a description of the options.
Expand All @@ -45,14 +49,16 @@
% model See MpauseFld for a description of the options.
% outData : char, 1xF, default='out'
% Directory to use for output of complex spectrum amplitudes.
% coeffPath : char, 1xE, default='modelCoeffs'
% figDir : char, 1xG, default='figures'
% Directory to use for output figures.
% coeffPath : char, 1xH, default='modelCoeffs'
% Directory containing model coefficients files.
% fPatternFT : char, 1xG, default='FTdata'
% fPatternFT : char, 1xI, default='FTdata'
% Pattern for file names of FFT spectrum data saved to disk.
% nptsApprox : int, default=12*365.25*3*24
% Desired number of points to use in time series for inversion. A whole number of the period of
% interest (typically synodic period, as it is the strongest oscillation) will ultimately be
% selected, which is why this number is approximate.
% fPatternTseries : char, 1xJ, default='TseriesData'
% Pattern for file names of time series data saved to disk.
% figXtn : char, 1xK, default='pdf'
% Extension to use for figures, which determines the file type.
% magPhase_deg : double, default=0
% Arbitrary offset in degrees by which to rotate the magnetospheric field evaluation.
%
Expand Down Expand Up @@ -86,12 +92,15 @@
if ~exist('DO_FFT', 'var'); DO_FFT = 0; end
if ~exist('DO_MPAUSE', 'var'); DO_MPAUSE = 0; end
if ~exist('LIVE_PLOTS', 'var'); LIVE_PLOTS = 0; end
if ~exist('nptsApprox', 'var'); nptsApprox = 12*365.25*3*24; end
if ~exist('specificModel', 'var'); specificModel = 0; end
if ~exist('specificMPmodel', 'var'); specificMPmodel = 0; end
if ~exist('outData', 'var'); outData = 'out'; end
if ~exist('figDir', 'var'); figDir = 'figures'; end
if ~exist('coeffpath', 'var'); coeffPath = 'modelCoeffs'; end
if ~exist('fPatternFT', 'var'); fPatternFT = 'FTdata'; end
if ~exist('nptsApprox', 'var'); nptsApprox = 12*365.25*3*24; end
if ~exist('fPatternTseries', 'var'); fPatternTseries = 'TseriesData'; end
if ~exist('figXtn', 'var'); figXtn = 'pdf'; end
if ~exist('magPhase_deg', 'var'); magPhase_deg = 0; end

parentName = LoadSpice(moonName);
Expand Down Expand Up @@ -210,7 +219,7 @@
[Bvec, Mdip_nT, Odip_km] = MagFldParent(parentName, rM_km, thetaM, phiM, ...
MagModel, CsheetModel, magPhase_deg, SPHOUT);
end
if DO_MPAUSE
if DO_MPAUSE && ~strcmp(MPmodel, 'None')

nSW_pcc = 4 / a_AU^2 * ones(1,npts);
vSW_kms = 400 * ones(1,npts);
Expand Down Expand Up @@ -247,7 +256,7 @@
end

BD = ICAdecomposition(moonName, parentName, t_h*3600, BvecMoon, magModelDescrip, ...
SPHOUT, 1, 1, LIVE_PLOTS);
SPHOUT, 1, 1, LIVE_PLOTS, figDir, figXtn);

T_h = 1 ./ BD.f / 3600;
npeaks = length(T_h);
Expand Down Expand Up @@ -286,9 +295,9 @@
nOsc = floor(approxDur_h / Tinterest_h);
end
[TfinalFFT_h, B0vecFFT, B1vecFFT] = ExcitationSpectrum(moonName, nOsc, rate, ...
Tinterest_h, outData, fPatternFT, fPatternTseries, magPhase_deg, 0, SPHOUT, ...
DO_MPAUSE);
PlotSpectrum(moonName, LIVE_PLOTS, outData, fPatternFT);
Tinterest_h, coordType, outData, fPatternFT, fPatternTseries, ...
magPhase_deg, 0, SPHOUT, DO_MPAUSE);
PlotSpectrum(moonName, LIVE_PLOTS, outData, figDir, fPatternFT, figXtn);
end

npts = length(t_h);
Expand All @@ -311,8 +320,8 @@
Bvec2Reprod(i,:) = real(B1vec2(i) * exp(-1i * omega_ph(i) * t_h));
Bvec3Reprod(i,:) = real(B1vec3(i) * exp(-1i * omega_ph(i) * t_h));
indB = find(BeModes == sortBeModes(i));
BeStrings(i) = sprintf('|B| amp: %.4e f_Hz: %.4e T_h: %.18f', BeModes(indB), ...
1/3600/T_h(indB), T_h(indB));
BeStrings(i) = sprintf('|B| amp: %.4e f_Hz: %.4e T_h: %.18f', ...
BeModes(indB), 1/3600/T_h(indB), T_h(indB));
end
disp(BeStrings);

Expand All @@ -325,27 +334,44 @@
Bv1lbl = 'B_x'; Bv2lbl = 'B_y'; Bv3lbl = 'B_z';
end

figure; hold on;
set(gcf,'Name', ['Bvec1 full time, ' magModelDescrip]);
plot(t_h, BvecMoon(1,:));
plot(t_h, Bvec1Tot);
xlabel('Time (hr)');
ylabel([Bv1lbl ' (nT)']);
legend([Bv1lbl ' model'], [Bv1lbl ' exc']);
figure; hold on;
xx = t_h;
yy = [BvecMoon(1,:); Bvec1Tot];
windowName = ['Bvec1 full time, ' magModelDescrip];
legendStrings = [string([Bv1lbl ' model']), string([Bv1lbl ' exc'])];
titleInfo = [Bv1lbl ' reproduced vs. ' magModelDescrip];
xInfo = 'Time (hr)';
yInfo = [Bv1lbl ' (nT)'];
fName = [moonName 'Bv1_fullRepro' magModelDescrip era coordType];
PlotGeneric(xx, yy, legendStrings, windowName, titleInfo, xInfo, yInfo, figDir, ...
fName, figXtn, LIVE_PLOTS)

fig = figure('Visible', figVis); hold on;
set(gcf,'Name', ['Bvec2 full time, ' magModelDescrip]);
plot(t_h, BvecMoon(2,:));
plot(t_h, Bvec2Tot);
xlabel('Time (hr)');
ylabel([Bv2lbl ' (nT)']);
legend([Bv2lbl ' model'], [Bv2lbl ' exc']);
figure; hold on;
if ~LIVE_PLOTS
outFig = fullfile(figDir, [moonName 'Bv2_fullRepro' magModelDescrip era ...
coordType '.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end

fig = figure('Visible', figVis); hold on;
set(gcf,'Name', ['Bvec3 full time, ' magModelDescrip]);
plot(t_h, BvecMoon(3,:));
plot(t_h, Bvec3Tot);
xlabel('Time (hr)');
ylabel([Bv3lbl ' (nT)']);
legend([Bv3lbl ' model'], [Bv3lbl ' exc']);
if ~LIVE_PLOTS
outFig = fullfile(figDir, [moonName 'Bv3_fullRepro' magModelDescrip era ...
coordType '.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end

Bvec1D = BvecMoon(1,:) - Bvec1Tot;
Bvec2D = BvecMoon(2,:) - Bvec2Tot;
Expand All @@ -368,7 +394,7 @@
Bvec2D1f = 2*Bvec2D1(iTmax:end);
Bvec3D1f = 2*Bvec3D1(iTmax:end);

figure; hold on;
fig = figure('Visible', figVis); hold on;
set(gcf,'Name', ['Residual FFT for ' era ' era, ' magModelDescrip]);
plot(Tfinal_h, abs(Bvec1D1f));
plot(Tfinal_h, abs(Bvec2D1f));
Expand All @@ -380,6 +406,12 @@
xlabel('Period (hr)');
ylabel('FT\{B_i - B_{i,exc}\} (nT)');
xlim([1 Tmax]);
if ~LIVE_PLOTS
outFig = fullfile(figDir, [moonName 'ResidualFFT' magModelDescrip era coordType ...
'.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end

Bdiff = sqrt(abs(Bvec1D1f).^2 + abs(Bvec2D1f).^2 + abs(Bvec3D1f).^2);
maxBdiff = sort(Bdiff, 'descend');
Expand All @@ -399,8 +431,8 @@
'Beth_Re(nT),Beth_Im(nT),Bephi_Re(nT),Bephi_Im(nT)'];
else
BeType = 'Be1xyz_';
header = ['period(hr),B0x(nT),B0y(nT),B0z(nT),Bex_Re(nT),Bex_Im(nT),Bey_Re(nT),' ...
'Bey_Im(nT),Bez_Re(nT),Bez_Im(nT)'];
header = ['period(hr),B0x(nT),B0y(nT),B0z(nT),Bex_Re(nT),Bex_Im(nT),' ...
'Bey_Re(nT),Bey_Im(nT),Bez_Re(nT),Bez_Im(nT)'];
end
outFname = fullfile(outData, [BeType moonName '_' era '_' fEnd '.txt']);
spectrumData = [T_h' B0vec1' B0vec2' B0vec3' real(B1vec1)' imag(B1vec1)' ...
Expand Down
25 changes: 20 additions & 5 deletions comparison/CompareEarModels.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function CompareEarModels(LIVE_PLOTS, scName, xtn, SEQUENTIAL)
function CompareEarModels(LIVE_PLOTS, scName, xtn, SEQUENTIAL, coeffPath, figDir, figXtn)
% Compare magnetic field measurements from spacecraft near Earth against each implemented magnetic
% field model.
%
Expand All @@ -17,10 +17,16 @@ function CompareEarModels(LIVE_PLOTS, scName, xtn, SEQUENTIAL)
% directory must exist with this name in the ``MAG`` directory within the top-level P\lanetMag
% directory. This directory will be searched for data files with the ``xtn`` extension, and each
% of these files will be loaded.
% xtn : char, 1xD, default='.tab'
% xtn : char, 1xC, default='.tab'
% File extension for data files found in ``fullfile('MAG', sc)``, beginning with ``'.'``.
% SEQUENTIAL : bool
% Whether to plot points by index or hours relative to a reference time.
% coeffPath : char, 1xD, default='modelCoeffs'
% Directory containing model coefficients files.
% figDir : char, 1xE, default='figures'
% Directory to use for output figures.
% figXtn : char, 1xF, default='pdf'
% Extension to use for figures, which determines the file type.

% Part of the PlanetMag framework for evaluation and study of planetary magnetic fields.
% Created by Corey J. Cochrane and Marshall J. Styczinski
Expand All @@ -32,6 +38,8 @@ function CompareEarModels(LIVE_PLOTS, scName, xtn, SEQUENTIAL)
if ~exist('scName', 'var'); scName = "Swarm"; end
if ~exist('xtn', 'var'); xtn = '.tab'; end
if ~exist('SEQUENTIAL', 'var'); SEQUENTIAL = 1; end
if ~exist('figDir', 'var'); figDir = 'figures'; end
if ~exist('figXtn', 'var'); figXtn = 'pdf'; end

parentName = 'Earth';
sc = char(scName);
Expand Down Expand Up @@ -83,22 +91,29 @@ function CompareEarModels(LIVE_PLOTS, scName, xtn, SEQUENTIAL)
xDescrip = 'Time relative to NY 2020 (h)';
end
lat_deg = 90 - rad2deg(theta);
figure; hold on;

fig = figure('Visible', figVis); hold on;
set(gcf,'Name', [char(scName) ' latitudes']);
plot(xx, lat_deg, 'DisplayName', 'latitude');
xlabel(xDescrip);
ylabel('Latitude (degrees)');
title('Swarm ABC latitudes on 2020-01-01');
legend();
if ~LIVE_PLOTS
outFig = fullfile(figDir, ['SwarmABClatitudes.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end

%% Plot and calculate products
nOpts = 1; nMPopts = 0;
opts = 1:nOpts;
MPopts = -1:-1;
for opt=opts
for MPopt=MPopts
GetBplotAndLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, scName, ...
parentName, spkParent, orbStr, opt, MPopt, SEQUENTIAL);
PlotBandLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, scName, ...
parentName, spkParent, orbStr, opt, MPopt, SEQUENTIAL, coeffPath, figDir, ...
figXtn, LIVE_PLOTS);
end
end

Expand Down
19 changes: 14 additions & 5 deletions comparison/CompareJupModels.m
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
function CompareJupModels(LIVE_PLOTS, scName, moonName, orbNum, moonProx_RP, PlanetMinDist_RP, ...
PlanetMaxDist_RP, dataDir, SEQUENTIAL, FULLORBITS, FLYBYS, JUNOTOO)
PlanetMaxDist_RP, dataDir, figDir, figXtn, SEQUENTIAL, coeffPath, FULLORBITS, FLYBYS, JUNOTOO)
% Compare magnetic field measurements from spacecraft near Jupiter against each implemented
% magnetic field model.
%
Expand Down Expand Up @@ -42,8 +42,14 @@ function CompareJupModels(LIVE_PLOTS, scName, moonName, orbNum, moonProx_RP, Pla
% the parent planet's radius.
% dataDir : char, 1xD, default='out'
% Name of directory where excitation moments are printed to disk.
% figDir : char, 1xE, default='figures'
% Directory to use for output figures.
% figXtn : char, 1xF, default='pdf'
% Extension to use for figures, which determines the file type.
% SEQUENTIAL : bool, default=1
% Whether to plot points by index or hours relative to a reference time (J2000).
% coeffPath : char, 1xG, default='modelCoeffs'
% Directory containing model coefficients files.
% FULLORBITS : bool, default=1
% Whether to evalate goodness of fit from full-orbit data in System III coordinates beyond some
% threshold distance ``rMinMoon_RP``. Independent of ``FLYBYS``.
Expand All @@ -69,6 +75,8 @@ function CompareJupModels(LIVE_PLOTS, scName, moonName, orbNum, moonProx_RP, Pla
if ~exist('PlanetMinDist_RP', 'var'); PlanetMinDist_RP = 10; end
if ~exist('PlanetMaxDist_RP', 'var'); PlanetMaxDist_RP = 20; end
if ~exist('dataDir', 'var'); dataDir = 'out'; end
if ~exist('figDir', 'var'); figDir = 'figures'; end
if ~exist('figXtn', 'var'); figXtn = 'pdf'; end
if ~exist('SEQUENTIAL', 'var'); SEQUENTIAL = 1; end
if ~exist('FULLORBITS', 'var'); FULLORBITS = 1; end
if ~exist('FLYBYS', 'var'); FLYBYS = 0; end
Expand Down Expand Up @@ -250,12 +258,13 @@ function CompareJupModels(LIVE_PLOTS, scName, moonName, orbNum, moonProx_RP, Pla
for opt=opts
for MPopt=MPopts
if FULLORBITS
GetBplotAndLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, scName, ...
parentName, spkParent, orbStr, opt, MPopt, SEQUENTIAL);
PlotBandLsq(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, scName, ...
parentName, spkParent, orbStr, opt, MPopt, SEQUENTIAL, coeffPath, figDir, ...
figXtn, LIVE_PLOTS);
else
GetBplotAndLsqMoon(fbets, fbt_h, fbr_km, fbtheta, fbphi, fbxyz_km, r_RM, BxSC, ...
PlotBandLsqMoon(fbets, fbt_h, fbr_km, fbtheta, fbphi, fbxyz_km, r_RM, BxSC, ...
BySC, BzSC, scName, parentName, spkParent, moonName, char(scName), fbStr, ...
opt, MPopt, SEQUENTIAL, dataDir, jt_h);
opt, MPopt, SEQUENTIAL, dataDir, figDir, figXtn, LIVE_PLOTS, jt_h);
end
end
end
Expand Down
35 changes: 27 additions & 8 deletions comparison/CompareNepModels.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL, coeffPath, figDir, figXtn)
% Compare magnetic field measurements from spacecraft near Neptune against each implemented
% magnetic field model.
%
Expand All @@ -18,6 +18,12 @@ function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
% SEQUENTIAL : bool, default=0
% Whether to plot points by index or hours relative to a reference time (typically closest
% approach).
% coeffPath : char, 1xC, default='modelCoeffs'
% Directory containing model coefficients files.
% figDir : char, 1xD, default='figures'
% Directory to use for output figures.
% figXtn : char, 1xE, default='pdf'
% Extension to use for figures, which determines the file type.

% Part of the PlanetMag framework for evaluation and study of planetary magnetic fields.
% Created by Corey J. Cochrane and Marshall J. Styczinski
Expand All @@ -28,6 +34,8 @@ function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
if ~exist('LIVE_PLOTS', 'var'); LIVE_PLOTS = 0; end
if ~exist('scName', 'var'); scName = "Voyager 2"; end
if ~exist('SEQUENTIAL', 'var'); SEQUENTIAL = 0; end
if ~exist('figDir', 'var'); figDir = 'figures'; end
if ~exist('figXtn', 'var'); figXtn = 'pdf'; end

cspice_kclear;
parentName = 'Neptune';
Expand Down Expand Up @@ -93,7 +101,8 @@ function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
x = r .* sin(th) .* cos(ph);
y = r .* sin(th) .* sin(ph);
z = r .* cos(th);
figure; hold on

fig = figure('Visible', figVis); hold on
dx = x - xyz_km(1,:);
dy = y - xyz_km(2,:);
dz = z - xyz_km(3,:);
Expand All @@ -106,25 +115,30 @@ function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
ylabel('Coordinate diff (km)');
title('Location difference in NLS frame, PDS - SPICE');
legend()
if ~LIVE_PLOTS
outFig = fullfile(figDir, ['Voyager2NeptuneFlybyPDSvsSPICE.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end

r_km = r; theta = th; phi = ph; xyz_km = [x; y; z];
end


%% Plot and calculate products
nOpts = 1; nMPopts = 0;
opts = 1:nOpts;
%MPopts = 1:(nMPopts + 1); % Add 1 to force noMP model in addition
MPopts = -1:-1;
for opt=opts
for MPopt=MPopts
GetBplotAndLsqNeptune(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, ...
scName, spkParent, orbStr, opt, MPopt, SEQUENTIAL, 1, 1, 1, 1);
PlotBandLsqNeptune(ets, t_h, r_km, theta, phi, xyz_km, BrSC, BthSC, BphiSC, ...
scName, spkParent, orbStr, opt, MPopt, SEQUENTIAL, coeffPath, figDir, figXtn, ...
LIVE_PLOTS, 1, 1, 1, 1);
end
end


%% Plot trajectory
figure(1000); clf(); hold on;
fig = figure(1000, 'Visible', figVis); clf(); hold on;
set(gcf,'Name', 'Trajectories');
the = linspace(0,pi,19); ph = linspace(0,2*pi,37);
[the2D, ph2D] = meshgrid(the,ph);
Expand Down Expand Up @@ -177,5 +191,10 @@ function CompareNepModels(LIVE_PLOTS, scName, SEQUENTIAL)
name{3} = "SPICE in NLS as defined in O8";

legend([scTraj{1} scTraj{2} scTraj{3}], [name{1}, name{2}, name{3}])


if ~LIVE_PLOTS
outFig = fullfile(figDir, ['Voyager2NeptuneFlybyTrajectories.' figXtn]);
saveas(fig, outFig)
disp(['Figure saved to ' outFig '.'])
end
end
Loading

0 comments on commit 4269ae9

Please sign in to comment.