diff --git a/PSPoperations/__pycache__/data_handling.cpython-38.pyc b/PSPoperations/__pycache__/data_handling.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..49509b62a75d10286c70b8515c840f6b0aa0494a
Binary files /dev/null and b/PSPoperations/__pycache__/data_handling.cpython-38.pyc differ
diff --git a/PSPoperations/__pycache__/data_quality.cpython-38.pyc b/PSPoperations/__pycache__/data_quality.cpython-38.pyc
index a318e2cbeadf6a52cd370cd7bc005d2f437c39e9..66ad3157df0b1fe28a53a232e9dd6b66ac18abab 100644
Binary files a/PSPoperations/__pycache__/data_quality.cpython-38.pyc and b/PSPoperations/__pycache__/data_quality.cpython-38.pyc differ
diff --git a/PSPoperations/__pycache__/data_transformation.cpython-38.pyc b/PSPoperations/__pycache__/data_transformation.cpython-38.pyc
index 6fc03a2c8fc606a3e9f7fc3c9073b61f16587210..8991004c53960278dba9ff7e59bf24cb70a95a15 100644
Binary files a/PSPoperations/__pycache__/data_transformation.cpython-38.pyc and b/PSPoperations/__pycache__/data_transformation.cpython-38.pyc differ
diff --git a/PSPoperations/__pycache__/plot_settings.cpython-38.pyc b/PSPoperations/__pycache__/plot_settings.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..2e9945b7d6381f47e2d8e13affba059927fd0854
Binary files /dev/null and b/PSPoperations/__pycache__/plot_settings.cpython-38.pyc differ
diff --git a/PSPoperations/data_handling.py b/PSPoperations/data_handling.py
new file mode 100644
index 0000000000000000000000000000000000000000..3772ca29c3fe956936c625488ba57f8e30934eb7
--- /dev/null
+++ b/PSPoperations/data_handling.py
@@ -0,0 +1,73 @@
+import os
+import typing as tp
+import numpy as np
+
+# CDF library (is needed to interface with the measurement data files)
+# see https://cdf.gsfc.nasa.gov/
+os.environ["CDF_LIB"] = "/usr/local/cdf/lib"
+
+
+def cdf_slice(cdf_file, key: str):
+ """
+ DOC
+
+ :param cdf_file:
+ :param key:
+ :return:
+ """
+ return cdf_file[key][...]
+
+
+def data_generation(cdf_file) -> tp.Dict:
+ """
+ DOC
+
+ :param cdf_file:
+ :return:
+ """
+ data = {
+ "dqf": cdf_slice(cdf_file, key="general_flag"),
+ "epoch": cdf_slice(cdf_file, key="Epoch"),
+ "pos": cdf_slice(cdf_file, key="sc_pos_HCI"),
+ "vr": cdf_slice(cdf_file, key="vp_moment_RTN")[:, 0],
+ "dvrhi": cdf_slice(cdf_file, key="vp_moment_RTN_deltahigh")[:, 0],
+ "dvrlo": cdf_slice(cdf_file, key="vp_moment_RTN_deltalow")[:, 0],
+ "wp": cdf_slice(cdf_file, key="wp_moment"),
+ "dwphi": cdf_slice(cdf_file, key="wp_moment_deltahigh"),
+ "dwplo": cdf_slice(cdf_file, key="wp_moment_deltalow")
+ }
+
+ return data
+
+
+def array_reduction(data_array: np.ndarray,
+ index_list: np.ndarray) -> np.ndarray:
+ """
+ Reduction of an array by a list of indices.
+
+ :param data_array: NDARRAY,
+ Array that is to be reduced
+ :param index_list: NDARRAY,
+ List of indices that should be deleted
+ :return: NDARRAY,
+ The reduced array
+ """
+ # Choose axis=0 to delete rows and not flatten the array
+ reduced_array = np.delete(data_array, index_list, axis=0)
+
+ return reduced_array
+
+
+def full_reduction(data_dict: tp.Dict,
+ bad_indices: np.array) -> tp.Dict:
+ """
+ DOC
+
+ :param data_dict:
+ :param bad_indices:
+ :return:
+ """
+ for key in data_dict.keys():
+ data_dict[key] = array_reduction(data_dict[key], bad_indices)
+
+ return data_dict
diff --git a/PSPoperations/data_quality.py b/PSPoperations/data_quality.py
index 416702ba192b9d04aaf31b84ecb11abdbedd7a54..bdb6be47e210f365487e76aa9a3f6cf1522d5806 100644
--- a/PSPoperations/data_quality.py
+++ b/PSPoperations/data_quality.py
@@ -1,4 +1,5 @@
import numpy as np
+import typing as tp
def general_flag(general_flag_array: np.ndarray) -> np.ndarray:
@@ -22,19 +23,34 @@ def general_flag(general_flag_array: np.ndarray) -> np.ndarray:
return index_list
-def array_reduction(data_array: np.ndarray,
- index_list: np.ndarray) -> np.ndarray:
+def meas_failed(meas_array: np.ndarray) -> np.ndarray:
"""
- Reduction of an array by a list of indices.
-
- :param data_array: NDARRAY,
- Array that is to be reduced
- :param index_list: NDARRAY,
- List of indices that should be deleted
+ Takes a 1D array and returns indices of failed measurements (as
+ indicated by a value of -1e30, SWEAP documentation)
+
+ :param meas_array: NDARRAY,
+ Measurement array
:return: NDARRAY,
- The reduced array
+ Index array
"""
- # Choose axis=0 to delete rows and not flatten the array
- reduced_array = np.delete(data_array, index_list, axis=0)
-
- return reduced_array
+ index_array = np.where(meas_array <= -0.5e30)[0]
+
+ return index_array
+
+
+def full_meas_eval(data_dict: tp.Dict):
+ """DOC"""
+ col_ind = []
+
+ for key in data_dict.keys():
+
+ # Skip unnecessary keys (maybe make this more dynamic?)
+ if key in ["pos", "dqf", "epoch"]:
+ continue
+
+ col_ind.append(meas_failed(data_dict[key]))
+
+ # Make an array of unique indices
+ col_ind = np.unique(np.concatenate(col_ind))
+
+ return col_ind
diff --git a/PSPoperations/data_transformation.py b/PSPoperations/data_transformation.py
index 776c806717e0953fb7634ead7962be722b59f21a..a3df0292cb5f080f0d595844cfdfd2644c382d08 100644
--- a/PSPoperations/data_transformation.py
+++ b/PSPoperations/data_transformation.py
@@ -1,5 +1,7 @@
import numpy as np
import typing as tp
+# The warnings (for me) here don't seem to matter
+from astropy.constants import k_B, m_p
def pos_cart_to_sph(position_matrix: np.ndarray) -> tp.Tuple:
@@ -28,30 +30,15 @@ def pos_cart_to_sph(position_matrix: np.ndarray) -> tp.Tuple:
return r, theta, phi
-def vel_cart_to_sph(vel_matrix, r, theta, phi):
+def wp_to_temp(thermal_speed: np.ndarray) -> np.ndarray:
"""
- ???
+ Calculate temperature from thermal speed of protons.
- :param vel_matrix:
- :param r:
- :param theta:
- :param phi:
- :return:
+ :param thermal_speed: NDARRAY,
+ Thermal speed (wp) of protons
+ :return: NDARRAY,
+ Temperature from wp = sqrt(2k_BT / m)
"""
- # Sanity checks
- assert vel_matrix.shape[1] == 3, \
- "PLEASE CHECK POSITION MATRIX FOR COORDINATE CONVERSION!"
- assert r.shape[0] == theta.shape[0] == phi.shape[0] == vel_matrix.shape[0], \
- "PLEASE CHECK DIMENSIONS OF R, THETA, PHI AND VEL_MATRIX!"
-
- vx = vel_matrix[:, 0]
- vy = vel_matrix[:, 1]
- vz = vel_matrix[:, 2]
-
- vr = np.sin(theta) * np.cos(phi) * vx + \
- np.sin(theta) * np.sin(phi) * vy + \
- np.cos(theta) * vz
-
- return vr
+ wp_si = thermal_speed * 1e3
-
\ No newline at end of file
+ return wp_si ** 2 * m_p.value / (2 * k_B.value)
diff --git a/PSPoperations/plot_settings.py b/PSPoperations/plot_settings.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a50ac8149629e0f3a58cfc0687d1861b336adf9
--- /dev/null
+++ b/PSPoperations/plot_settings.py
@@ -0,0 +1,63 @@
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+
+def rc_setup():
+ """Generalized plot attributes"""
+ mpl.rcParams["xtick.direction"] = "in"
+ mpl.rcParams["ytick.direction"] = "in"
+ mpl.rcParams["xtick.minor.visible"] = "True"
+ mpl.rcParams["ytick.minor.visible"] = "True"
+
+
+def plot_r_vr(r, vr, v_ulim, v_lolim, save_ind: str = "no"):
+ """
+ Plots radial velocity wrt heliocentric distance
+
+ :param r: NDARRAY,
+ Distance (must be in km)
+ :param vr: NDARRAY,
+ Radial velocity (must be in km/s)
+ :param v_ulim, v_lolim: NDARRAY,
+ Upper (lower) uncertainty of vr
+ :param save_ind: STR,
+ "yes" if plot should be saved (default: no)
+ """
+ fig, ax = plt.subplots(figsize=(10, 6))
+ plt.plot(r / 6.957e5, vr)
+ plt.fill_between(r / 6.957e5, vr + v_ulim, vr - v_lolim, color="grey",
+ alpha=0.25)
+ ax.set(
+ xlabel="r [R$_\\odot$]", ylabel="v$_r$ [km s$^{-1}$]",
+ ylim=(0, 800)
+ )
+ ax.grid(True, alpha=0.5)
+
+ if save_ind.lower() == "yes":
+ plt.savefig("plot_vr.png")
+
+
+def plot_r_temp(r, T, T_ulim, T_lolim, save_ind: str = "no"):
+ """
+ Plots logarithmic temperature wrt heliocentric distance
+
+ :param r: NDARRAY,
+ Heliocentric distance (mst be in km)
+ :param T: NDARRAY,
+ logarithmic Temperature (Kelvin)
+ :param T_ulim, T_lolim: NDARRAY,
+ Upper (lower) uncertainty of logT
+ :param save_ind: STR,
+ "yes" if plot should be saved (default: no)
+ """
+ fig, ax = plt.subplots(figsize=(10, 6))
+ plt.plot(r / 6.957e5, T)
+ plt.fill_between(r / 6.957e5, T_ulim, T - T_lolim,
+ color="grey", alpha=0.25)
+ ax.set(
+ xlabel="r [R$_\\odot$]", ylabel="T [K]", yscale="log"
+ )
+ ax.grid(True, alpha=0.5)
+
+ if save_ind.lower() == "yes":
+ plt.savefig("plot_logT.png")
diff --git a/TESTING/plot_logT.png b/TESTING/plot_logT.png
new file mode 100644
index 0000000000000000000000000000000000000000..a33aca3b4dcee0013851f48a9a29c71eedb6d0f9
Binary files /dev/null and b/TESTING/plot_logT.png differ
diff --git a/TESTING/plot_vr.png b/TESTING/plot_vr.png
new file mode 100644
index 0000000000000000000000000000000000000000..004d704acaa00458951b0ab61d35407892117aa0
Binary files /dev/null and b/TESTING/plot_vr.png differ
diff --git a/TESTING/psp_swp_spc_l3i_20211117_v02.cdf b/TESTING/psp_swp_spc_l3i_20211117_v02.cdf
new file mode 100644
index 0000000000000000000000000000000000000000..adb641f3f4584932c96daa492fb563b0579deacc
Binary files /dev/null and b/TESTING/psp_swp_spc_l3i_20211117_v02.cdf differ
diff --git a/TESTING/testing.py b/TESTING/testing.py
new file mode 100644
index 0000000000000000000000000000000000000000..00c593eaee9097d1d4de98808947f6dbfead8129
--- /dev/null
+++ b/TESTING/testing.py
@@ -0,0 +1,49 @@
+import os
+import sys
+from unipath import Path
+import numpy as np
+import matplotlib.pyplot as plt
+from spacepy import pycdf
+
+# Append custom Python modules from parent directory
+sys.path.append(Path(sys.path[0]).parent)
+from PSPoperations import data_quality as dq
+from PSPoperations import data_transformation as dt
+from PSPoperations import plot_settings as ps
+from PSPoperations import data_handling as dh
+
+
+# Plot set up
+ps.rc_setup()
+
+# CDF library (is needed to interface with the measurement data files)
+# see https://cdf.gsfc.nasa.gov/
+os.environ["CDF_LIB"] = "/usr/local/cdf/lib"
+
+# open CDF file and generate faulty index array
+cdf_data = pycdf.CDF("psp_swp_spc_l3i_20211117_v02.cdf")
+data = dh.data_generation(cdf_data)
+
+# Indices of non-usable data from general flag + reduction
+bad_ind = dq.general_flag(data["dqf"])
+data = dh.full_reduction(data, bad_ind)
+
+# Additional reduction from failed measurement indices
+mf_ind = dq.full_meas_eval(data)
+data = dh.full_reduction(data, mf_ind)
+
+# Add to data dictionary (spherical coordinates)
+data["r"], data["theta"], data["phi"] = dt.pos_cart_to_sph(data["pos"])
+
+# Add to data dictionary (temperature)
+data["T"] = dt.wp_to_temp(data["wp"])
+data["dThi"] = dt.wp_to_temp(data["dwphi"])
+data["dTlo"] = dt.wp_to_temp(data["dwplo"])
+
+# vr plot
+ps.plot_r_vr(data["r"], data["vr"], data["dvrhi"], data["dvrlo"],
+ save_ind="yes")
+ps.plot_r_temp(data["r"], data["T"], data["dThi"], data["dTlo"],
+ save_ind="yes")
+
+plt.show()
diff --git a/main.py b/main.py
index 075b2b79b578f7023a4c17e302d3c79813c8bb1e..469a9cf193c86178748f61aa57913b1099e04b15 100644
--- a/main.py
+++ b/main.py
@@ -1,3 +1,20 @@
+import os
+import sys
+import matplotlib.pyplot as plt
+import numpy as np
+from spacepy import pycdf
+from PSPoperations import data_quality as dq
+from PSPoperations import data_transformation as dt
+from PSPoperations import plot_settings as ps
+
+# Plot parameters
+ps.rc_setup()
+
+# CDF library (is needed to interface with the measurement data files)
+# see https://cdf.gsfc.nasa.gov/
+os.environ["CDF_LIB"] = "/usr/local/cdf/lib"
+
+
"""
PSEUDOCODE
@@ -18,8 +35,71 @@ BACKGROUND
OPERATIONS
- Sort through "general flag" (only use where set to 0)
+ -- Also sort through each array to find entries with -1e30, which
+ marks failed measurements
- Calculate spherical heliocentric coordinates (HIC)
- Transform to usable data parameters (vr, log_T, log_rho)
+- Generate log-file with important information
+ -- Total number of data points
+ -- Reduced number of data points
+ -- Epoch from start to finish
"""
+ENCOUNTER_NUM = "encounter_5"
+DATA_LOCATION = sys.path[0]+"/DATA/"+ENCOUNTER_NUM
+
+# SANITY CHECK IF DATA LOCATION EXISTS
+if not os.path.isdir(DATA_LOCATION):
+ print(f"\n{DATA_LOCATION} IS NOT A VALID DIRECTORY!\n")
+ sys.exit(0)
+
+
+def main():
+ # Total array setup
+ r = theta = phi = np.array([])
+ vr = dvhi = dvlo = np.array([])
+
+ for file in os.listdir(DATA_LOCATION):
+ cdf_data = pycdf.CDF(DATA_LOCATION+"/"+file)
+
+ dqf = cdf_data["general_flag"][...]
+ bad_ind = dq.general_flag(dqf)
+
+ # generate reduced epoch array
+ time_file = dq.array_reduction(cdf_data["Epoch"][...], bad_ind)
+
+ # Position
+ pos_file = dq.array_reduction(cdf_data["sc_pos_HCI"][...], bad_ind)
+ r_file, theta_file, phi_file = dt.pos_cart_to_sph(pos_file)
+
+ # Velocity
+ vr_file = dq.array_reduction(
+ cdf_data["vp_moment_RTN"][...][:, 0],
+ bad_ind)
+ dvhi_file = dq.array_reduction(
+ cdf_data["vp_moment_RTN_deltalow"][...][:, 0],
+ bad_ind)
+ dvlo_file = dq.array_reduction(
+ cdf_data["vp_moment_RTN_deltahigh"][...][:, 0],
+ bad_ind)
+
+ # Update total arrays
+ r = np.append(r, r_file)
+ vr = np.append(vr, vr_file)
+
+ plt.figure()
+ plt.plot(r, vr)
+ plt.ylim(0, 1000)
+
+ plt.show()
+
+
+def plot_setup():
+ fig_vr, ax_vr = plt.subplots()
+ ax_vr.set(
+ xlabel="r [R$_\\odot$]", ylabel="v$_r$ [km/s$^{-1}$]",
+ ylim=(0, 800)
+ )
+
-DATA_LOCATION =
+if __name__ == "__main__":
+ main()