#include <iomanip>

#include "ParticleRotor.hh"

namespace PRM {
  ParticleRotor::ParticleRotor(Nilsson &nil, Pairing &pair) :
    fNilsson(nil),
    fPair(pair),
    fEfac(pair.fEfac),
    fGamma(nil.fGamma),
    fA(nil.fA),
    fZ(nil.fZ),
    fN(nil.fN),
    fGam(false),
    fVMI(false),
    fEta(1.0),
    fMinSpin(1),
    fMaxSpin(17),
    fCoriolis(1),
    fHb2B(0)
  {
    if ( fGamma != 0) { fGam = true; }
    else { fGam = false; }
    
    int dim = fNilsson.fEnergies.size();
    for (int i=0; i<dim; ++i) {
      double nilen = fNilsson.fEnergies[i];
      if (nilen < fPair.fEnergies[fPair.fImin]/fEfac) {
        fU.push_back(0.0);
        fV.push_back(1.0);
        fSQP.push_back(nilen*fEfac);
      }
      else if (nilen > fPair.fEnergies[fPair.fImax]/fEfac) {
        fU.push_back(1.0);
        fV.push_back(0.0);
        fSQP.push_back(nilen*fEfac);
      }
      else {
        fSQP.push_back(std::sqrt(std::pow(nilen*fEfac - fPair.fEferm, 2) +
                                 std::pow(fPair.fDelta, 2)));
        fU.push_back(std::sqrt(0.5*(1.0 + (nilen*fEfac - fPair.fEferm)/fSQP[i])));
        fV.push_back(std::sqrt(0.5*(1.0 - (nilen*fEfac - fPair.fEferm)/fSQP[i])));
      }
    }

    if (fNilsson.fDecoupleLevels.size() == 0) { return; }
    
    fPUV = arma::mat(dim,dim,arma::fill::zeros);
    for (int ini=0; ini<fNilsson.fDecoupleLevels.size(); ++ini) {
      int ni = fNilsson.fDecoupleLevels[ini];
      for (int inj=0; inj<fNilsson.fDecoupleLevels.size(); ++inj) {
        int nj = fNilsson.fDecoupleLevels[inj];
        fPUV(ni,nj) = std::pow(fU[ni]*fU[nj] + fV[ni]*fV[nj], fEta);

        fNilsson.fABCDE(0)(ni,nj) *= fPUV(ni,nj);
        fNilsson.fABCDE(-1)(ni,nj) *= fPUV(ni,nj);
        fNilsson.fABCDE(1)(ni,nj) *= fPUV(ni,nj);
        fNilsson.fABCDE(2)(ni,nj) *= fPUV(ni,nj);

        fNilsson.fJz(0)(ni,nj) *= fPUV(ni,nj);

        fNilsson.fJz2(0)(ni,nj) *= fPUV(ni,nj);
      }
    }
    
    if (!fVMI) {
      fNilsson.fJz2(0).zeros();
      fNilsson.fABCDE(0).zeros();
      fNilsson.fABCDE(2).zeros();

      for (int i=0; i<dim; ++i) {
        for (int ini=0; ini<fNilsson.fDecoupleLevels.size(); ++ini) {
          int ni = fNilsson.fDecoupleLevels[ini];
          for (int inj=ini; inj<fNilsson.fDecoupleLevels.size(); ++inj) {
            int nj = fNilsson.fDecoupleLevels[inj];
            double pfac = ((fU[i]*fU[ni] + fV[i]*fV[ni])*
                           (fU[i]*fU[nj] + fV[i]*fV[nj]) -
                           (fU[i]*fV[ni] - fV[i]*fU[ni])*
                           (fU[i]*fV[nj] - fV[i]*fU[nj]));

            fNilsson.fABCDE(0)(ni,nj) += (fNilsson.fJr(1)(ni,i)*fNilsson.fJr(1)(nj,i) +
                                         fNilsson.fJr(-1)(ni,i)*fNilsson.fJr(-1)(nj,i))*pfac/2.;
            fNilsson.fABCDE(2)(ni,nj) += (fNilsson.fJr(1)(ni,i)*fNilsson.fJr(-1)(nj,i) +
                                         fNilsson.fJr(-1)(ni,i)*fNilsson.fJr(1)(nj,i))*pfac;
            fNilsson.fJz2(0)(ni,nj) += fNilsson.fJr(0)(ni,i)*fNilsson.fJr(0)(nj,i)*pfac;
          }
        }
      }
      
      for (int ini=0; ini<fNilsson.fDecoupleLevels.size(); ++ini) {
        int ni = fNilsson.fDecoupleLevels[ini];
        for (int inj=ini; inj<fNilsson.fDecoupleLevels.size(); ++inj) {
          int nj = fNilsson.fDecoupleLevels[inj];
          fNilsson.fABCDE(0)(nj,ni) = fNilsson.fABCDE(0)(ni,nj);
          fNilsson.fABCDE(2)(nj,ni) = fNilsson.fABCDE(2)(ni,nj);
          fNilsson.fJz2(0)(nj,ni) = fNilsson.fJz2(0)(ni,nj);
        }
      }
    }

  };

