#ifndef PRM_CALCULATION_HH
#define PRM_CALCULATION_HH

#include <vector>
#include <iostream>

#include <armadillo>

#include "BasisState.hh"

#define PI 3.14159265358979323846

namespace PRM {
  class TensorOperator {
    //this is a holder class for the evaluated decoupling factors
  public:
    TensorOperator() {};    
    TensorOperator(int dim, int rank);
    void Set(int dim, int rank);
    int fRank;
    std::vector<arma::mat> values;

    arma::mat &GetValue(int rank);
    arma::mat &operator()(int rank);
  };
  
  class Nilsson {
  public:
    std::vector<BasisState> fBasisStates;
    int fZ;
    int fA;
    int fN;
    int fNNeut;
    int fNProt;
    int fOdd;
    int fNOdd;
    int fMaxOmg;

    std::vector<double> fKappa;
    std::vector<double> fMu;

    double fEpsP;
    double fGamma;
    double fEps2;
    double fEps22;
    double fEps4;
    double fEps6;

    double fOmRot;
    double fOmOm;
    double fFac;
    double fHac;

    int fCoupledShells;

    std::vector<int> fAvL2;
    std::vector<double> fGammaFunc;

    arma::mat fHamiltonian;
    arma::vec fEnergies;
    std::vector<int> fOmegas;
    arma::mat fCompositions;

    //decopuling factors
    std::vector<int> fDecoupleLevels;
    TensorOperator fABCDE; //this houses what are labeled a+b,c,d,e in gampn
    //a+b=1,c=-1,d=0,e=2
    TensorOperator fJz;
    TensorOperator fJz2;
    TensorOperator fS;

    TensorOperator fJr;
    TensorOperator fSr;
    TensorOperator fFr;
    
    TensorOperator fR2J;
    TensorOperator fR4J;
    TensorOperator fR2S;
    TensorOperator fR4S;
    TensorOperator fR2D;
    TensorOperator fR4D;
    
    TensorOperator fNorm;

    arma::vec fQ20s;
    arma::vec fQ22s;
    arma::vec fQ20;
    arma::vec fQ22;
    arma::vec fRho2;
    arma::vec fR2;
    arma::vec fRJz;    

    Nilsson(int Z, int A);
    double Gamma(int i);
    void SetNMax(int nn, int np);
    void MakeBasisStates(int parity);
    void PrintBasisStates(std::ostream &out, int ncols=5);

    void MakeHamiltonian();
    void Diagonalize();
    void Solve();

    void SetDecoupleLevels(std::vector<int> levels);
    void Decouple();
    void Decouple(std::vector<int> levels);

    void DiagQMom();      
    
    void PrintWavefunctions(std::ostream &out, int maxrow=5);
    void PrintDecouple(std::ostream &out);
    void PrintDiagQMom(std::ostream &out);
    void WriteEnergies(std::ostream &out);
    int GetMaxOmega();

    void SetKappas(std::vector<double> kappas);
    void SetMus(std::vector<double> mus);
    void SetDeformation(double epsp, double gamma, double eps4, double eps6);
    void VolumeConservation();

    void SetCranking(double omrot);

    double ClebschGordan(int j1, int j2, int J, int m1, int m2, int M);
    double GetRadialME(int Ni, int li, int Nj, int lj, int T);
    double RedMED1(int l2, int s2, int j2, int l2p, int s2p, int j2p);
    

  };
}   

#endif
