diff --git a/Cleopatra/Isotope.h b/Cleopatra/Isotope.h index 49f388f..966ac3c 100644 --- a/Cleopatra/Isotope.h +++ b/Cleopatra/Isotope.h @@ -1,536 +1,536 @@ -/*********************************************************************** - * - * This is Isotope.h, To extract the isotope mass from massXX.txt - * - *------------------------------------------------------- - * created by Ryan (Tsz Leung) Tang, Nov-18, 2018 - * email: goluckyryan@gmail.com - * ********************************************************************/ - - -#ifndef ISOTOPE_H -#define ISOTOPE_H - -#include -#include -#include -#include -#include -#include -#include "constant.h" // amu -#include //atoi -#include -using namespace std; - -string massData="/opt/Ptolemy/Cleopatra/mass20.txt"; - -// about the mass**.txt -// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu -// mass excess uncertaintly -// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons -// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons - -class Isotope { -public: - int A, Z; - double Mass, MassError, BEA; - string Name, Symbol; - string dataSource; - - Isotope(){findHeliosPath();}; - Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; - Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; - - void SetIso(int a, int z); - void SetIsoByName(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 FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; - void FindMassByName(string name); // give Z, mass, massError, BEA; - - int TwoJ(int nShell); - 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; - - - void setFileLines(){ - fileStartLine = 37; - fileEndLine = 3594; - - lineMass050_099 = 466; - lineMass100_149 = 1160; - lineMass150_199 = 1994; - lineMass200 = 2774; - } - - char * heliosPath; - bool isFindOnce; - - void findHeliosPath(){ - dataSource = massData; -// heliosPath = getenv("HELIOSSYS"); -// if( heliosPath ){ -// dataSource = heliosPath; -// dataSource += "/analysis" + massData; -// }else{ -// dataSource = ".." + massData; -// } - } - - - - -}; - -void Isotope::SetIso(int a, int z){ - this->A = a; - this->Z = z; - FindMassByAZ(a,z); -} - -void Isotope::SetIsoByName(string name){ - FindMassByName(name); -} - -void Isotope::FindMassByAZ(int A, int Z){ - string line; - int lineNum=0; - int list_A, list_Z; - - ifstream myfile; - int flag=0; - - setFileLines(); - - 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 ) numLineStart = lineMass200; - - myfile.open(dataSource.c_str()); - - if (myfile.is_open()) { - while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ - list_Z = atoi((line.substr(10,5)).c_str()); - list_A = atoi((line.substr(15,5)).c_str()); - - if ( A == list_A && Z == list_Z) { - this->BEA = atof((line.substr(54,11)).c_str()); - this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A; - this->MassError = atof((line.substr(65,7)).c_str()); - string str = line.substr(20,2); - str.erase(remove(str.begin(), str.end(), ' '), str.end()); - this->Symbol = str; - - ostringstream ss; - ss << A << this->Symbol; - this->Name = ss.str(); - flag = 1; - }else if ( list_A > A) { - this->BEA = -404; - this->Mass = -404; - this->MassError = -404; - this->Symbol = "non"; - this->Name = "non"; - break; - } - - } - } - - 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"; - - myfile.close(); - }else { - printf("Unable to open %s\n", dataSource.c_str()); - } -} - -void Isotope::FindMassByName(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; - return; - } - if( name == "p" ) name = "1H"; - if( name == "d" ) name = "2H"; - if( name == "t" ) name = "3H"; - if( name == "a" ) name = "4He"; - - string temp = name; - int lastDigit = 0; - - for(int i=0; temp[i]; i++){ - if(temp[i] == '0') lastDigit = i; - if(temp[i] == '1') lastDigit = i; - if(temp[i] == '2') lastDigit = i; - if(temp[i] == '3') lastDigit = i; - if(temp[i] == '4') lastDigit = i; - if(temp[i] == '5') lastDigit = i; - if(temp[i] == '6') lastDigit = i; - if(temp[i] == '7') lastDigit = i; - if(temp[i] == '8') lastDigit = i; - if(temp[i] == '9') lastDigit = i; - } - - this->Symbol = temp.erase(0, lastDigit +1); - //check is Symbol is 2 charaters, if not, add " " at the end - if( this->Symbol.length() == 1 ){ - this->Symbol = this->Symbol + " "; - } - - - temp = name; - int len = temp.length(); - temp = temp.erase(lastDigit+1, len); - - this->A = atoi(temp.c_str()); - //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); - - // find the nucleus in the data - string line; - int lineNum=0; - int list_A; - string list_symbol; - - ifstream myfile; - int flag=0; - - setFileLines(); - - 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 ) numLineStart = lineMass200; - - myfile.open(dataSource.c_str()); - - if (myfile.is_open()) { - while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ - list_symbol = line.substr(20,2); - list_A = atoi((line.substr(15,5)).c_str()); - - //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str()); - - if ( this->A == list_A && this->Symbol == list_symbol) { - this->Z = atoi((line.substr(10,5)).c_str()); - this->BEA = atof((line.substr(54,11)).c_str()); - this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A; - this->MassError = atof((line.substr(65,7)).c_str()); - - string str = line.substr(20,2); - str.erase(remove(str.begin(), str.end(), ' '), str.end()); - this->Symbol = str; - - ostringstream ss; - ss << this->A << this->Symbol; - this->Name = ss.str(); - flag = 1; - }else if ( list_A > this->A) { - this->BEA = -404; - this->Mass = -404; - this->MassError = -404; - this->Symbol = "non"; - this->Name = "non"; - 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( nucleusD.Mass != -404){ - return nucleusD.Mass + Nn*mn + Np*mp - this->Mass; - }else{ - return -404; - } -} - -double Isotope::CalSp2(int a, int z){ - Isotope nucleusD(A - a , Z - z); - Isotope nucleusS(a,z); - - if( nucleusD.Mass != -404 && nucleusS.Mass != -404){ - return nucleusD.Mass + nucleusS.Mass - this->Mass; - }else{ - return -404; - } -} - -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; -} - -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 = 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 = 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){ - - findHeliosPath(); - - printf(" using mass data : %s \n", dataSource.c_str()); - 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, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.); - printf(" total BE : %12.5f MeV\n",BEA*A/1000.); - printf(" mass in amu : %12.5f u\n",Mass/amu); - printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu); - 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 +/*********************************************************************** + * + * This is Isotope.h, To extract the isotope mass from massXX.txt + * + *------------------------------------------------------- + * created by Ryan (Tsz Leung) Tang, Nov-18, 2018 + * email: goluckyryan@gmail.com + * ********************************************************************/ + + +#ifndef ISOTOPE_H +#define ISOTOPE_H + +#include +#include +#include +#include +#include +#include +#include "constant.h" // amu +#include //atoi +#include +using namespace std; + +string massData="../Cleopatra/mass20.txt"; + +// about the mass**.txt +// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu +// mass excess uncertaintly +// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons +// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons + +class Isotope { +public: + int A, Z; + double Mass, MassError, BEA; + string Name, Symbol; + string dataSource; + + Isotope(){findHeliosPath();}; + Isotope(int a, int z){ findHeliosPath();SetIso(a,z); }; + Isotope(string name){ findHeliosPath(); SetIsoByName(name); }; + + void SetIso(int a, int z); + void SetIsoByName(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 FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol; + void FindMassByName(string name); // give Z, mass, massError, BEA; + + int TwoJ(int nShell); + 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; + + + void setFileLines(){ + fileStartLine = 37; + fileEndLine = 3594; + + lineMass050_099 = 466; + lineMass100_149 = 1160; + lineMass150_199 = 1994; + lineMass200 = 2774; + } + + char * heliosPath; + bool isFindOnce; + + void findHeliosPath(){ + dataSource = massData; +// heliosPath = getenv("HELIOSSYS"); +// if( heliosPath ){ +// dataSource = heliosPath; +// dataSource += "/analysis" + massData; +// }else{ +// dataSource = ".." + massData; +// } + } + + + + +}; + +void Isotope::SetIso(int a, int z){ + this->A = a; + this->Z = z; + FindMassByAZ(a,z); +} + +void Isotope::SetIsoByName(string name){ + FindMassByName(name); +} + +void Isotope::FindMassByAZ(int A, int Z){ + string line; + int lineNum=0; + int list_A, list_Z; + + ifstream myfile; + int flag=0; + + setFileLines(); + + 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 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_Z = atoi((line.substr(10,5)).c_str()); + list_A = atoi((line.substr(15,5)).c_str()); + + if ( A == list_A && Z == list_Z) { + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + break; + } + + } + } + + 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"; + + myfile.close(); + }else { + printf("Unable to open %s\n", dataSource.c_str()); + } +} + +void Isotope::FindMassByName(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; + return; + } + if( name == "p" ) name = "1H"; + if( name == "d" ) name = "2H"; + if( name == "t" ) name = "3H"; + if( name == "a" ) name = "4He"; + + string temp = name; + int lastDigit = 0; + + for(int i=0; temp[i]; i++){ + if(temp[i] == '0') lastDigit = i; + if(temp[i] == '1') lastDigit = i; + if(temp[i] == '2') lastDigit = i; + if(temp[i] == '3') lastDigit = i; + if(temp[i] == '4') lastDigit = i; + if(temp[i] == '5') lastDigit = i; + if(temp[i] == '6') lastDigit = i; + if(temp[i] == '7') lastDigit = i; + if(temp[i] == '8') lastDigit = i; + if(temp[i] == '9') lastDigit = i; + } + + this->Symbol = temp.erase(0, lastDigit +1); + //check is Symbol is 2 charaters, if not, add " " at the end + if( this->Symbol.length() == 1 ){ + this->Symbol = this->Symbol + " "; + } + + + temp = name; + int len = temp.length(); + temp = temp.erase(lastDigit+1, len); + + this->A = atoi(temp.c_str()); + //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A); + + // find the nucleus in the data + string line; + int lineNum=0; + int list_A; + string list_symbol; + + ifstream myfile; + int flag=0; + + setFileLines(); + + 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 ) numLineStart = lineMass200; + + myfile.open(dataSource.c_str()); + + if (myfile.is_open()) { + while (/*! myfile.eof() &&*/ flag == 0 && lineNum = numLineStart ){ + list_symbol = line.substr(20,2); + list_A = atoi((line.substr(15,5)).c_str()); + + //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str()); + + if ( this->A == list_A && this->Symbol == list_symbol) { + this->Z = atoi((line.substr(10,5)).c_str()); + this->BEA = atof((line.substr(54,11)).c_str()); + this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A; + this->MassError = atof((line.substr(65,7)).c_str()); + + string str = line.substr(20,2); + str.erase(remove(str.begin(), str.end(), ' '), str.end()); + this->Symbol = str; + + ostringstream ss; + ss << this->A << this->Symbol; + this->Name = ss.str(); + flag = 1; + }else if ( list_A > this->A) { + this->BEA = -404; + this->Mass = -404; + this->MassError = -404; + this->Symbol = "non"; + this->Name = "non"; + 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( nucleusD.Mass != -404){ + return nucleusD.Mass + Nn*mn + Np*mp - this->Mass; + }else{ + return -404; + } +} + +double Isotope::CalSp2(int a, int z){ + Isotope nucleusD(A - a , Z - z); + Isotope nucleusS(a,z); + + if( nucleusD.Mass != -404 && nucleusS.Mass != -404){ + return nucleusD.Mass + nucleusS.Mass - this->Mass; + }else{ + return -404; + } +} + +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; +} + +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 = 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 = 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){ + + findHeliosPath(); + + printf(" using mass data : %s \n", dataSource.c_str()); + 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, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.); + printf(" total BE : %12.5f MeV\n",BEA*A/1000.); + printf(" mass in amu : %12.5f u\n",Mass/amu); + printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu); + 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