#include <unistd.h>
#include <sstream>
#include <signal.h>
#include <iostream>
#include <iomanip>
#include <string.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <stdexcept>
#include <math.h>
#include <fstream>
#include <time.h>
#include <vector>
 
double pi = 3.14159265359;
 
int main(){
 
  using namespace std;
 
  // ----------------------------------- //
  cout << "Title    : HOD (Harmonic Oscillator Density calculation)"  << endl;
  cout << "Author   : G.Takayama"  << endl;
  cout << "Date     : 2023. 01. 26 Thu"  << endl;
  cout << "Contents : Harmonic oscillator density calculation "  << endl;
  cout << "Reference: "  << endl;
 
  double A, Z, a[2], r0[2], rho0[2], rho[2];
  string file;
 
  cout << "Mass of nuclide            : A  ?" << endl;
  cin >> A;
  cout << "Atomic number of nuclide   : Z  ?" << endl;
  cin >> Z;
  cout << "Proton's width parameter   : a  ?" << endl;
  cin >> a[0];
  cout << "Proton's radius parameter  : r0 ?" << endl;
  cin >> r0[0];
  cout << "Neutron's width parameter  : a  ?" << endl;
  cin >> a[1];
  cout << "Neutron's radius parameter : r0 ?" << endl;
  cin >> r0[1];
  cout << "output file's name : XXX.txt ?" << endl;
  cin >> file;
 
  double N = A-Z;
  double tempA[2] = {Z,N};
 
  for(int i=0;i<2;i++){
    rho0[i] = pow(r0[i]*sqrt(pi),3) * (1+1.5*a[i]);
    rho0[i] = tempA[i]/rho0[i];
  }
 
  string csname;
  csname = "./CSV/HOD/"+file+".csv";
  ofstream File(csname.c_str(),std::ios::app);
 
  cout << "Harmonic Oscillator Density calculation" << endl;
  File << "Harmonic Oscillator Density calculation" << endl;
  cout << "File name : " << file << endl;
  File << "File name : " << file << endl;
  cout << "Proton,          " <<          ",         " << "Neutron,         " << ",          " << endl;
  File << "Proton,          " <<          ",         " << "Neutron,         " << ",          " << endl;
  cout << "radius parameter," << r0[0] << ",         " << "radius parameter," << r0[1] << ", " << endl;
  File << "radius parameter," << r0[0] << ",         " << "radius parameter," << r0[1] << ", " << endl;
  cout << "width parameter ," <<  a[0] << ",         " << "width parameter ," <<  a[1] << ", " << endl;
  File << "width parameter ," <<  a[0] << ",         " << "width parameter ," <<  a[1] << ", " << endl;
  cout << "r [fm]"    << ", " << "rho [/fm3] " << ", " << "r [fm]"    << ", " << "rho [/fm3] " << endl;
  File << "r [fm]"    << ", " << "rho [/fm3] " << ", " << "r [fm]"    << ", " << "rho [/fm3] " << endl;
 
  double r;
  double rSum[2];
  double rms[2];
 
  for(int j=0;j<300;j++){
    for(int i=0;i<2;i++){
      r = double(j)/10;
      rho[i] = rho0[i] * (1+a[i]*(r/r0[i])*(r/r0[i])) * exp(-(r/r0[i])*(r/r0[i]));
      rSum[i] += rho[i] * 4*pi*r*r*0.1;
      rms[i]  += rho[i] * 4*pi*pow(r,4)*0.1/tempA[i];
      cout << r << ",     " << rho[i] << ", ";
      File << r << ",     " << rho[i] << ", ";
    }
    cout << endl;
    File << endl;
  }
 
  cout << rSum[0] << endl;
  cout << rSum[1] << endl;
 
  cout << pow(rms[0],0.5) << endl;
  cout << pow(rms[1],0.5) << endl;
 
}