/* Author:   Clare Quarman
 * Date:     Sat 8th Nov 2003
 * Purpose:  (1) rewriting of quprog.f as C++, therefore...
 *           (2) practice calling fortran functions & routines from C++
 */

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <string>
#include "QCDnum.cpp"

using namespace QCDnum;
using namespace FortranIO;
using std::cout;
using std::endl;
using std::string;

double xgluons(double x){
  double result = 1.01*pow(x,-0.35)*pow((1.0-x),5.22);
  cout << "xgluons \t" << result << endl;
  return result;
}

double xquarks(double x){
  return 1.01*pow(x,-0.35)*pow((1.0-x),5.22);
}

int main(){

  //-1-  initialisation
  cout << "1" <<endl;
  qninit();
  string version(""); // you have to know to initialise it to "".
  char versch[9]; // you have to know it should be 9 long
  string versstr="blar";
  string *versstr2;
  bool ifDP=true, ifDP2;
  int nxmax=0, nqmax=0;
  int nxmax2, nqmax2;
  qnvers1(version,  ifDP, nxmax, nqmax);
  qnvers2(versch,   ifDP, nxmax, nqmax);
  qnvers3(versstr,  ifDP, nxmax, nqmax);
  //qnvers0(versstr2, ifDP2, nxmax2, nqmax2);
  cout << "QCDNUM version is " << version  << endl;
  cout << "QCDNUM version is " << versch   << endl;
  cout << "QCDNUM version is " << versstr  << endl;
  string empty;
  bool ifDPempty;
  int iempty,jempty;
  qnvers4(empty,  ifDPempty, iempty, jempty);
  cout << "QCDNUM version is " << empty    << "\t qnvers4"<< endl;
  //  cout << "QCDNUM version is " << *versstr2 << endl;
  qntime("start");  


  //-2-  Setup grid in x and Q2
  cout << "2" <<endl;
  double q0=7.0;
  grxdef(80, 4.0e-4);
  grqdef(59, 4.0e0, 5000.0e0);
  int iarg3 = 1;
  double* pq0 = &q0;
  grqinp(pq0, 1); /* grqinp expects an array so have to give it a pointer
		     if only this one value - this is contrary to the
		     fortran example (quprog.f) as fortran treats a single
		     number as essentially an array of numbers with lenth 1 anyway.*/

  // test if I've done grxinp wrapper correctly - yep!
  double extraX[] = {0.5, 0.7};
  grxinp(extraX, 2);

  int nx, nq;
  double xmi, xma, qmi, qma;
  grgive(nx, xmi, xma, nq, qmi, qma);

  //-3-  Create weight tables
  cout << "3" <<endl;
  qnfilw(0, 0);

  //-3a- test the writing of weight table to file and reading from file
  cout << "3a" << endl;;
  int lun = 24;
  string file = "weights.txt";
  openFileWithLun(lun, file);
  qndump(lun);
  int istop = 0, ierr = 0; // see top of pg 49 of manual for istop & ierr meaning
  qnread(lun, istop, ierr);
  closeFileWithLun(lun);
  cout << "   ierr = " << ierr << endl;

  //-4-  Set charm, bottom thresholds at Q2 = 4, 25 GeV2
  cout << "4" <<endl;
  double cthres = 4.0e0, bthres = 25.0e0;
  qthres(cthres, bthres);

  //-5-  Set input value for alphas
  cout << "5" <<endl;
  double mz = 91.2e0, a_s = 0.118e0;
  qnrset("alfas", a_s);
  qnrset("alfq0", mz*mz);

  //-6-  Input singlet and gluon at Q2 = 7 Gev2
  //     quarks(x) and gluons(x) are user defined functions
  cout << "6" <<endl; 
  int iq0 = iqfromq(q0);
  for (int ix = 1; ix <= nx; ++ix){
    double x = xfromix(ix);
    qnpset("singlet", ix, iq0, xquarks(x));
    qnpset("gluon"  , ix, iq0, xgluons(x));
  }

  //-7-  Singlet/gluon evolution
  cout << "7" <<endl; 
  cout << "nq:  " << nq << endl;
  cout << "qma: " << qma << endl;
  evolsg(iq0, 1, nq);

  //-8-  Get singlet, gluon, F2 etc. at any x and Q2
  cout << "8" <<endl; 
  double x    = 0.005;
  double q2   = 20.0;
  int iflag;
  double qval = qpdfxq(      "singlet", x, q2, iflag);
  cout << iflag <<endl;
  double gval = qpdfxq(      "gluon"  , x, q2, iflag);
  cout << iflag <<endl;
  double fval = qstfxq("F2", "singlet", x, q2, iflag);
  cout << iflag <<endl;

  cout << "####################" << endl;
  cout << "x    " << x           << endl;
  cout << "q2   " << q2          << endl;
  cout << "qval " << qval        << endl;
  cout << "gval " << gval        << endl;
  cout << "fval " << fval        << endl;
  cout << "####################" << endl;

  for (q2 = q0; q2<20.0; q2 += 1.0){
    qval = qpdfxq(      "singlet", x, q2, iflag);
    gval = qpdfxq(      "gluon"  , x, q2, iflag);
    fval = qstfxq("F2", "singlet", x, q2, iflag);
    cout << "####################" << endl;
    cout << "x    " << x           << endl;
    cout << "q2   " << q2          << endl;
    cout << "qval " << qval        << endl;
    cout << "gval " << gval        << endl;
    cout << "fval " << fval        << endl;
    cout << "####################" << endl; 

  }

  //-9-  Some printout
  cout << "9" <<endl; 
  int output = 6;
  qprint(output,"all");


  //-10- Finished with QCDnum
  cout << "10" <<endl;
  endQCDnum();
   
  return EXIT_SUCCESS;
}