text: ////////////////////////////////////////////////////////////////////////
  //
  // Implementation of oscillations of neutrinos in matter in a
  // three-neutrino framework. 
  //
  // coelho@lal.in2p3.fr
  ////////////////////////////////////////////////////////////////////////
  
  #include "PMNS_Fast.h"
  #include "MatrixDecomp/zheevh3.h"
  
  #include <iostream>
  #include <cassert>
  #include <stdlib.h>
  
  using namespace OscProb;
  
  //......................................................................
  ///
  /// Constructor. \sa PMNS_Base::PMNS_Base
  ///
  /// This class is restricted to 3 neutrino flavours.
  ///
  PMNS_Fast::PMNS_Fast() : PMNS_Base(), fHam() {}
  
  //......................................................................
  ///
  /// Nothing to clean.
  ///
  PMNS_Fast::~PMNS_Fast(){}
  
  //......................................................................
  ///
  /// Set all mixing parameters at once.
  /// @param th12    - The value of the mixing angle theta_12
  /// @param th23    - The value of the mixing angle theta_23
  /// @param th13    - The value of the mixing angle theta_13
  /// @param deltacp - The value of the CP phase delta_13
  ///
  void PMNS_Fast::SetMix(double th12, double th23, double th13, double deltacp) 
  {
  
    SetAngle(1,2, th12);
    SetAngle(1,3, th13);
    SetAngle(2,3, th23);
    SetDelta(1,3, deltacp);
  
  }
  
  //......................................................................
  ///
  /// Set both mass-splittings at once.
  ///
  /// These are Dm_21 and Dm_32 in eV^2.\n
  /// The corresponding Dm_31 is set in the class attributes.
  ///
  /// @param dm21 - The solar mass-splitting Dm_21 
  /// @param dm32 - The atmospheric mass-splitting Dm_32 
  ///
  void PMNS_Fast::SetDeltaMsqrs(double dm21, double dm32) 
  {
  
    SetDm(2, dm21);
    SetDm(3, dm32 + dm21);
  
  }
  
  //......................................................................
  ///
  /// Build the full Hamiltonian in matter.
  ///
  /// Here we divide the mass squared matrix Hms by the 2E
  /// to obtain the vacuum Hamiltonian in eV. Then, the matter
  /// potential is added to the electron component.
  ///
  void PMNS_Fast::UpdateHam()
  {
  
    double lv = 2 * kGeV2eV*fEnergy;     // 2E in eV 
  
    double kr2GNe = kK2*M_SQRT2*kGf;
    kr2GNe *= fPath.density * fPath.zoa; // Matter potential in eV
  
    // Finish building Hamiltonian in matter with dimension of eV
    for(int i=0;i<fNumNus;i++){
      fHam[i][i] = fHms[i][i]/lv;
      for(int j=i+1;j<fNumNus;j++){
        if(!fIsNuBar) fHam[i][j] = fHms[i][j]/lv;
        else          fHam[i][j] = conj(fHms[i][j])/lv;
      }
    }
    if(!fIsNuBar) fHam[0][0] += kr2GNe;
    else          fHam[0][0] -= kr2GNe;
  
  }
  
  //......................................................................
  ///
  /// Solve the full Hamiltonian for eigenvectors and eigenvalues.
  ///
  /// This is using a method from the GLoBES software available at
  /// http://www.mpi-hd.mpg.de/personalhomes/globes/3x3/ \n
  /// We should cite them accordingly
  ///
  void PMNS_Fast::SolveHam()
  {
  
    // Do vacuum oscillation in low density
    if(fPath.density < 1.0e-6){
      SetVacuumEigensystem();
      return;
    }
  
    // Build Hamiltonian
    BuildHms();
  
    // Check if anything changed  
    if(fGotES) return;
  
    // Try caching if activated
    if(TryCache()) return;
  
    UpdateHam();
  
    double fEvalGLoBES[3];
    complexD fEvecGLoBES[3][3];
  
    // Solve Hamiltonian for eigensystem using the GLoBES method
    zheevh3(fHam,fEvecGLoBES,fEvalGLoBES);
  
    // Fill fEval and fEvec vectors from GLoBES arrays  
    for(int i=0;i<fNumNus;i++){
      fEval[i] = fEvalGLoBES[i];
      for(int j=0;j<fNumNus;j++){
        fEvec[i][j] = fEvecGLoBES[i][j];
      }
    }
    
    fGotES = true;
  
    // Fill cache if activated
    FillCache();
  
  }
  
  //.....................................................................
  ///
  /// Set the eigensystem to the analytic solution in vacuum.
  ///
  /// We know the vacuum eigensystem, so just write it explicitly
  ///
  void PMNS_Fast::SetVacuumEigensystem()
  {
  
    double  s12, s23, s13, c12, c23, c13;
    complexD idelta(0.0, fDelta[0][2]);
    if(fIsNuBar) idelta = conj(idelta);
  
    s12 = sin(fTheta[0][1]);  s23 = sin(fTheta[1][2]);  s13 = sin(fTheta[0][2]);
    c12 = cos(fTheta[0][1]);  c23 = cos(fTheta[1][2]);  c13 = cos(fTheta[0][2]);
  
    fEvec[0][0] =  c12*c13;
    fEvec[0][1] =  s12*c13;
    fEvec[0][2] =  s13*exp(-idelta);
  
    fEvec[1][0] = -s12*c23-c12*s23*s13*exp(idelta);
    fEvec[1][1] =  c12*c23-s12*s23*s13*exp(idelta);
    fEvec[1][2] =  s23*c13;
  
    fEvec[2][0] =  s12*s23-c12*c23*s13*exp(idelta);
    fEvec[2][1] = -c12*s23-s12*c23*s13*exp(idelta);
    fEvec[2][2] =  c23*c13;
  
    fEval[0] = 0;
    fEval[1] = fDm[1] / (2 * kGeV2eV*fEnergy);
    fEval[2] = fDm[2] / (2 * kGeV2eV*fEnergy);
  
  }
  
  ////////////////////////////////////////////////////////////////////////
  