////////////////////////////////////////////////////////////////////////////
//
//    This  program  is  a  C/C++   code   provided  as  an   accessory   to
//    interrogate  the  provided datasets,  which  consist  of emission  and
//    recombination  coefficients of  hydrogen  calculated for  Case B  with
//    electron  energies described  by  kappa distributions.   The data  are
//    tabulated  on  a  grid   of  electron  number  density,  Ne,  electron
//    temperature, Te, and kappa.   If necessary the program interpolates to
//    the  user  defined  values  of  the three  parameters  using  Lagrange
//    interpolation.  Interpolated  values  are  given  for  polynomials  of
//    orders: 1, 3 and 5. The  purpose of multiple interpolations is to test
//    the convergence and hence check the reliability of the interpolation.
//
//    Two main input data files are  required to run the program: e1bk.d and
//    t1bk.d.  The first file contains emission coefficients for transitions
//    between states with  upper quantum number nu and  lower quantum number
//    nl,  where the quantum  numbers range  between 2-99  with nl<nu,  as a
//    function of Ne, Te and  kappa. The second file contains total hydrogen
//    recombination  coefficients and recombination  coefficients to  the 2s
//    state as a  function of Ne, Te and kappa. The  structure of these data
//    files  is fully  explained in  the  documentation of  the dataset  and
//    associated paper (see ReadMe.dat file  and [1] in the References; also
//    refer to [2]).
//
//    During the  execution of the program,  options are given  to guide the
//    user and obtain the required values.  All input and output data are in
//    the cgs system.  There are three main options:
//
//    1)  Obtaining the emission coefficient in  units  of   erg.cm^3.s^{-1}
//    for a given transition. The required inputs are the  upper  and  lower
//    quantum numbers  of  the transition, as well as Ne, Te and kappa.
//
//    2) Obtaining the total recombination  coefficient of H for a given set
//    of physical parameters. The required inputs are Ne, Te and kappa.
//
//    3) Obtaining the total recombination coefficient to the 2s state of  H
//    for a given set of physical parameters. The required inputs are Ne, Te
//    and kappa.
//
//    As indicated already, the raw  and processed data are calculated based
//    on the  assumption of kappa  electron energy distribution  rather than
//    Maxwell-Boltzmann   (MB)   distribution.    However,  data   with   MB
//    distribution can also be obtained  by using large values of kappa. For
//    this purpose, the  use of the maximum available  kappa value, which is
//    10^6, is recommended.
//
//    ====================================================
//
//    Program data sheet:
//    -Name: intrat3D
//    -Version: 1.0
//    -Date: 2014
//    -Language: C/C++
//    -Tested compilers: g++, Dev-C++ 5.7.1, Microsoft Visual Studio 6.0
//    -Compilation from command line (e.g. for g++): $g++ intrat3D.cpp
//    -Tested platforms: Ubuntu 12.04, Windows XP and Windows 8 (64-bit)
//
//    ====================================================
//
//    Acknowledgement:
//
//    -P.J. Storey: University College London, pjs@star.ucl.ac.uk
//    -Taha Sochi:  University College London, t.sochi@ucl.ac.uk
//
//    ====================================================
//
//    Notes:
//
//    (1) The  intrat code and documentation  of Hummer and  Storey [2] were
//    consulted  and  used  as  a   model  for  writing  intrat3D  code  and
//    documentation.
//
//    (2)  Error  reports,  suggestions  and  enquiries should  be  sent  to
//    P.J. Storey or Taha Sochi.
//
//    (3) The users of this  code  should  avoid using  the double precision
//    scientific  d-notation  which is used in  fortran in their  input. For
//    instance, to input 1000 in  scientific notation they  should  use  1e3
//    rather than 1d3 to avoid failure.
//
//    ====================================================
//
//    References:
//
//    [1]  P.J.  Storey,  Taha   Sochi  (2014)  Emission  and  recombination
//    coefficients for  hydrogen with kappa  electron  energy  distribution.
//    MNRAS X(X): X-X.
//
//    [2] P.J.   Storey, D.G.  Hummer (1995)  Recombination line intensities
//    for  hydrogenic  ions  -  IV.  Total  recombination  coefficients  and
//    machine-readable tables for Z=1 to 8. MNRAS 272(1): 41-48.
//
////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <cmath>