  void ParticleRotor::SetBGrodzin() {
    fHb2B = 1225./std::pow(fA,7./3.);
  }
  
  void ParticleRotor::SetB(double E2plus) {
    /*
    double gamma = fGamma;
    //this is what is in asyrmo.f
    if (gamma == 0.0) { gamma = 0.1*PI/180.0; }
    if (gamma == PI/3.0) { gamma = 59.9*PI/180.0; }
    fHb2B = 6. * std::pow(fNilsson.fEpsP,2) * E2plus * std::pow(std::sin(3.*gamma),2) /
      (4.5*(3.-std::sqrt(9.-8.*std::pow(std::sin(3.*gamma),2))));
    */

    //analytically equivalent to
    fHb2B = 4./3. * std::pow(fNilsson.fEpsP,2) * E2plus * ((3. + std::sqrt(9. - 8.*std::pow(std::sin(3.*fGamma),2)))/8.);
    //which has no instability at gamma = 0 or gamma = 60
    

  }

  void ParticleRotor::SetBSimple(double E2plus) {
    fHb2B = E2plus/std::pow(fA,2);
  }
  
  void ParticleRotor::SetMomentsOfInertia(double A1, double A2, double A3) {
    fA1 = A1;
    fA2 = A2;
    fA3 = A3;
  }

  void ParticleRotor::SetMomentsOfInertia() {
    if (std::abs(fGamma*180.0/PI) < 1.0) {
      fA1 = 1./3. * fHb2B/2. / std::pow(fNilsson.fEpsP,2);
      fA2 = fA1;
      fA3 = 1000.*fA2; //really large I guess?
    }
    else if (std::abs(fGamma*180.0/PI-60.0) < 1.0) {
      fA2 = 1./3. * fHb2B/2. / std::pow(fNilsson.fEpsP,2);
      fA3 = fA2;
      fA1 = 1000.*fA2; //really large I guess?
    }
    else {
      fA1 = 1./4./std::pow(std::sin(fGamma+2.*PI/3.),2) * fHb2B/2. / std::pow(fNilsson.fEpsP,2);
      fA2 = 1./4./std::pow(std::sin(fGamma-2.*PI/3.),2) * fHb2B/2. / std::pow(fNilsson.fEpsP,2);
      fA3 = 1./4./std::pow(std::sin(fGamma),2) * fHb2B/2. / std::pow(fNilsson.fEpsP,2);
    }
  }

  void ParticleRotor::PrintMomentsOfInertia(std::ostream &out) {
    out << "hbar^2/B = " << fHb2B*fA*fA << " /A^2 MeV" << std::endl;
    out << std::setprecision(5) << "   A1 = " << fA1 << ", A2 = " << fA2 << ", A3 = " << fA3 << std::endl;    
  }
  
  double ParticleRotor::Inertia(int Ir2, double A00, double stiff) {
    double r = Ir2*(Ir2+1)/4./stiff + 1./216./std::pow(A00, 3);
    double t = r + 1./216./std::pow(A00,3);
    t *= Ir2*(Ir2+1)/4./stiff;
    t = std::sqrt(t);
    double retval = std::pow(r+t, 1./3.) + std::pow(r-t, 1./3.) + 1./6./A00;
    retval = 1./(2*retval);
    return retval;
  }

