gradient.py
7.15 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
#! /bin/env python
# coding: utf8
import numpy as np
import mms_tools_lib as mms ## Should this be modified to "import read_data ?"
import os
import matplotlib.pyplot as plt
class calcGradient(object):
"""
Object which purpose is to facilitate the use of the gradient object, in particular for MMS data
NOTE: The file / data management could be merged with the differential timing method
Use the donnee object to provide the interpolated data required to compute the gradient
Raise a warning if only one data point is provided before the interpolation
"""
def __init__(self, orbit_files = None, data_files = None, orbit_vars = None, data_var = None, test = False ):
#self.datapath = '/home/vernisse/Etudes_MMS/statistique/All_CS_MultiSat/time_table_20150908_101859_20150908_101933'
self.info = {}
if test : self._initDefaultvars() #For Testing only
if orbit_files is not None : self.info['orbit_files'] = orbit_files
if data_files is not None : self.info['data_files'] = data_files
if orbit_vars is not None : self.info['orbit_vars'] = orbit_vars
if data_var is not None : self.info['data_var'] = data_var
self._initializeVars()
def _initializeVars(self, dictionary_ = {'rx' : None, 'ry' : None, 'rz' : None}):
"""
Input what are the x, y, z positions of the spacecraft
for mms it is usually mec_rx_gse
if no dictionnary is input, will look for rx, ry, rz, itself
"""
self.list_ = {}
## Shall be rewriten using regular expression
if any([i is None for i in dictionary_.values()]):
print 'researching position variables'
for var in self.info["orbit_vars"] :
if 'rx' in var : self.list_['rx'] = var
elif 'ry' in var : self.list_['ry'] = var
elif 'rz' in var : self.list_['rz'] = var
print 'Match found : ', [i +' :'+ self.list_[i] for i in ['rx', 'ry', 'rz']]
else :
self.list.update(dictionary_)
def _interpolateAll(self):
"""
Interpolate all rx, ry, rz and data for all satellites on satellite 1 time axis
"""
self.orbitData = []
for file_ in self.info['orbit_files'] : self.orbitData.append(mms.donnee(file_))
self.data = []
for file_ in self.info['data_files'] : self.data.append(mms.donnee(file_))
listSats = ['Sat1', 'Sat2', 'Sat3', 'Sat4']
listPositions = ['rx', 'ry', 'rz']
for i, sat in enumerate(listSats) :
for composante in listPositions:
self.data[0].addField(self.orbitData[i].get_data(self.list_[composante]), self.orbitData[i].get_time(), sat + '_' + composante)
self.data[0].addField(self.data[i].get_data(self.info['data_var']), self.data[i].get_time(), sat + '_' + self.info['data_var'])
def _initDefaultvars(self):
"""
Initialize default values (for testing)
"""
DefaultOrbitFiles = ['obj7.cef', 'obj8.cef', 'obj9.cef', 'obj10.cef']
#DefaultDataFiles = ['obj3.cef', 'obj4.cef', 'obj5.cef', 'obj6.cef']
DefaultDataFiles = ['obj15.cef', 'obj16.cef', 'obj17.cef', 'obj18.cef']
DefaultPath = '/home/vernisse/Etudes_MMS/statistique/AllDataBrst_MultiSat/time_table_20150908_100000_20150908_113000'
self.info['orbit_files'] = [os.path.join(DefaultPath, file_) for file_ in DefaultOrbitFiles]
self.info['data_files'] = [os.path.join(DefaultPath, file_) for file_ in DefaultDataFiles]
self.info['orbit_vars'] = ['mec_rx_gsm', 'mec_ry_gsm', 'mec_rz_gsm']
self.info['data_var'] = 'dis_numberdensity_brst'
def evaluateGradient(self):
self._interpolateAll()
self.gradient = Gradient(self.data[0].donnee, datakey = self.info['data_var'])
self.gradient.getGradient()
def plot(self):
plt.plot_date(self.data[0].get_time(), self.gradient.Data['gradient_'+self.info['data_var']]['x'], 'red')
plt.plot_date(self.data[0].get_time(), self.gradient.Data['gradient_'+self.info['data_var']]['y'], 'green')
plt.plot_date(self.data[0].get_time(), self.gradient.Data['gradient_'+self.info['data_var']]['z'], 'blue')
def show(self):
plt.show()
def tocef(self, filename = 'test.cef'):
gx = 'gradient_'+self.info['data_var']+'_x'
gy = 'gradient_'+self.info['data_var']+'_y'
gz = 'gradient_'+self.info['data_var']+'_z'
mms.cefwrite({'epoch' : self.data[0].get_time(), gx : self.gradient.Data['gradient_'+self.info['data_var']]['x'], gy : self.gradient.Data['gradient_'+self.info['data_var']]['y'], gz : self.gradient.Data['gradient_'+self.info['data_var']]['z']}, filename = filename)
class Gradient(object) :
"""
Gradient calculation method
Require 16 inputs : x,y,z positions of 4 satellites
and the scalar field on 4 satellites
Return the gradient of the scalar field for each data point
All input fields MUST BE already interpolated to the same time axis
"""
def __init__(self, Data = {'Sat1_rx' : None, 'Sat1_ry' : None, 'Sat1_rz' : None, 'Sat2_rx' : None, 'Sat2_ry' : None, 'Sat2_rz' : None, 'Sat3_rx' : None, 'Sat3_ry' : None, 'Sat3_rz' : None, 'Sat4_rx' : None, 'Sat4_ry' : None, 'Sat4_rz' : None, 'Sat1_data' : None, 'Sat2_data' : None, 'Sat3_data' : None, 'Sat4_data' : None}, datakey = None):
self.datakey = datakey
self._checkKeys(Data)
self.Data = {}
self.setData(Data)
#self.getGradient()
def _checkKeys(self, dictionary_):
self.satkeys = ['Sat1', 'Sat2','Sat3', 'Sat4']
self.poskeys = ['rx','ry','rz']
allkeys = []
for sat in self.satkeys :
for pos in self.poskeys : allkeys.append(sat + '_' + pos)
allkeys.append(sat + '_' + self.datakey)
if all([i in dictionary_ for i in allkeys]) : print 'Keys checked'
else : print 'WARNING, MISSING DATA TO COMPUTE GRADIENT'
def setData(self, dictionary_):
"""
Two methods :
Update Gradient dictionary_ data.
"""
self.Data.update(dictionary_)
self.AbscSize = len(self.Data['Sat1_'+self.datakey])
def getMesoCentre(self):
#1 : construction de l'inverse du tenseur volumétrique
rmeso = np.zeros([3, self.AbscSize]) #Position du mesocentre
self._DeltaOrbites = np.zeros([4,3, self.AbscSize])
for i, pos in enumerate(sorted(self.poskeys)) :
for sat in self.satkeys : rmeso[i,:] += 0.25 * self.Data[sat +'_' + pos][:] #calculating Mesocentre position
for num, sat in enumerate(self.satkeys) : self._DeltaOrbites[num,i,:] = np.copy(self.Data[sat + '_' + pos][:] - rmeso[i,:]) #delta entre Sat et mesocentre
# Constructing volumetric tensor and its inverse
R = np.zeros([3,3, self.AbscSize])
self._Rinv = np.zeros([3,3,self.AbscSize])
for data in range(0, self.AbscSize) :
for i in range(0,3) :
for j in range(0,3) :
R[i,j,data] = np.mean(self._DeltaOrbites[:,i,data] * self._DeltaOrbites[:,j,data])
self._Rinv[:,:,data] = np.linalg.inv(R[:,:,data])
return self._Rinv
def getGradient(self):
Rinv = self.getMesoCentre()
gradient = {}
for comp in ['x','y','z'] : gradient[comp] = np.zeros(self.AbscSize)
for data in range(0, self.AbscSize) :
for numcomp, comp in enumerate(['x','y','z']) :
for k in range(0,3) :
for s1n, s1 in enumerate(self.satkeys) :
for s2n, s2 in enumerate(self.satkeys[s1n + 1:], start = s1n + 1) :
gradient[comp][data] += 1./16. * ((self.Data[s1 + '_' + self.datakey][data] - self.Data[s2 + '_' + self.datakey][data]) * (self._DeltaOrbites[s1n, k, data] - self._DeltaOrbites[s2n, k, data])) * self._Rinv[k,numcomp,data]
self.Data['gradient_' + self.datakey] = gradient
return gradient