using namespace std;

const int mxl = 5000;
const int mxn = 20;
const int mxt = 20;
const int mxk = 50;

struct _theory
{
   double densx[mxn];
   double tempx[mxt];
   double kappax[mxk];
   double skap[mxk];
   double ex[mxl][mxn][mxt][mxk];
   double a[mxn][mxt][mxk];
   double a2s[mxn][mxt][mxk];
   int ntempx, ndensx, nkappax, nmin, ncut;
} theory;

double phi(double *xiq, int lg, int itau, double x)
{
   int l;
   double fk = 1.0;
   for (l = 1; l <= lg; ++l)
      if (l != itau)
         fk *= x - xiq[l-1];
   return fk;
}

void lagr(double *f, double *x, double *fv, double xv, int np, int ni, int ndim, int *ier)
{
   double xiq[10];
   int lg = 2*np;
   *ier = 0;
   int i,in,ilow,itau;

   if (xv < x[0] || xv > x[ni-1])
   {
      *ier = 1;
      return;
   }

   //  LOCATE XV IN X
   for (i = 2; i <= ni; ++i)
   {
      if (x[i-1] >= xv)
      {
         //  SELECT SUBSET OF X FOR INTERPOLATION, TAKING CARE AT ENDS
         in = ni - lg + 1;
         ilow = i - np;

         if (ilow < 1)
            ilow = 1;

         if (ilow > in)
            ilow = in;

         for ( itau = 1; itau <= lg; ++itau)
            xiq[itau-1] = x[ilow + itau - 2];

         //  CONSTRUCT INTERPOLATION COEFFICIENTS
         *fv = 0.0;
         for ( itau = 1; itau <= lg; ++itau)
         {
            double fl = phi(xiq, lg, itau, xv) / phi(xiq, lg, itau, xiq[itau - 1]);
            *fv += fl * f[ilow + itau - 2];
         }
         return;
      }
   }
}

double gammln(double xx)
{
   static const double cof[6] = { 76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155,.001208650973866179, -5.395239384953e-6 };
   static const double stp = 2.5066282746310005;
   double x, y, ser, tmp;
   int j;

   x = xx;
   y = x;
   tmp = x + 5.5;
   tmp = (x + 0.5) * log(tmp) - tmp;
   ser = 1.000000000190015;

   for ( j = 0; j < 6; ++j)
   {
      y += 1.0;
      ser += cof[j] / y;
   }

   return tmp + log(stp * ser / x);
}

double onedk(int il, int in, int jt, double zk, int np, int iopt)
{
   double fk[mxk];
   double xk[mxk];
   int ik;

   for ( ik = 0; ik < theory.nkappax; ++ik)
   {
      xk[ik] = theory.kappax[ik];
      if (iopt == 1)
         fk[ik] = theory.ex[il-1][in-1][jt-1][ik] * theory.skap[ik];

      else if (iopt == 2)
         fk[ik] = theory.a[in-1][jt-1][ik] *	theory.skap[ik];

      else if (iopt == 3)
         fk[ik] = theory.a2s[in-1][jt-1][ik] *theory.skap[ik];

      fk[ik] = log10(fk[ik]);
   }

   int ier;
   double fvk;
   double zkl = log10(zk);
   lagr(fk,xk,&fvk,zkl,np,theory.nkappax,mxk,&ier);

   return fvk;
}

double onedt(int il, int in, double yt, double zk, int np, int iopt)
{
   double ft[mxt];
   double xt[mxt];

   for (int it = 1; it <= theory.ntempx; ++it)
   {
      xt[it-1] = theory.tempx[it-1];
      ft[it-1] = onedk(il, in, it, zk, np, iopt);
   }

   double ytl = log10(yt);
   int ier;
   double fvt;

   lagr(ft,xt,&fvt,ytl,np,theory.ntempx,mxt,&ier);

   return fvt;
}