  void ParticleRotor::PrintSPLevels(std::ostream &out) {
    if (fNilsson.fDecoupleLevels.size() == 0) { return; }
    out << " S.P. Level     EQP       U       V     (UU + VV)^" << std::setprecision(2) << fEta << std::endl;
    for (int ini=0; ini<fNilsson.fDecoupleLevels.size(); ++ini) {
      int ni = fNilsson.fDecoupleLevels[ini];
      out << std::fixed << std::setw(11) << ni+1 << " ";
      out << std::fixed << std::setprecision(3) << std::setw(7) << fSQP[ni] << " ";
      out << std::fixed << std::setprecision(3) << std::setw(7) << fU[ni] << " ";
      out << std::fixed << std::setprecision(3) << std::setw(7) << fV[ni] << "   ";
      for (int inj=0; inj<fNilsson.fDecoupleLevels.size(); ++inj) {
        int nj = fNilsson.fDecoupleLevels[inj];
        out << std::fixed << std::setprecision(3) << std::setw(7) << fPUV(ni,nj) << " ";
      }
      out << std::endl;
    }
  }

  /*
  void ParticleRotor::CalcCoreStates(std::ostream &out) {
    if (!fVMI) { fA0s[0] = 1.0; }
    
    for (int ir = 0; ir < 18+fMaxSpin; ir+=2) {

      if (fVMI) {fA0s[ir+1] = Inertia(ir,fA0s[0],fStiff);}
      else { fA0s[ir+1] = 1.0; }

      
      
    }
    }*/
  
  void ParticleRotor::MakeHamiltonians(std::ostream &out) {    
    int nlevels = fNilsson.fDecoupleLevels.size();
    if (nlevels == 0) { return; }

    if (!fGam) {
      out << "Axial symmetry => K restricted to Omega" << std::endl;
      fKs.clear();
      for (int i=0; i<nlevels; ++i) {
        int ni = fNilsson.fDecoupleLevels[i];
        if (fNilsson.fJz(0)(ni,ni) > 0) { fKs.push_back((int)(2*fNilsson.fJz(0)(ni,ni)+0.1)); }
        else { fKs.push_back((int)(2*fNilsson.fJz(0)(ni,ni)-0.1)); } //get K of the nilsson state
        out << std::fixed << std::setw(25) << "SP orbital : " << ni+1 << std::setw(6) << " " << fKs[i] << "/2" << std::endl;
      }
    }
    
    for (int spin = fMinSpin; spin <= fMaxSpin; spin+=2) {
      int spindim = (spin+1)/2;
      if (!fGam) { spindim = 1; }
      int matdim = spindim*nlevels;
      fHamiltonians[spin] = arma::mat(matdim, matdim, arma::fill::zeros);
      int kkmax = spin;
      int n = (spin+1)/2;
      if ((spin-1)%4 != 0) { kkmax = spin - 2; }
      fNStates[spin] = std::vector<int>((fMaxSpin-fMinSpin)/2+1,matdim);

      //kkmax also defined here, not sure what it is or if I need it.

      if (fVMI) { //variable moment of inertia

      }
      else { //not variable moment of inertia
        for (int ini = 0; ini < nlevels; ++ini) {
          int ni = fNilsson.fDecoupleLevels[ini];
          for (int inj = 0; inj <= ini; ++inj) {
            int nj = fNilsson.fDecoupleLevels[inj];
            double delta = 0.0;
            if (ni == nj) { delta = 1.0; }

            for (int i = 0; i<spindim; ++i) {
              int Kp = 0;
              if (!fGam) { Kp = fKs[ini]; }
              else {//something here
                Kp = kkmax - 4*i;
              }

              for (int j = 0; j<spindim; ++j) {
                int K = 0;
                if (!fGam) { K = fKs[inj]; }
                else {//something here
                  K = kkmax - 4*j;
                }

                int mati = i+spindim*ini;
                int matj = j+spindim*inj;

                if (matj > mati) { continue; }

                int ktest = std::max(std::abs(K), std::abs(Kp));
                if (ktest > spin && matj == mati && ni == nj) {
                  fHamiltonians[spin](mati,matj) = 50.*(fA1+fA2)*(spin*(spin+2.)+K*K)+fSQP[ni];
                  fHamiltonians[spin](matj,mati) = 50.*(fA1+fA2)*(spin*(spin+2.)+K*K)+fSQP[ni];
                  fNStates[spin][(spin-1)/2] -= 1;
                }
                if (ktest > spin) { continue; }

                //main matrix elements
                double me = 0;
                if (K == 1) {
                  if (Kp == 1) {
                    me = (fA1 + fA2)/2. * ((spin*(spin+2.)/4. - 1./4.) * delta +
                                           fNilsson.fABCDE(1)(ni,nj) * n * std::pow(-1, n)) +
                      (fA1 - fA2)/2. * fNilsson.fABCDE(-1)(ni,nj)*n*std::pow(-1, n) +
                      fA3*(std::pow(K/2.0, 2) * fNilsson.fNorm(0)(ni,nj) - K * fNilsson.fJz(0)(ni,nj) + fNilsson.fJz2(0)(ni,nj)) + 
                      (fA1 + fA2)/2. * fNilsson.fABCDE(0)(ni,nj) +
                      (fA1 - fA2)/4. * fNilsson.fABCDE(2)(ni,nj);
                  }
                  if (Kp == -3) {
                    me = ((fA1 - fA2)*((spin + 1.)/4.*delta +
                                              fNilsson.fABCDE(1)(ni,nj) * std::pow(-1, n)) +
                                 (fA1 + fA2)*fNilsson.fABCDE(-1)(ni,nj) * std::pow(-1, n))*
                      std::sqrt((spin-1.)*(spin+3.)/4.)/2.;
                  }                    
                }
                else {
                  if (Kp == K-4) {
                    me = (fA1 - fA2)/4.*delta*std::sqrt((spin+K)*(spin+K-2.)*(spin-K+2.)*(spin-K+4.)/16.);
                  }
                  if (Kp == -K + 2) {
                    me = ((fA1+fA2)*fNilsson.fABCDE(1)(ni,nj) + (fA1-fA2)*fNilsson.fABCDE(-1)(ni,nj))/2. *
                      std::pow(-1, n) * std::sqrt(spin * (spin+2.)/4. - K*(K-2.)/4.);
                  }
                  if (Kp == - K - 2) {
                    me = ((fA1 - fA2) * fNilsson.fABCDE(1)(ni,nj) + (fA1 + fA2) * fNilsson.fABCDE(-1)(ni,nj))/2. *
                      std::pow(-1, n) * std::sqrt(spin * (spin+2.)/4. - K*(K+2.)/4.);
                  }
                  if (Kp == K) {
                    me = (fA1 + fA2)/2. * (spin *(spin+2.)/4. - K*K/4.) * delta +
                      fA3 * (std::pow((double)K/2., 2) * fNilsson.fNorm(0)(ni,nj) - (double)K*fNilsson.fJz(0)(ni,nj) + fNilsson.fJz2(0)(ni,nj)) +
                      (fA1 + fA2)/2. * fNilsson.fABCDE(0)(ni,nj) + (fA1 - fA2)/4. * fNilsson.fABCDE(2)(ni,nj);

                  }
                }
                if (ni != nj) { me *= fCoriolis; }
                if (ni == nj && mati == matj) {
                  me += fSQP[ni];
                }
                fHamiltonians[spin](mati,matj) = me;
                fHamiltonians[spin](matj,mati) = me;
              } //end j loop
            } //end i loop
          } //end nj loop
        } //end ni loop        
      } //end VMI if
    } // end spin loop
    return;
  }

