///////////////////////////////////////////////////////////////////////////////
//
// Implementation of oscillations of neutrinos with scalar NSI
// three-neutrino framework.
//
// jcoelho@apc.in2p3.fr and urahaman@km3net.de
///////////////////////////////////////////////////////////////////////////////

#include "PMNS_SNSI.h"

using namespace OscProb;

//.............................................................................
///
/// Constructor. \sa PMNS_NSI::PMNS_NSI
///
/// This class is restricted to 3 neutrino flavours.
///
PMNS_SNSI::PMNS_SNSI() : PMNS_NSI() { SetLowestMass(0); }

//.............................................................................
///
/// Nothing to clean.
///
PMNS_SNSI::~PMNS_SNSI() {}

//.............................................................................
///
/// Set lightest neutrino mass
///
/// @param m - The lightest mass in eV
///
void PMNS_SNSI::SetLowestMass(double m)
{
  // Check if value is actually changing
  fBuiltHms *= (fM == m);

  fM = m;
}

//.............................................................................
///
/// Get lightest neutrino mass
///
/// @return The lightest mass in eV
///
double PMNS_SNSI::GetLowestMass() { return fM; }

//.............................................................................
///
/// Reimplement vacuum Hms update
///
/// \sa PMNS_Base::BuildHms
///
void PMNS_SNSI::BuildHms()
{
  // Check if anything changed
  if (fBuiltHms) return;

  // Tag to recompute eigensystem
  fGotES = false;

  // Start with m1 = 0
  double m1 = 0;
  // Make sure all masses are positive
  for (int i = 1; i < fNumNus; i++) {
    if (fDm[i] + m1 * m1 < 0) m1 = sqrt(fabs(fDm[i]));
  }
  // Add minimum mass
  m1 += fM;

  // Build M - m1
  fHms[0][0] = 0;
  for (int j = 1; j < fNumNus; j++) {
    // Set mass difference m_j - m_1
    fHms[j][j] = sqrt(fDm[j] + m1 * m1) - m1;
    // Reset off-diagonal elements
    for (int i = 0; i < j; i++) { fHms[i][j] = 0; }
    // Rotate j neutrinos
    for (int i = 0; i < j; i++) { RotateH(i, j, fHms); }
  }

  // Add back m1
  for (int i = 0; i < fNumNus; i++) { fHms[i][i] += m1; }

  // Tag as built
  fBuiltHms = true;
}

//.............................................................................
///
/// Square a Hermitian matrix
///
/// @param A - input matrix A
///
void HermitianSquare(complexD (&A)[3][3])
{
  complexD tmp[3][3];
  for (int i = 0; i < 3; i++) {
    for (int j = i; j < 3; j++) {
      tmp[i][j] = A[i][j];
      if (i < j) tmp[j][i] = conj(A[i][j]);
      A[i][j] = 0;
    }
  }

  for (int i = 0; i < 3; i++) {
    for (int j = i; j < 3; j++) {
      for (int k = 0; k < 3; k++) { A[i][j] += tmp[i][k] * tmp[k][j]; }
    }
  }
}

//.............................................................................
///
/// Reimplement Hamiltonian update for scalar NSI
///
/// 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_SNSI::UpdateHam()
{
  double sqrtlv = sqrt(2 * kGeV2eV * fEnergy); // sqrt(2*E) in eV^1/2

  // Effective density of fermions in eV * MeV^2
  double nsiCoup = 1e6 * kK2 * fPath.density * GetZoACoup();

  // Add the dM matrix
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) {
      fHam[i][j] = (fHms[i][j] + nsiCoup * fEps[i][j]) / sqrtlv;
      if (fIsNuBar) fHam[i][j] = conj(fHam[i][j]);
    }
  }

  // Square fHam
  HermitianSquare(fHam);

  // Add matter potential in eV
  double kr2GNe = kK2 * M_SQRT2 * kGf * fPath.density * fPath.zoa;

  if (!fIsNuBar)
    fHam[0][0] += kr2GNe;
  else
    fHam[0][0] -= kr2GNe;
}

///////////////////////////////////////////////////////////////////////////////