double onedn(int il, double xn, double yt, double zk, int np, int iopt)
{
   double an[mxn], fn[mxn];

   for (int in = 1; in <= theory.ndensx; ++in)
   {
      fn[in - 1] = onedt(il,in,yt,zk,np,iopt);
      an[in - 1] = theory.densx[in - 1];
   }

   int  ier;
   double fvn;
   double xnl = log10(xn);

   lagr(fn,an,&fvn,xnl,np,theory.ndensx,mxn,&ier);

   return fvn;
}


double threedntk(int il, double xn, double yt, double zk, int np, int iopt)
{
   double t = pow(10.0, onedn(il, xn, yt, zk, np, iopt) );

   // remove kappa scaling
   double gg = 1.5*log(zk - 1.5) + gammln(zk - 0.5) - gammln(zk + 1.0);
   gg = exp(gg) / (1.0 - 3.0 / (zk * 2.0));

   return t/gg;
}


void getinput(void)
{
   static const int ni[3] = {1, 2, 3};
   static const int maxx = 3;
   static const char opts[5] = {'e', 'i', 'n', 'c', 'r'};
   int il,i,k;
   double ri[3];
   double zk, xn, yt;
   double lower, upper;
   int iopt, istop;

   static const string emis_coef = "Emission coefficient";
   static const string re_coef_t = "Recombination coefficient (total)";
   static const string re_coef_s = "Recombination coefficient (2s)";
   static const string dens = "    density";
   static const string temp = "temperature";
   static const string kppa = "      kappa";
   static const string inor = "interp. order =";

   ofstream out16("intrat3D.d" );

   cout << " Input and units (in cgs):\n"
        << " Electron number density: cm^{-3}\n"
        << " Electron temperature: K\n"
        << " Kappa: no units\n\n"
        << " Input options (iopt)\n"
        << " iopt = 1 for emission coefficient for a specific transition\n"
        << " iopt = 2 for total recombination coefficient for H in Case B\n"
        << " iopt = 3 for recombination coefficient to the H 2s state\n\n"
        << " Output and units (in cgs):\n"
        << " iopt = 1  emission coefficient in erg.cm^3.s^{-1}\n"
        << " iopt = 2 and 3 recombination coefficients in cm^3.s^{-1}\n\n";

L2:
   cout << " \n Enter iopt\n";
   cin >> iopt;

   int nu = 0;
   int nl = 0;

L3:
   if (iopt == 1)
   {
      cout << " Enter upper and lower n, between " << theory.ncut << " and 2\n";
      cin >> nu >> nl;

      if (nu <= nl)
      {
         cout << " upper n may not be less or equal lower n\n";
         goto L3;
      }

      if (nu > theory.ncut || nl > theory.ncut || nu < 3 || nl < 2 && nl < nu)
      {
         cout << " upper n or lower n outside allowed range (" << theory.ncut << ", 2)" << endl;
         goto L3;
      }
   }
   else if (iopt != 2 && iopt != 3)
   {
      cout << " Invalid value of iopt\n";
      goto L2;
   }

L4:
   cout << " Enter electron number density, temperature and kappa\n";
   cin >> xn >> yt >> zk;

   istop = 0;
   lower = pow(10.0, theory.densx[0]);
   upper = pow(10.0, theory.densx[theory.ndensx - 1]);
   if (xn < lower || xn > upper )
   {
      cout << " Density outside allowable range " << lower << " to " 	<< upper << endl;
      istop = 1;
   }

   lower = pow(10.0, theory.tempx[0]);
   upper = pow(10.0, theory.tempx[theory.ntempx - 1]);
   if (yt < lower || yt > upper )
   {
      cout << " Temperature outside allowable range " << lower << " to " << upper << endl;
      istop = 1;
   }

   lower = pow(10.0, theory.kappax[0]);
   upper = pow(10.0, theory.kappax[theory.nkappax - 1]);
   if (zk < lower || zk > upper)
   {
      cout << " kappa outside allowable range " << lower << " to " << upper << endl;
      istop = 1;
   }

   if (istop == 1)
      goto L8;

   il = (theory.ncut - nu) * (theory.ncut + nu - 1) / 2 + nl;


   for ( k = 0; k < maxx; ++k)
      ri[k] = threedntk(il, xn, yt, zk, ni[k], iopt);


   if (iopt == 1)
   {
      cout  << "\n\n" << setw(53+emis_coef.length()) << emis_coef << endl;
      out16 << "\n\n" << setw(53+emis_coef.length()) << emis_coef << endl;

      cout << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << 2*ni[i] - 1;
      cout << endl;

      out16 << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << 2*ni[i] - 1;
      out16 << endl;

      cout << setw(4) << nu
           << setw(3) << nl
           << setw(10+dens.length()-7) << setprecision(3) << scientific << uppercase << xn
           << setw(5+temp.length())    << setprecision(3) << scientific << uppercase << yt
           << setw(1+kppa.length())    << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      cout << " \n" << endl;

      out16 << setw(4) << nu
            << setw(3) << nl
            << setw(10+dens.length()-7) << setprecision(3) << scientific << uppercase << xn
            << setw(5+temp.length())    << setprecision(3) << scientific << uppercase << yt
            << setw(1+kppa.length())    << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      out16 << "\n" << endl;
   }
   else if (iopt == 2)
   {
      cout  << "\n\n" << setw(53+re_coef_t.length()) << re_coef_t << endl;
      out16 << "\n\n" << setw(53+re_coef_t.length()) << re_coef_t << endl;

      cout << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << 2*ni[i] - 1;
      cout << endl;

      out16 << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << 2*ni[i] - 1;
      out16 << endl;


      cout << setw(10+dens.length()) << setprecision(3) << scientific << uppercase << xn
           << setw(5+temp.length())  << setprecision(3) << scientific << uppercase << yt
           << setw(1+kppa.length())  << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      cout << " \n" << endl;

      out16 << setw(10+dens.length()) << setprecision(3) << scientific << uppercase << xn
            << setw(5+temp.length())  << setprecision(3) << scientific << uppercase << yt
            << setw(1+kppa.length())  << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      out16 << "\n" << endl;
   }
   else if (iopt == 3)
   {
      cout  << "\n\n" << setw(53+re_coef_s.length()) << re_coef_s << endl;
      out16 << "\n\n" << setw(53+re_coef_s.length()) << re_coef_s << endl;

      cout << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << 2*ni[i] - 1;
      cout << endl;

      out16 << setw(10+dens.length()) << dens << setw(5+temp.length()) << temp << setw(1+kppa.length()) << kppa << setw(4+inor.length()) << inor;
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << 2*ni[i] - 1;
      out16 << endl;



      cout << setw(10+dens.length()) << setprecision(3) << scientific << uppercase << xn
           << setw(5+temp.length())  << setprecision(3) << scientific << uppercase << yt
           << setw(1+kppa.length())  << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         cout << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      cout << " \n" << endl;

      out16 << setw(10+dens.length()) << setprecision(3) << scientific << uppercase << xn
            << setw(5+temp.length())  << setprecision(3) << scientific << uppercase << yt
            << setw(1+kppa.length())  << setprecision(3) << scientific << uppercase << zk << setw(4+inor.length()) << " ";
      for ( i = 0; i < maxx; ++i)
         out16 << setw(11) << setprecision(3) << scientific << uppercase << ri[i];
      out16 << "\n" << endl;
   }

L8:
   if (iopt == 1)
      cout <<	 "exit(e), new iopt(i), new transition(n), new physical conditions(c)" << endl;
   else if (iopt == 2 || iopt == 3)
      cout << " exit(e), new iopt(i), new physical conditions(c)" << endl;

   string key;
   cin >> key;

   if (key[0] == opts[0])
   {
      out16.close();
      return;
   }

   if (key[0] == opts[1])
      goto L2;

   if (key[0] == opts[2])
      goto L3;

   if (key[0] == opts[3])
      goto L4;

   cout << " not recognised" << endl;
   goto L8;
}