  void ParticleRotor::Diagonalize() {
    if (fNilsson.fDecoupleLevels.size() == 0) { return; }
    for (int spin = fMinSpin; spin <= fMaxSpin; spin+=2) {
      int dim = fHamiltonians[spin].n_rows;
      fEnergies[spin] = arma::vec(dim,arma::fill::zeros);
      fCompositions[spin] = arma::mat(dim,dim,arma::fill::zeros);
      
      arma::eig_sym(fEnergies[spin], fCompositions[spin], fHamiltonians[spin]);
    }
  }

  void ParticleRotor::PrintResults(std::ostream &out, double maxEn /*= -1*/) {
    if (fNilsson.fDecoupleLevels.size() == 0) { return; }
    for (int spin = fMinSpin; spin <= fMaxSpin; spin +=2) {
      //for (int ind = 1; ind <= fEnergies[spin].size(); ++ind) {
      int spindim = (spin+1)/2;
      int nlevels = fNilsson.fDecoupleLevels.size();
      int matdim = spindim*nlevels;
      int kkmax = spin;
      int n = (spin+1)/2;
      if ((spin-1)%4 != 0) { kkmax = spin - 2; }

      for (int ind = 1; ind <= fNStates[spin][(spin-1)/2]; ++ind) {
        double sum = 0;
        //if (fEnergies[spin][ind-1] > 10) { continue; }
        if (maxEn > 0 && fEnergies[spin][ind-1] > maxEn) { continue; }
        if (fGam) {
        out << "            Energy    K=   ";
          for (int i=0; i<spindim; ++i) {
            out << std::fixed << std::setw(3) << kkmax-4*i <<  "/2   ";
          }
          out << std::endl;
        }
        out << std::fixed << std::setw(2) << spin << "/2 (" << ind << ")    ";
        out << std::fixed << std::setw(6) << std::setprecision(3) << fEnergies[spin][ind-1] << "    ";
        arma::vec wf = fCompositions[spin].col(ind-1);
        int prct = 0;
        if (!fGam) {
          for (int iwf = 0; iwf < wf.size(); ++iwf) {
            sum += std::pow(wf[iwf],2);
            if (std::abs(wf[iwf])>1e-4) {
              if (prct > 0) { out << std::setw(22) << " "; }
              out << fNilsson.fDecoupleLevels[iwf]+1 << " : " << wf[iwf] << std::endl;
              ++prct;
            }
          }
          out << std::endl;
        }
        else {
          for (int isp = 0; isp < fNilsson.fDecoupleLevels.size(); ++isp) {
            double maxwf = 0.0;            
            for (int i=0; i<spindim; ++i) {
              sum += std::pow(wf[isp*spindim+i],2);
              if (std::abs(wf[isp*spindim+i]) > maxwf) {
                maxwf = std::abs(wf[isp*spindim+i]);
              }              
            }
            if (maxwf > 0.05) {
              if (prct > 0) { out << std::setw(22) << " "; }
              out << fNilsson.fDecoupleLevels[isp]+1 << " : " ;
              for (int i=0; i<spindim; ++i) {
                out << std::fixed << std::setw(6) << std::setprecision(3) << wf[isp*spindim+i] << "  ";
              }
              ++prct;
              out << std::endl;
            }
          }
        }
        out << "sum = " << sum << std::endl;
        out << "KDom = " << GetKDom(spin,ind) << "   SPDom =" << GetSPDom(spin,ind) << std::endl;
        for (int isp = 0; isp < fNilsson.fDecoupleLevels.size(); ++isp) {
          out << "sum("<<fNilsson.fDecoupleLevels[isp] << ") = " << GetSPSum(spin,ind,fNilsson.fDecoupleLevels[isp])<<"   ";
        }
        out << std::endl;
        for (int ik = 0; ik < spindim; ++ik) {
          int k = kkmax-ik*4;
          out << "sum("<<k << ") = " << GetKSum(spin,ind,k)<<"   ";
        }        
        out << std::endl;
        out << "Max amplitude = (" << GetMaxKSP(spin,ind).first << ","<< GetMaxKSP(spin,ind).second << ")" << std::endl;

        out << std::endl;
      }      
    }
  }

