#ifndef PRM_PARTICLE_ROTOR_HH
#define PRM_PARTICLE_ROTOR_HH

#include <iostream>
#include <map>

#include <armadillo>

#include "Nilsson.hh"
#include "Pairing.hh"

namespace PRM {
  class ParticleRotor {
  public:
    Nilsson &fNilsson;
    Pairing &fPair;
    int fA,fZ,fN;

    double fGamma;    

    int fMinSpin;
    int fMaxSpin;

    std::vector<int> fKs;
    std::map<int,arma::mat> fHamiltonians;
    std::map<int,arma::vec> fEnergies;
    std::map<int,arma::mat> fCompositions;
    std::map<int,std::vector<int> > fNStates;

    double fHb2B;
    double fA1,fA2,fA3;
    double fCoriolis;
    
    double fEfac;
    bool fGam;
    bool fVMI;
    double fEta;
    std::vector<double> fU;
    std::vector<double> fV;
    std::vector<double> fSQP;
    arma::mat fPUV;

    ParticleRotor(Nilsson &nil, Pairing &pair);

    void SetB(double E2plus);
    void SetBGrodzin();
    void SetBSimple(double E2plus);
    void SetMomentsOfInertia(double A1, double A2, double A3);
    void SetMomentsOfInertia();
    
    double Inertia(int Ir2, double A00, double stiff);
    void MakeHamiltonians(std::ostream &out);
    void Diagonalize();

    double GetKSum(int spin,int indx, int k);
    double GetSPSum(int spin, int indx, int ni);
    int GetKDom(int spin, int indx);
    int GetSPDom(int spin, int indx);
    std::pair<int,int> GetMaxKSP(int spin, int indx);
   
    void PrintMomentsOfInertia(std::ostream &out);
    void PrintSPLevels(std::ostream &out);
    void PrintResults(std::ostream &out, double maxEn = -1);
    void PrintAmplitudes(std::ostream &out , std::vector<std::pair<int, int> > states);
    
  };
}

#endif