void hdatk(void)
{
   cout << " reading dataset e1bk.d\n";

   ifstream in15( "e1bk.d" );
   ifstream in16( "t1bk.d" );

   // read e1bk.d dataset
   string line;

   // read the first line
   getline( in15, line );
   istringstream ist( line.c_str() );
   ist >> theory.ndensx >> theory.ntempx >> theory.nkappax;
   ist.clear();

   int ik,ia,ib;

   for ( ik = 0; ik < theory.nkappax; ++ik)
   {
      for (ia = 0; ia < theory.ntempx; ++ia)
      {
         for ( ib = 0; ib < theory.ndensx; ++ib)
         {
            getline( in15, line );
            ist.str(line);
            string str1, str2;

            ist >> str1 >> theory.densx[ib] >> theory.tempx[ia] >> theory.kappax[ik] >> str2 >> theory.nmin >> theory.ncut;

            int ne = theory.ncut * (theory.ncut - 1) / 2;
            ist.clear();

            int j=0;
            int numline = ne/8;
            for (int i=0; i < numline; ++i)
            {
               getline( in15, line );
               ist.str(line);
               for (int k=0; k < 8; ++k)
               {
                  ist >> theory.ex[j][ib][ia][ik];
                  ++j;
               }
               ist.clear();
            }

            // read the remaining data , if any
            if (ne % 8 > 0)
            {
               getline( in15, line );
               ist.str(line);
               double a;
               while (ist >> a)
               {
                  theory.ex[j][ib][ia][ik] = a;
                  ++j;
                  if (j == ne)
                     break;
               }
               ist.clear();
            }
         }
      }
   }

   cout << " read completed" << endl;

   int i,j,k;

   // generate kappa scale factors
   for ( ik = 0; ik < theory.nkappax; ++ik)
   {
      double zk = pow(10.0, theory.kappax[ik]);
      theory.skap[ik] = 1.5*log(zk - 1.5) + gammln(zk - 0.5) - gammln(zk + 1.0);
      theory.skap[ik] = exp(theory.skap[ik]) / (1.0 - 1.5 / zk);
   }

   // read t1bk.d dataset
   cout << " reading dataset t1bk.d\n";

   for ( j = 0; j < theory.nkappax; ++j)
      for ( k = 0; k < theory.ntempx; ++k)
         for ( i = 0; i < theory.ndensx; ++i)
            in16 >> theory.a[i][k][j];

   for ( j = 0; j < theory.nkappax; ++j)
      for ( k = 0; k < theory.ntempx; ++k)
         for ( i = 0; i < theory.ndensx; ++i)
            in16 >> theory.a2s[i][k][j];

   cout << " read completed\n" << endl;

   in15.close();
   in16.close();

   cout << " Minimum principal quantum number = 2\n";
   cout << " Maximum principal quantum number = " << theory.ncut << "\n\n";

   cout << " Range of densities:    "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.densx[0])
        << " - "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.densx[theory.ndensx - 1])   << endl;

   cout << " Range of temperatures: "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.tempx[0])
        << " - "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.tempx[theory.ntempx - 1])   << endl;

   cout << " Range of kappa:        "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.kappax[0])
        << " - "
        << setprecision(3) << scientific << uppercase << pow(10.0, theory.kappax[theory.nkappax - 1]) << endl;
   cout << " \n\n";
}

