/***********************************************************************
 * 
 *  This is ClassIsotope.h, To extract the isotope mass from nubase_4.mas20.txt
 * 
 *-------------------------------------------------------
 *  created by Ryan (Tsz Leung) Tang, Sept-6, 2024
 *  email: goluckyryan@gmail.com
 * ********************************************************************/


#ifndef ISOTOPE_H
#define ISOTOPE_H

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include "constant.h" // amu
#include <stdlib.h>  //atoi
#include <algorithm>

std::string massData="nubase_4.mas20.txt";

class Isotope {
public:
  int A, Z;
  double Mass, MassError, BEA; // BEA in keV
  short twoSpin; // 2*J
  short parity;
  std::string Name, Symbol; // Name = 6Li, Symbol = Li
  std::string dataSource;
  
  Isotope(){dataSource = massData;};
  Isotope(int a, int z){ dataSource = massData; SetIso(a,z);  };
  Isotope(std::string name){ dataSource = massData; SetIsoByName(name); };

  void SetMassTablePath(std::string path){ dataSource = path;}

  void SetIso(int a, int z);
  void SetIsoByName(std::string name);

  double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal 
  double CalSp2(int a, int z); // this is for (a,z) nucleus removal

  double CalBeta(double T){
    double Etot = Mass + T;
    double gamma = 1 + T/Mass;
    double beta = sqrt(1 - 1 / gamma / gamma ) ;
    return beta;
  }

  void Print();
  void ListShell();
   
private:

  void InvalidIso();
  void BreakName(std::string nameStr);
  void BreakJpi(std::string jPiStr);
  void DigestLine(std::string line, int A, int Z);

  void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol;
  void FindMassByName(std::string name); // give Z, mass, massError, BEA;

  int TwoJ(int nShell);
  std::string Orbital(int nShell);
  int magic(int i){
    switch (i){
      case 0: return   2; break;
      case 1: return   8; break;
      case 2: return  20; break;
      case 3: return  28; break;
      case 4: return  40; break;
      case 5: return  50; break;
      case 6: return  82; break;
      case 7: return 128; break;
    }
    return 0;
  }

  int magicShellID(int i){
    switch (i){
      case 0: return   0; break;
      case 1: return   2; break;
      case 2: return   5; break;
      case 3: return   6; break;
      case 4: return   9; break;
      case 5: return  10; break;
      case 6: return  15; break;
      case 7: return  21; break;
    }
    return 0;
  }
  
  int fileStartLine;
  int fileEndLine;
  int lineMass050_099;
  int lineMass100_149;
  int lineMass150_199;
  int lineMass200_249;
  int lineMass250;

  int colA[2];
  int colZ[2];
  int colName[2];
  int colMassExcess[2];
  int colMassError[2];
  int colSpin[2];
  
  void setFileLinesAndCol(){
    fileStartLine = 26;
    fileEndLine = 5869;
    
    lineMass050_099 = 638;
    lineMass100_149 = 1718;
    lineMass150_199 = 3123;
    lineMass200_249 = 4591;
    lineMass250     = 5495;

    colA[0] = 0;
    colA[1] = 3;

    colZ[0] = 4;
    colZ[1] = 8;

    colName[0] = 11;
    colName[1] = 16;

    colMassExcess[0] = 18;
    colMassExcess[1] = 31;

    colMassError[0] = 31;
    colMassError[1] = 42;

    colSpin[0] = 88;
    colSpin[1] = 97;
  }
  
  bool isFindOnce;
  
};

void Isotope::SetIso(int a, int z){
  FindMassByAZ(a,z); 
}

void Isotope::SetIsoByName(std::string name){
  FindMassByName(name); 
}

void Isotope::InvalidIso(){
  this->BEA  = std::nan("");
  this->Mass = std::nan("");
  this->MassError = std::nan("");
  this->Symbol = "non";
  this->Name   = "non";
  this->twoSpin = 0;
  this->parity = 0;
}

void Isotope::BreakName(std::string nameStr){
  std::string::size_type lastDigit = nameStr.find_last_of("0123456789");

  this->Symbol = nameStr.substr(lastDigit + 1);
  if (this->Symbol.length() == 1) {
    this->Symbol += " ";
  }

  this->A = std::stoi(nameStr.substr(0, lastDigit + 1));

}

