mva.py
3.13 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
#! /bin/env python
# coding: utf8
import numpy as np
import matplotlib.pyplot as plt
class mva(object):
"""
Methode de calcul du minimum variance analysis
Input field : magnetic field or any other 3D field to apply the MVA method in NUMPY ARRAY
X,Y,Z in rows and data in columns such as B[:,0] = all data, x component
input field2 : table 3x2 first row is magnetosphere vector, second is magnetosheath then MVA is determined using cross product between magnetosheath and magnetosphere, such as B[0,:] = MSP vector en NUMPY ARRAY
TO DO : option to input only 1 timeserie and MSP, MSH, current sheet time interval using decorator?)
Results are provided in self.results
MVA provide a transport matrix, with L,M,N in each colums of the matrix, such as :
new_vect = old_vect([x,y,z]) * matrix, or new_vect = inv(matrix) * old_vect[[x], [y], [z]]
"""
def __init__(self, field, field2 = None, noprint = False):
self.field = field #np.array([[x,]]) #en entrée : champ avec 3 composantes sur n mesures
self.results = {}
self.calc(noprint = noprint)
#self.noprint = noprint
if field2 != None :
self.field_cross = field2 # en entré : deux vecteurs de 3 composantes déjà moyenné
self.calc_cross(noprint = noprint)
def calc(self, noprint = False) :
self.matrice = np.zeros([3,3])
for i in range(3) :
for j in range(3) :
self.matrice[i,j] = np.mean(self.field[:,i] * self.field[:,j]) - self.field[:,i].mean() * self.field[:,j].mean()
eigVal, eigVect = np.linalg.eig(self.matrice)
Vec_N, Vec_M, Vec_L = [eigVect[:,i] for i in np.argsort(abs(eigVal))]
self.results['eigen_Val'] = eigVal[np.argsort(abs(eigVal))]
self.results['eigen_Vect'] = {'l' : Vec_L, 'm' : Vec_M, 'n' : Vec_N}
self.results['Matrice_Passage'] = np.mat([Vec_L, Vec_M, Vec_N]).T
if not noprint :
print 'Resultat par diagonalisation de Matrice : ', '\n', 'L : ', Vec_L, 'Eigen val :', self.results['eigen_Val'][2] , '\n', 'M : ', Vec_M, 'Eigen val :', self.results['eigen_Val'][1], '\n', 'N : ', Vec_N, 'Eigen val :', self.results['eigen_Val'][0]
def calc_cross(self, noprint = False) :
Vec_MSP, Vec_MSH = [self.field_cross[0,:], self.field_cross[1,:]]
Vec_N2 = np.cross(Vec_MSP, Vec_MSH) / np.linalg.norm(np.cross(Vec_MSP,Vec_MSH))
tmp = np.cross(Vec_N2, self.results['eigen_Vect']['l'])
Vec_M2 = (tmp / np.linalg.norm(tmp)) * np.vdot(tmp, Vec_MSP) / abs(np.vdot(tmp, Vec_MSP))
Vec_L2 = np.cross(Vec_M2, Vec_N2)
self.results['eigen_Vect2'] = {'l' : Vec_L2, 'm' : Vec_M2, 'n' : Vec_N2}
self.results['Matrice_Passage2'] = np.mat([Vec_L2, Vec_M2, Vec_N2]).T
if not noprint :
print 'Resultat par produit croisé SH / SP : \n', 'L : ', Vec_L2, '\n', 'M : ', Vec_M2, '\n', 'N : ', Vec_N2
def hodogram(self, system = 1) :
"""
project field in MVA new coordinate system and plot the L vs N and M vs N graphs
"""
if system == 1 : matrice = self.results['Matrice_Passage']
elif system == 2 : matrice = self.results['Matrice_Passage2']
else : print 'Choose MVA 1 or 2'
proj = np.array([vect *matrice for vect in self.field])[:,0]
f = plt.figure()
plt.plot(proj[:,0], proj[:,2])
plt.plot(proj[:,1], proj[:,2])
plt.show()