PMNS_Sterile.cxx
void PMNS_Sterile::SolveHam()
{

  // Build Hamiltonian
  BuildHms();

  // Check if anything changed
  if(fGotES) return;
  
  // Try caching if activated  
  if(TryCache()) return;

  double rho = fPath.density;
  double zoa = fPath.zoa;

  double lv = 2 * kGeV2eV*fEnergy;                 // 2E in eV 
  double kr2GNe = kK2*M_SQRT2*kGf * rho*zoa;       // Electron matter potential in eV
  double kr2GNn = kK2*M_SQRT2*kGf * rho*(1-zoa)/2; // Neutron matter potential in eV

  // Finish building Hamiltonian in matter with dimension of eV
  for(size_t i=0;i<size_t(fNumNus);i++){
    complexD buf = fHms[i][i]/lv;
    *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_rect(real(buf),0);
    for(size_t j=i+1;j<size_t(fNumNus);j++){
      buf = fHms[i][j]/lv;
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
    }
    if(i>2){
      // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example. 
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) ,  kr2GNn);
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) , -kr2GNn);;
    }
  }
  // Add nue CC coherent forward scattering from sterile neutrinos. 
  if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) ,  kr2GNe);
  else          *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) , -kr2GNe);

  // Solve Hamiltonian for eigensystem
  gsl_eigen_hermv(fGSL.H_GSL, fGSL.fEvalGSL, fGSL.fEvecGSL, fGSL.W_GSL);

  // Fill fEval and fEvec vectors from GSL objects
  for(int i=0;i<fNumNus;i++){
    fEval[i] = gsl_vector_get(fGSL.fEvalGSL,i);
    for(int j=0;j<fNumNus;j++){
      gsl_complex buf = gsl_matrix_complex_get(fGSL.fEvecGSL,i,j);
      fEvec[i][j] = complexD( GSL_REAL(buf), GSL_IMAG(buf) );
    }
  }
  
  fGotES = true;
  
  // Fill cache if activated  
  FillCache();
  
}
//......................................................................
///
/// Build and solve the full Hamiltonian for eigenvectors and eigenvalues.
///
/// This uses GSL to solve the eigensystem.
///
/// 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 and sterile components.
///
/// The matter potential is implemented as in arXiv:hep-ph/0606054v3 (eq. 3.15).
///
Build and solve the full Hamiltonian for eigenvectors and eigenvalues.

This uses GSL to solve the eigensystem.

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 and sterile components.

The matter potential is implemented as in arXiv:hep-ph/0606054v3 (eq.
Build and solve the full Hamiltonian for eigenvectors and eigenvalues . This uses GSL to solve the eigensystem . 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 and sterile components . The matter potential is implemented as in arXiv : hep - ph / 0606054v3 ( eq .
void PMNS_Sterile::SolveHam()
{
  BuildHms();
  if(fGotES) return;
  if(TryCache()) return;
  double rho = fPath.density;
  double zoa = fPath.zoa;
  double lv = 2 * kGeV2eV*fEnergy;                 
  double kr2GNe = kK2*M_SQRT2*kGf * rho*zoa;       
  double kr2GNn = kK2*M_SQRT2*kGf * rho*(1-zoa)/2; 
  for(size_t i=0;i<size_t(fNumNus);i++){
    complexD buf = fHms[i][i]/lv;
    *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_rect(real(buf),0);
    for(size_t j=i+1;j<size_t(fNumNus);j++){
      buf = fHms[i][j]/lv;
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
    }
    if(i>2){
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) ,  kr2GNn);
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) , -kr2GNn);;
    }
  }
  if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) ,  kr2GNe);
  else          *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) , -kr2GNe);
  gsl_eigen_hermv(fGSL.H_GSL, fGSL.fEvalGSL, fGSL.fEvecGSL, fGSL.W_GSL);
  for(int i=0;i<fNumNus;i++){
    fEval[i] = gsl_vector_get(fGSL.fEvalGSL,i);
    for(int j=0;j<fNumNus;j++){
      gsl_complex buf = gsl_matrix_complex_get(fGSL.fEvecGSL,i,j);
      fEvec[i][j] = complexD( GSL_REAL(buf), GSL_IMAG(buf) );
    }
  }
  fGotES = true;
  FillCache();
}
void PMNS_Sterile :: SolveHam ( ) { BuildHms ( ) ; if ( fGotES ) return ; if ( TryCache ( ) ) return ; double rho = fPath . density ; double zoa = fPath . zoa ; double lv = 2 * kGeV2eV * fEnergy ; double kr2GNe = kK2 * M_SQRT2 * kGf * rho * zoa ; double kr2GNn = kK2 * M_SQRT2 * kGf * rho * ( 1 - zoa ) / 2 ; for ( size_t i = 0 ; i < size_t ( fNumNus ) ; i ++ ) { complexD buf = fHms [ i ] [ i ] / lv ; * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_rect ( real ( buf ) , 0 ) ; for ( size_t j = i + 1 ; j < size_t ( fNumNus ) ; j ++ ) { buf = fHms [ i ] [ j ] / lv ; if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , j , i ) = gsl_complex_rect ( real ( buf ) , - imag ( buf ) ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , j , i ) = gsl_complex_rect ( real ( buf ) , imag ( buf ) ) ; } if ( i > 2 ) { if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , i , i ) , kr2GNn ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , i , i ) , - kr2GNn ) ; ; } } if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , 0 , 0 ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , 0 , 0 ) , kr2GNe ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , 0 , 0 ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , 0 , 0 ) , - kr2GNe ) ; gsl_eigen_hermv ( fGSL . H_GSL , fGSL . fEvalGSL , fGSL . fEvecGSL , fGSL . W_GSL ) ; for ( int i = 0 ; i < fNumNus ; i ++ ) { fEval [ i ] = gsl_vector_get ( fGSL . fEvalGSL , i ) ; for ( int j = 0 ; j < fNumNus ; j ++ ) { gsl_complex buf = gsl_matrix_complex_get ( fGSL . fEvecGSL , i , j ) ; fEvec [ i ] [ j ] = complexD ( GSL_REAL ( buf ) , GSL_IMAG ( buf ) ) ; } } fGotES = true ; FillCache ( ) ; }
......................................................................
......................................................................
// Build Hamiltonian // Check if anything changed // Try caching if activated   // 2E in eV  // Electron matter potential in eV // Neutron matter potential in eV // Finish building Hamiltonian in matter with dimension of eV // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example.  // Add nue CC coherent forward scattering from sterile neutrinos.  // Solve Hamiltonian for eigensystem // Fill fEval and fEvec vectors from GSL objects // Fill cache if activated  