void Isotope::BreakJpi(std::string jPiStr){

  jPiStr.erase(std::remove(jPiStr.begin(), jPiStr.end(), '('), jPiStr.end());
  jPiStr.erase(std::remove(jPiStr.begin(), jPiStr.end(), ')'), jPiStr.end());

  std::string::size_type lastDigit = jPiStr.find_last_of("0123456789");

  std::string spinStr = jPiStr.substr(0, lastDigit + 1);

  std::string::size_type lastSlash = spinStr.find_last_of('/');

  if (lastSlash != std::string::npos) {
    twoSpin = std::stoi(spinStr.substr(0, lastSlash + 1));
  } else {
    twoSpin = 2 * std::stoi(spinStr);
  }

  this->parity = (jPiStr.substr(lastDigit + 1) == "+") ? 1 : -1;

}

void Isotope::DigestLine(std::string line, int A, int Z){

  double massExcess = atof((line.substr(colMassExcess[0],colMassExcess[1]-colMassExcess[0])).c_str());

  Mass = massExcess/1000. + A * amu - Z * me;
  MassError =  atof((line.substr(colMassError[0],colMassError[1]-colMassError[0])).c_str()) / 1000.;

  BEA = (Z * mp + (A-Z) * mn - Mass)/A *1000;

  std::string str = line.substr(colName[0],colName[1]-colName[0]);
  str.erase(remove(str.begin(), str.end(), ' '), str.end());
  Name = str;

  BreakName(Name);

  if( this->Name == "1H" ) this->Name = "p";
  if( this->Name == "2H" ) this->Name = "d";
  if( this->Name == "3H" ) this->Name = "t";
  if( this->Name == "4He" ) this->Name = "a";

  str = line.substr(colSpin[0], colSpin[1]-colSpin[0]);
  str.erase(remove(str.begin(), str.end(), ' '), str.end());
  str.erase(std::remove(str.begin(), str.end(), '*'), str.end());
  str.erase(std::remove(str.begin(), str.end(), '#'), str.end());
  std::string Jpi = str;

  BreakJpi(Jpi);

  // printf("A %d  Z: %d | %s | %f | %f |%s|\n", A, Z, Name.c_str(), massExcess, Mass, Jpi.c_str());
}


void Isotope::FindMassByAZ(int A, int Z){
  std::string line;
  int    lineNum=0;

  std::ifstream myfile;
  int    flag=0;

  setFileLinesAndCol();

  int numLineStart = fileStartLine;
  int numLineEnd  = fileEndLine;

  this->A = A;
  this->Z = Z;

  if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
  if ( A >=100 && A < 150) numLineStart = lineMass100_149;
  if ( A >=150 && A < 200) numLineStart = lineMass150_199;
  if ( A >=200 && A < 250) numLineStart = lineMass200_249;
  if ( A >=250           ) numLineStart = lineMass250;

  myfile.open(dataSource.c_str());

  if (myfile.is_open()) {
    while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
      lineNum ++ ;
      // printf("%3d  ",lineNum);
      getline (myfile,line);

      if (lineNum >= numLineStart ){
        int list_Z     = std::stoi((line.substr(colZ[0],colZ[1]-colZ[0])));
      	int list_A     = std::stoi((line.substr(colA[0],colA[1]-colA[0])));

      	if ( A == list_A && Z * 10 == list_Z) {
    
          DigestLine(line, A, Z);          

     		  flag = 1;
      	}else if ( list_A > A) {
          InvalidIso();
          break;
        }

      }
    }
    
    myfile.close();
  }else {
    printf("Unable to open %s\n", dataSource.c_str());
  }
}

