// ------------------------------------------------------------------
//
//  GenericMinVar.h
//  Qdos2 - Generic Constrained Minimum Variance Algorithm for. Based on 
//  the original MinVariance.h.
//
//  Janet Barnes
//  May 2000
//
//  This version of Minimum Variance contains an additional projection
//  matrix P, which, when constrain is enabled, is projected from both 
//  sides on the covariance  matrix, thus creating a conversion matrix 
//  (eigenvecs) which  ensures that the time series used to constrain has 
//  a zero average along the direction corresponding to the normal 
//  direction. See blue ISSI book on Multispacecraft Analysis, p193. 
//
//  Cloned & modified by 
//  SEH,    May 2002 - cloned from original MinVariance.h
//  SEH,  6 Sep 2003 - Generic version version, 
//  SEH, 11 May 2005 - corrected RH correction; 
//                     
// ------------------------------------------------------------------

#ifndef QDOS_GENERICMINVAR_H
#define QDOS__GENERICMINVAR_H

#define VERBOSE              1

#include "Algorithm.h"
#include "Matrix.cc"

namespace QDOS{

class GenericMinVar
{
    Matrix<double> _eigvals;       // Matrix Dim(3) of eigenvalues
    Matrix<double> _eigvecs;       // Matrix Dim(3,3) of eigenvectors (in rows)
    Matrix<double> _eigqual;       // Matrix Dim(3) eigenvalue quality
    double _orgqual;               // orthogonality quality
    Matrix<double> _covariance;    // covariance matrix
    Matrix<double> Q;              // MFR Q matrix for which eigenvals _ eigenvecs are found
    Matrix<double> temp;           // temporary storage matrix
    int i,j;
    
public:
    GenericMinVar(int constrained, Matrix<double> Q, Matrix<double> P):_eigvals(Matrix<double>(Dim(3))),\
        _eigvecs(Matrix<double>(Dim(3,3))),_covariance(Matrix<double>(Dim(3,3))),\
      _eigqual(Matrix<double>(Dim(3)))
    {
 
        _covariance.resize(Dim(3,3));

        _covariance = Q;        
   

        _covariance[1][0] = _covariance[0][1];
        _covariance[2][0] = _covariance[0][2];
        _covariance[2][1] = _covariance[1][2];
        
	// 
	// Constrained MVA : find eigenvec of [P # C # P] * n = lamda * n
	// instead of [C] * n = lambda * n. 
	//
	//   0 - standard MVA
	//   1 - constrain so that Bn == 0 (see p 209 in blue ISSI book)
        //       P matrix calculated here
	//   2 - constrain using given P matrix
	//
	// SEH, 14 May 2002
	//      28 May 2002, MFR version
	//      26 Mar 2003, MFL version
	//

	if (constrained >= 1) {

          temp.resize(Dim(3,3));

	  temp = _covariance*P;
	  _covariance = P*temp;

  	// force symmetric covariance matrix again
	  _covariance[1][0] = _covariance[0][1];
	  _covariance[2][0] = _covariance[0][2];
	  _covariance[2][1] = _covariance[1][2];
	}

        // get eigenvectors and values  
        MatrixEig eig(_covariance);
        _eigvecs=eig.eigenvectors().seq_to_matrix();
        _eigvals=eig.eigenvalues().seq_to_matrix();


	// If constrained : eigval 0 == 0. Swap 1st and 2nd eigenvec

        if (constrained >= 1) {
          _eigvals[0] = _eigvals[1];
          _eigvals[1] = _eigvals[2];
          _eigvals[2] = 0;                                // by definition
          temp = _eigvecs;
          for (i=0; i<3; i++) _eigvecs[0][i] = temp[1][i];
          for (i=0; i<3; i++) _eigvecs[1][i] = temp[2][i];
          for (i=0; i<3; i++) _eigvecs[2][i] = temp[0][i];
        }

        
        // check the positive sense of the eigenvector associated with
        // the minimum eigenvalue and if necessary reverse the signs
        // of all three components i.e. vint x vmax = vector perpendicular
        // to both vint and vmax. If vmin is perpendicular to vint 
        // and vmax (i.e. system is ortogonal), then (vinx x vmax).vmin=1
        // If < 0 then not a right-handed system, so reverse signs of vmin
        double rh_factor=(_eigvecs[1].cross(_eigvecs[2])).dot(_eigvecs[0]);
        if(rh_factor<0.)
            _eigvecs[0]*=-1;
        _orgqual=std::abs(std::abs(rh_factor)-1.0);
                
        // ensure that the sign of the average minimum variance component
        // is positive by computing average of sum of B-field data.min
        // eigenvector. If less than zero negate both vmin and vint (so
        // preserving right-handed system)


        // SEH 11 May 2005 : this was already commented out; a positive normal comp
	// does not make sense for us.

        _eigqual=_eigvals/_eigvals[2];
     
        // SEH 2003_03_26: for MFL, MFR, pr definition, the eigenvector describing
	// the normal component should have a positive X-value; If not,
	// flip the eigenvector matrix, conserve RH system

	// SEH 11 May 2005 : Correction; to preserve RH sys, only two eigvals should be flipped.
	if (_eigvecs[0][0] < 0.0){
            _eigvecs[0]*=-1.0;
            _eigvecs[1]*=-1.0;
	    //            _eigvecs[2]*=-1.0;
	}
       
    }



    // methods for returning members and calculating variance
    Matrix<double>& covariance() 
    { return _covariance; }
    
    Matrix<double>& eigenvalues() 
    { return _eigvals;}
     
    Matrix<double>& eigenvectors()
    { return _eigvecs;}   
    
    Matrix<double>& eigenquality() 
    { return _eigqual;}
    
    double& orthogonality_quality()
    { return _orgqual;}    // orthogonality quality     

        
};
    
}  // end namespace QDOS
#endif