int main(int argc, char *argv[])
{
   cout << " \n\n           ****** Welcome to Intrat3D *****\n";
   cout << " This code is provided as an accessory to the data associated" << endl;
   cout << " with the following paper:" << endl;
   cout << " P.J. Storey, Taha Sochi. Emission and recombination" << endl;
   cout << " coefficients for hydrogen with kappa electron energy" << endl;
   cout << " distribution. MNRAS XXX." << endl << endl;
   cout << " Error reports, suggestions and enquiries should be sent to" << endl;
   cout << " P.J. Storey (UCL, pjs@star.ucl.ac.uk) or" << endl;
   cout << " Taha Sochi  (UCL, t.sochi@ucl.ac.uk)." << endl << endl;
   cout << " Interpolation of emission and recombination coefficinets" << endl;
   cout << " for hydrogen calculated in Case B with kappa-distributed" << endl;
   cout << " electron energies. Maxwell-Boltzmann results can be" << endl;
   cout << " regained by specifying a large value of kappa (e.g. 10^6)." << endl << endl;
   cout << " Results are given for first, third and fifth order" << endl;
   cout << " Lagrange interpolation in the grid of computed" << endl;
   cout << " coefficients so that convergence can be monitored." << endl << endl;
   cout << " Results are a function of electron number density," << endl;
   cout << " electron temperature and kappa" << endl << endl;

   hdatk();
   getinput();
   return 0;
}