void Isotope::FindMassByName(std::string name){

  // done seperate the Mass number and the name 
  if( name == "n" ) {
    this->Name = "1n";
    this->BEA       = 0;
    this->Mass      = mn;
    this->MassError = 0;
    this->Name      = "n";
    this->A         = 1;
    this->Z         = 0;
    this->twoSpin   = 1;
    this->parity    = 1;
    return;
  }

  if( name == "p" ) name = "1H";
  if( name == "d" ) name = "2H";
  if( name == "t" ) name = "3H";
  if( name == "a" ) name = "4He";

  this->Name = name;
  BreakName(name);

  // find the nucleus in the data
  std::string line;
  int    lineNum=0;

  std::ifstream myfile;
  int    flag=0;

  setFileLinesAndCol();

  int numLineStart = fileStartLine;
  int numLineEnd  = fileEndLine;

  if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
  if ( A >=100 && A < 150) numLineStart = lineMass100_149;
  if ( A >=150 && A < 200) numLineStart = lineMass150_199;
  if ( A >=200 && A < 250) numLineStart = lineMass200_249;
  if ( A >=250           ) numLineStart = lineMass250;

  myfile.open(dataSource.c_str());

  if (myfile.is_open()) {
    while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
      lineNum ++ ;
      //printf("%3d  ",lineNum);
      getline (myfile,line);

      if (lineNum >= numLineStart ){

        std::string nameStr = line.substr(colName[0],colName[1]-colName[0]);
        nameStr.erase(remove(nameStr.begin(), nameStr.end(), ' '), nameStr.end());

        int list_A     = std::stoi((line.substr(colA[0],colA[1]-colA[0])));
        int list_Z     = std::stoi((line.substr(colZ[0],colZ[1]-colZ[0])));

        // printf("A %d Z %d | %s \n", list_A, list_Z, nameStr.c_str());

        if ( list_A == A && list_Z % 10 == 0  &&  Name == nameStr) {

          Z = list_Z/10;

          DigestLine(line, list_A, Z);

          flag = 1;
        }else if ( list_A > this->A) {
          InvalidIso();
          break;
        }
      }
    }
    myfile.close();
  }else {
    printf("Unable to open %s\n", dataSource.c_str());
  }
}

double Isotope::CalSp(int Np, int Nn){
  Isotope nucleusD(A - Np - Nn, Z - Np);  

  if( !std::isnan(nucleusD.Mass) ){
    return nucleusD.Mass + Nn*mn + Np*mp - this->Mass;
  }else{
    return std::nan("");
  }
}

double Isotope::CalSp2(int a, int z){
  Isotope nucleusD(A - a , Z - z);
  Isotope nucleusS(a,z);  

  if( !std::isnan(nucleusD.Mass) && !std::isnan(nucleusS.Mass) ){
    return nucleusD.Mass + nucleusS.Mass - this->Mass;
  }else{
    return std::nan("");
  }
}

int Isotope::TwoJ(int nShell){

  switch(nShell){
    case  0: return  1; break; // 0s1/2
    case  1: return  3; break; // 0p3/2
    case  2: return  1; break; // 0p1/2 -- 8
    case  3: return  5; break; // 0d5/2
    case  4: return  1; break; // 1s1/2
    case  5: return  3; break; // 0d3/2 -- 20
    case  6: return  7; break; // 0f7/2 -- 28
    case  7: return  3; break; // 1p3/2
    case  8: return  1; break; // 1p1/2
    case  9: return  5; break; // 0f5/2 -- 40
    case 10: return  9; break; // 0g9/2 -- 50
    case 11: return  7; break; // 0g7/2
    case 12: return  5; break; // 1d5/2
    case 13: return 11; break; // 0h11/2
    case 14: return  3; break; // 1d3/2
    case 15: return  1; break; // 2s1/2 -- 82
    case 16: return  9; break; // 0h9/2
    case 17: return  7; break; // 1f7/2
    case 18: return 13; break; // 0i13/2
    case 19: return  3; break; // 2p3/2
    case 20: return  5; break; // 1f5/2
    case 21: return  1; break; // 1p1/2 -- 126
    case 22: return  9; break; // 1g9/2
    case 23: return 11; break; // 0i11/2
    case 24: return 15; break; // 0j15/2
    case 25: return  5; break; // 2d5/2
    case 26: return  1; break; // 3s1/2
    case 27: return  3; break; // 2d3/2
    case 28: return  7; break; // 1g7/2
  }

  return 0;
}

