#! /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()