  double ParticleRotor::GetSPSum(int spin, int indx, int ni) {
    if (fGam) {
      for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
        if (ni != fNilsson.fDecoupleLevels[i]) { continue; }
        int spindim = (spin+1)/2;
        int nlevels = fNilsson.fDecoupleLevels.size();
        int matdim = spindim*nlevels;
        int kkmax = spin;
        if ((spin-1)%4 != 0) { kkmax = spin - 2; }

        double sum = 0.0;
        arma::vec wf = fCompositions[spin].col(indx-1);
        for (int j=0; j<spindim; ++j) { //loop through K
          sum += std::pow(wf[i*spindim+j],2);
        }
        return sum;
      }
      std::cerr << "Level " << ni << " has no decoupling factors calculated for it" << std::endl;
      exit(1);
    }
    else {
      for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
        if (ni != fNilsson.fDecoupleLevels[i]) { continue; }
        arma::vec wf = fCompositions[spin].col(indx-1);
        return std::pow(wf[i],2);
      }
      std::cerr << "Level " << ni << " has no decoupling factors calculated for it" << std::endl;
      exit(1);
    }        
  }


  double ParticleRotor::GetKSum(int spin, int indx, int k) {
    if (fGam) {
      int spindim = (spin+1)/2;
      int nlevels = fNilsson.fDecoupleLevels.size();
      int matdim = spindim*nlevels;
      int kkmax = spin;
      if ((spin-1)%4 != 0) { kkmax = spin - 2; }

      for (int i=0; i<spindim; ++i) {
        if (std::abs(k) != std::abs(kkmax-4*i)) { continue; }
        double sum = 0.0;
        arma::vec wf = fCompositions[spin].col(indx-1);
        for (int j=0; j<nlevels; ++j) { //loop through SP levels
          sum += std::pow(wf[j*spindim+i],2);
        }
        return sum;
      }
      std::cerr << "Level K=" << k << " has no decoupling factors calculated for it" << std::endl;
      exit(1);
    }
    else {
      double sum = 0;
      for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
        if (std::abs(k) != std::abs(fKs[i])) { continue; }
        arma::vec wf = fCompositions[spin].col(indx-1);
        sum += std::pow(wf[i], 2);
      }
      return sum;
    }        
  }
  
  int ParticleRotor::GetSPDom(int spin, int indx) {
    if (fNilsson.fDecoupleLevels.size() == 0) { return -1; }
    
    if (fGam) {
      double maxsp = 0;
      int imaxsp = 0;
      for (int isp = 0; isp < fNilsson.fDecoupleLevels.size(); ++isp) {
        double sum = GetSPSum(spin,indx,fNilsson.fDecoupleLevels[isp]);        
        if (sum > maxsp) {
          maxsp = sum;
          imaxsp = fNilsson.fDecoupleLevels[isp];
        }
      }        
      return imaxsp;
    }
    else {
    double maxcomp = 0.0;
    int maxind = -1;
    
    for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
      double comp = std::abs(fCompositions[spin](i,indx-1));
      if (comp > maxcomp) {
        maxcomp = comp;
        maxind = i;
      }
    }

    if (maxind == -1) { std::cerr << "Error! Max composition not found!" << std::endl; return -1; }

    return fNilsson.fDecoupleLevels[maxind];
    }
  }
  
  int ParticleRotor::GetKDom(int spin, int indx) {
    if (fNilsson.fDecoupleLevels.size() == 0) { return -1; }

    int spindim = (spin+1)/2;
    int nlevels = fNilsson.fDecoupleLevels.size();
    int matdim = spindim*nlevels;
    int kkmax = spin;
    if ((spin-1)%4 != 0) { kkmax = spin - 2; }

    double maxk = 0;
    int imaxk = 0;        
    for (int ik=0; ik < spindim; ++ik) {        
      double sum = GetKSum(spin,indx,kkmax-4*ik);
      if (sum > maxk) {
        maxk = sum;
        imaxk = kkmax-4*ik;
      }
    }        
    return imaxk;
  }

  std::pair<int,int> ParticleRotor::GetMaxKSP(int spin, int indx) {
    if (!fGam) {
      int sp = GetSPDom(spin,indx);
      for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
        if (sp == fNilsson.fDecoupleLevels[i]) {
          return {fKs[i], sp};
        }
      }
      std::cerr << "Error in finding max K/SP! " << std::endl;
      exit(0);      
    }
    else {
      int spindim = (spin+1)/2;
      int nlevels = fNilsson.fDecoupleLevels.size();
      int matdim = spindim*nlevels;
      int kkmax = spin;
      if ((spin-1)%4 != 0) { kkmax = spin - 2; }
      arma::vec wf = fCompositions[spin].col(indx-1);

      double maxwf = 0.0;
      int maxk = 0;
      int maxsp = 0;
      for (int isp = 0; isp < fNilsson.fDecoupleLevels.size(); ++isp) {
        for (int i=0; i<spindim; ++i) {
          if (std::abs(wf[isp*spindim+i]) > maxwf) {            
            maxwf = std::abs(wf[isp*spindim+i]);
            maxk = kkmax - i*4;
            maxsp = fNilsson.fDecoupleLevels[isp];
          }
        }
      }
      return {maxk,maxsp};
    }
  }

  void ParticleRotor::PrintAmplitudes(std::ostream &out, std::vector<std::pair<int, int> >states) {
    if (fNilsson.fDecoupleLevels.size() == 0) { return; }    

    out << "          Energy       ";
    for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
      out << std::setw(3) << fNilsson.fDecoupleLevels[i]+1 << std::setw(6) << " ";
    }
    out << std::endl;
    for (auto &state : states) {
      int spin = state.first;
      int indx = state.second;
      out << " " << std::setw(3) << spin << "/2_" << indx << "  " << 
        std::fixed << std::setprecision(3) << std::setw(6) << fEnergies[spin][indx-1] << std::setw(4) << " ";
      for (int i=0; i<fNilsson.fDecoupleLevels.size(); ++i) {
        if (!fGam) {
        out << std::fixed << std::setprecision(3) << std::setw(6) << fCompositions[spin](i, indx-1) << "   ";
        }
        else {
          out << std::fixed << std::setprecision(3) << std::setw(6) << GetSPSum(spin,indx,fNilsson.fDecoupleLevels[i]) << "   ";
        }
      }
      out << std::endl;               
    }
  }  
}