std::string Isotope::Orbital(int nShell){

  switch(nShell){
    case  0: return  "0s1 "; break;  //
    case  1: return  "0p3 "; break;  //
    case  2: return  "0p1 "; break;  //-- 8
    case  3: return  "0d5 "; break;  //
    case  4: return  "1s1 "; break;  //
    case  5: return  "0d3 "; break;  //-- 20
    case  6: return  "0f7 "; break;  //-- 28
    case  7: return  "1p3 "; break;  //
    case  8: return  "1p1 "; break;  //
    case  9: return  "0f5 "; break;  //-- 40
    case 10: return  "0g9 "; break;  //-- 50
    case 11: return  "0g7 "; break;  //
    case 12: return  "1d5 "; break;  //
    case 13: return  "0h11"; break;  //
    case 14: return  "1d3 "; break;  //
    case 15: return  "2s1 "; break;  //-- 82
    case 16: return  "0h9 "; break;  //
    case 17: return  "1f7 "; break;  //
    case 18: return  "0i13"; break;  //
    case 19: return  "2p3 "; break;  //
    case 20: return  "1f5 "; break;  //
    case 21: return  "1p1 "; break;  //-- 126
    case 22: return  "1g9 "; break;  //
    case 23: return  "0i11"; break;  //
    case 24: return  "0j15"; break;  //
    case 25: return  "2d5 "; break;  //
    case 26: return  "3s1 "; break;  //
    case 27: return  "2d3 "; break;  //
    case 28: return  "1g7 "; break;  //
  }

  return "nan";
}

void Isotope::ListShell(){

  if( Mass < 0 ) return;

  int n = A-Z;
  int p = Z;

  int k = std::min(n,p);
  int nMagic = 0;
  for( int i = 0; i < 7; i++){
    if( magic(i) < k && k <= magic(i+1) ){
      nMagic = i;
      break;
    }
  }

  int coreShell = magicShellID(nMagic-1);
  int coreZ1 = magic(nMagic-1);
  int coreZ2 = magic(nMagic);

  Isotope core1( 2*coreZ1, coreZ1);
  Isotope core2( 2*coreZ2, coreZ2);

  printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str());
  printf("         || ");
  int t = std::max(n,p);
  int nShell = 0;
  do{
    int occ = TwoJ(nShell)+1;
    if( nShell > coreShell) {
      printf("%4s", Orbital(nShell).c_str());
         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
        printf("|");
      }else{
        printf(",");
      } 
    }
    t = t - occ;
    nShell++;
  }while( t > 0  && nShell < 29);
  for( int i = 1; i <= 6; i++){
    if (nShell < 28) {
      printf("%4s,", Orbital(nShell).c_str());
    }else if( nShell == 28 ) {
      printf("%4s", Orbital(nShell).c_str());
    }
    nShell ++;
  }
  if (nShell < 29) printf("%4s", Orbital(nShell).c_str());
  printf("\n");


  printf(" Z = %3d || ", p);
  nShell = 0;
  do{
    int occ = TwoJ(nShell)+1;
    if( nShell > coreShell ){
      if( p > occ ) {
         printf("%-4d", occ);
         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
          printf("|");
        }else{
          printf(",");
        } 
      }else{
        printf("%-4d", p);
      }
    }
    p = p - occ;
    nShell++;
  }while( p > 0  && nShell < 29);
  printf("\n");

  printf(" N = %3d || ", n);
  nShell = 0;  
  do{
    int occ = TwoJ(nShell)+1;
    if ( nShell > coreShell ){
      if( n > occ ) {
         printf("%-4d", occ);
         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
          printf("|");
        }else{
          printf(",");
        } 
      }else{
        printf("%-4d", n);
      }
    }
    n = n - occ;
    nShell++;
  }while( n > 0  && nShell < 29);
  printf("\n");

  printf("------------------ \n");
}



void Isotope::Print(){

  if (Mass > 0){  
    
    printf(" using mass data : %s \n", dataSource.c_str());

    std::string jPiStr = "";
    if( twoSpin % 2  == 0) {
      jPiStr = std::to_string(twoSpin/2).c_str();
    }else{
      jPiStr = std::to_string(twoSpin) + "/2";
    }

    printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, J-pi \e[31m%s%s\e[m\n", 
            Name.c_str(), 
            Z, 
            A, 
            Mass,
            jPiStr.c_str(),
            parity == 1 ? "+" : "-");     

    printf(" total BE    : %12.5f MeV\n",BEA*A/1000.);
    printf(" total BE/A  : %12.5f MeV\n",BEA/1000.);
    printf(" mass in amu : %12.5f u\n",Mass/amu);
    printf(" mass excess : %12.5f MeV +- %12.5f\n", Mass + Z*0.510998950 - A*amu, MassError);
    printf("-------------- Seperation energy \n");
    printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1));
    printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1));
    printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n",               CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2));
    printf(" S4p: %8.4f| S4n: %8.4f| \n",                             CalSp(4, 0), CalSp(0, 4));

  }else{
    printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A);
  }


}

#endif