1+ import logging
2+ import re
3+ import numpy as np
4+ from pyspedas .tplot_tools import get_data
5+
6+ logging .captureWarnings (True )
7+ logging .basicConfig (format = '%(asctime)s: %(message)s' , datefmt = '%d-%b-%y %H:%M:%S' , level = logging .INFO )
8+
9+ def cluster_get_codif_dist (tname , probe = None ):
10+ """
11+ Returns 3D particle data structures containing Cluster CODIF
12+ data for use with pySPEDAS particle routines.
13+
14+ Input
15+ ----------
16+ tname: str
17+ tplot variable name containing the CODIF distribution data
18+ Example: 'ions_3d__C1_CP_CIS_CODIF_HS_H1_PSD'
19+
20+ Parameters
21+ ----------
22+ probe: str
23+ Spacecraft probe #
24+
25+ Returns
26+ ----------
27+ 3D particle data structure(s) containing Cluster CODIF distribution functions
28+ """
29+ data = get_data (tname )
30+
31+ if data is None :
32+ logging .error ('Problem extracting the CODIF distribution data.' )
33+ return
34+
35+ if any (data [i ].ndim == 2 for i in (2 , 3 , 4 )):
36+ raise ValueError (
37+ "data[2], data[3], and data[4] must all be 1-D arrays"
38+ )
39+
40+ # check which species by looking at second last part of input string
41+ # this could be a function
42+ parts = tname .split ("_" )
43+ if re .fullmatch (r"[A-Z][a-z]?\d" , parts [- 2 ]):
44+ species = parts [- 2 ].lower ()
45+ else :
46+ raise ValueError (f"No valid species found in: { tname } )
47+
48+ # possible masses in amu
49+ if species == 'h1' :
50+ mass = 1.04535e-2
51+ charge = 1.0
52+ elif species == 'he1' :
53+ mass = 4.18138e-2
54+ charge = 1.0
55+ elif species == 'o1' :
56+ mass = 0.167255
57+ charge = 1.0
58+ else :
59+ logging .error ('Invalid species: ' + species )
60+ #return
61+
62+ out = {'project_name' : 'Cluster' ,
63+ 'spacecraft' : probe ,
64+ 'species' : species ,
65+ 'data_name' : 'CODIF ' + species ,
66+ 'charge' : charge ,
67+ 'units_name' : 'df_km' ,
68+ 'mass' : mass }
69+
70+ # --- Reformat data
71+ # Shuffle the output to be [time, energy, phi, theta]
72+ # (data[1] is currently [time, theta, phi, energy])
73+ out_data = data [1 ].transpose ([0 , 3 , 2 , 1 ]) #* 1.0E-18 * 1.0E-12 # convert from s^3 km^-6 to s^3 cm^-6
74+ out_bins = np .zeros (out_data .shape ) + 1
75+ out_dphi = np .zeros (out_data .shape ) + 22.5 # 16 pixels in 360 degree window
76+ out_dtheta = np .zeros (out_data .shape ) + 22.5 # 8 pixels in 180 degree window
77+ out_denergy = np .zeros (out_data .shape )
78+
79+ # theta
80+ # elevations are constant across time
81+ theta_reform = np .asarray (data [2 ]) # shape: (theta,)
82+ theta_reform = theta_reform [None , None , :] # shape: (1, 1, theta)
83+ theta_rebinned = np .broadcast_to (
84+ theta_reform ,
85+ (len (data [4 ]), len (data [3 ]), theta_reform .shape [2 ])
86+ ) # shape: (energy, phi, theta)
87+ out_theta = np .broadcast_to (
88+ theta_rebinned ,
89+ (len (data [0 ]),) + theta_rebinned .shape
90+ ) # shape: (time, energy, phi, theta)
91+
92+ # energy
93+ energy_len = len (data [4 ])
94+ energy_reform = np .asarray (data [4 ])[:, None , None ] # shape: (energy, 1, 1)
95+ energy_table = np .broadcast_to (
96+ energy_reform ,
97+ (len (data [4 ]), len (data [3 ]), len (data [2 ]))
98+ ) # shape: (energy, phi, theta)
99+ out_energy = np .broadcast_to (
100+ energy_table ,
101+ (len (data [0 ]),) + energy_table .shape
102+ ) # shape: (time, energy, phi, theta)
103+
104+ # phi
105+ phi_reform = np .asarray (data [3 ])[None , :, None ] # shape: (1, phi, 1)
106+ phi_rebinned = np .broadcast_to (
107+ phi_reform ,
108+ (len (data [4 ]), len (data [3 ]), len (data [2 ]))
109+ ) # shape: (energy, phi, theta)
110+ out_phi = np .broadcast_to (
111+ phi_rebinned ,
112+ (len (data [0 ]),) + phi_rebinned .shape
113+ ) # shape: (time, energy, phi, theta)
114+
115+ out_list = []
116+ for time_idx , time in enumerate (data [0 ]):
117+ out_table = {** out }
118+ out_table ['data' ] = out_data [time_idx , :]
119+ out_table ['bins' ] = out_bins [time_idx , :]
120+ out_table ['theta' ] = out_theta [time_idx , :]
121+ out_table ['phi' ] = out_phi [time_idx , :]
122+ out_table ['energy' ] = out_energy [time_idx , :]
123+ out_table ['dtheta' ] = out_dtheta [time_idx , :]
124+ out_table ['dphi' ] = out_dphi [time_idx , :]
125+ out_table ['denergy' ] = out_denergy [time_idx , :]
126+ out_table ['n_energy' ] = energy_len
127+ out_table ['n_theta' ] = len (data [3 ])
128+ out_table ['n_phi' ] = len (data [2 ])
129+ out_table ['start_time' ] = time # note: assumes the data weren't centered
130+ out_table ['end_time' ] = time + 4 # one spin = 4 seconds
131+ out_list .append (out_table )
132+
133+ return out_list
0 commit comments