#ifndef __pest3__
#define __pest3__

/* Pest3 class. A. Pletzer Feb-23-2000
   From various equilibrium file formats compute the
   Delta' stability */

#ifndef __CRAY
#include <string>
using std::string;
#else
#include <mstring.h>
#endif

#include <iostream>
using std::cerr;
#include <stdlib.h>
#include <math.h>
#include "EZcdf.hh"

const int pest3_kmax = 10; // max no of convergence runs
const int pest3_msinmax = 10; // max no rational surfaces
const int pest3_nosingmax = 10; // max no of singular surfaces

/* PEST3 access routines */

extern "C" {

#include "pest3_api.h"

}

const string pest3_inputError     = "--Pest3 Input Error--";
const string pest3_execError      = "--Pest3 Execution Error--";
const string pest3_unknownError   = "--Pest3 Unknown Error--";

class pest3Error{
public:
  pest3Error(string mes){
     cout << mes.c_str() << endl;
  }
};

class pest3 {
public:
  int nErrors;
  string execName;
  int ms; /* resonant poloidal mode */
  int n ; /* toroidal mode */
  int mm[pest3_kmax]; /* number of radial nodes */
  int msin[pest3_msinmax]; /* singular surface selector */
  int msin_len; /* actual size of msin */
  int lmax1; /* Fourier modes -lmax1, ..., lmax1 */
  double b; /* wall distance from magnetic axis is  a*(1+b) */
  double newQ; /* new safety factor value at edge, or if negative on axis */
  int mesh; /* set this to zero to get packed mesh */
  int inputFormat; /*-1=input from binary file INP1 generated by CHEASE
		      0=input TRANSP UFILE
		      1=input from netCDF file inp1.cdf generated by CHEASE
		      2=input from netCDF file eqdsk.cdf generated by JSOLVER 
		      3=input from G-EQDSK file 
		      4=input from G-EQDSK file thru ESC-Q
		      5= Menard's psipqgRZ format in MKS (netCDF)
		      6= Belova's freeqbe format (netCDF)
	      ... other 
		*/
  
  int np1; // no of poloial points + 1 (directed equilibrium input only)
  int nr; // no of radial points (directed equilibrium input only)
  int akima; // =1 if Jacobian should be Akima splined.
  double separatrix; // How close to the separatrix (directed equilibrium input only)?
  double x2axis; // By how much should some profles be extrapolated to the axis.
  double x2edge; // By how much should some profles be extrapolated to the edge.


  int nosurf, mth;

  int nosing;

  double majorRadius, minorRadius, TroyonG, betaPoloidal, betaToroidal;
  double fourPiInductance, totalToroidalCurrent, elongation, triangularity;
  double rMagnetic, b0SquareCentre;
  double mu[pest3_nosingmax];
  double dr[pest3_nosingmax];
  double qs[pest3_nosingmax];
  double q0, qmin, qa, qmax;
  double psi0, psia;
  double psis[pest3_nosingmax];
  double re_deltap[pest3_nosingmax][pest3_nosingmax];
  double im_deltap[pest3_nosingmax][pest3_nosingmax];
  double err_deltap[pest3_nosingmax][pest3_nosingmax];
  double err_gamtap[pest3_nosingmax][pest3_nosingmax];
  double re_gamtap[pest3_nosingmax][pest3_nosingmax];
  double im_gamtap[pest3_nosingmax][pest3_nosingmax];
  int nIdealInstabilities;
  int solver;
  double wLambda;
  double dlayb;
  int lsymhi;
  double betmsh;
  string inputPath; // Can be MDSPLus path,...
  double equilTime;
  int inputPath_len;
  // factor to convert Delta' from psi to minor radius r
  double rdpdr[pest3_nosingmax]; 
  // next order term in large Frobenius expansion in zero beta limit
  double lambdas[pest3_nosingmax]; 
  
  /* c'tor */
  pest3(){

  int i;

  nErrors=0;
  ms = 3;
  n  = 2;
  for(i=1; i<pest3_kmax; i++) mm[i] = 0;
  mm[0]  = 100;

  for(i=1; i<pest3_msinmax; i++) msin[i] = 0;
  msin[0] = ms;

  msin_len=1;
  lmax1 = 10;
  b = 0.0;
  newQ = 0.0;
  mesh = 0;
  inputFormat = +1;

  akima = 0;
  np1 = 129;
  nr = 201;
  separatrix = 0.9999;
  x2axis = 0.01;
  x2edge = 0.00;

  nosurf = 0;
  mth = 0;
  q0 = 0.0;
  qmin = 0.0;
  qa = 0.0;
  qmax = 0.0;
  psi0 = 0.0;
  psia = 0.0;
  nIdealInstabilities = 0;
  majorRadius = 0.0;
  minorRadius = 0.0;
  TroyonG = 0.0;
  betaPoloidal = 0.0;
  betaToroidal = 0.0;;
  fourPiInductance = 0.0;
  totalToroidalCurrent = 0.0;
  elongation = 0.0;
  triangularity = 0.0;
  rMagnetic = 0.0;
  b0SquareCentre = 0.0;
//    mjac = 2;
//    njac = 0;
  solver =0;
  wLambda = 0.0;
  dlayb = 0.9; 
  lsymhi = 0;
  betmsh = 1.0;
  inputPath = "inp1.cdf";
  inputPath_len = (int)inputPath.length();
  equilTime = 0.0;

  //mu = 0.0;
  //dr = 0.0;
  //rdpdr = 0.0; 
  //lambdas = 0.0; 
  //psis = 0.0;
//   re_deltap = 0.0;
//   im_deltap = 0.0;
//   err_deltap = 0.0;
//   re_gamtap = 0.0;
//   im_gamtap = 0.0;
//   err_gamtap = 0.0;
  }

  /* d'tor */
  ~pest3(){
  }

  void usage(void);      
  int parseArgs(int narg, char *argc[]);
  int checkArgs(void);
  int isQsInPlasma(void);
  int run(void);
  
};


  void pest3::usage(void){
    cout << "Usage: pest3 [Options]\n";
    cout << "Options:\n";
    cout << "-n<ns> Set the toroidal mode number (2)\n";
    cout << "-m<ms> Set the resonant poloidal mode (3). Alternatively, \n";
    cout << "\t the user can also give a string of the form \"..x.x\" to \n";
    cout << "\t specify the rational surfaces where the Delta' are to be \n";
    cout << "\t computed: an 'x' for a selected surface and a '.' for surfaces\n";
    cout << "\t that should be skipped. Each symbol ('.' or 'x') corresponds\n";
    cout << "\t to a surface where ns*q is an integer, counting from the axis.\n";
    cout << "-k<Ns> Set number of radial finite elements (100)\n";
    cout << "\t To perform a convergence study type -k\"Ns1 Ns2 ...\" where\n";
    cout << "\t Ns1, Ns2, .. are a range of radial number of finite elements.\n";
    cout << "\t The matching data so obtained are the extrapolation to inifinte\n";
    cout << "\t resolution (Ns->oo) and an estimate of the finite element error\n";
    cout << "\t is returned assuming a convergence in 1/Ns^2. Note that that the \n";
    cout << "\t Frobebius coefficients of the large solution are only computed\n";
    cout << "\t once, for Ns1. A typical choice is -k\"100 70 140 200\".\n";
    cout << "-l<Lmax1> Poloidal Fourier modes span -|Lmax1|,...+|Lmax1| (10)\n";
    cout << "-b<Bw> Relative wall distance from plasma surface normalized to \n";
    cout << "\t minor radius (0). For instance -b0.2 for wall position at 1.2\n";
    cout << "\t times minor radius. Bw=0 is wall on plasma surface whereas Bw>=10 is\n";
    cout << "\t wall at infinity.\n";
    cout << "-S for symmetric large solution support. Default is asymmetric.\n";
    cout << "\t This option is typically used for higher m modes, which tend to \n";
    cout << "\t be localized about their resonant surface\n";
    cout << "-d<rw> Relative extension of large solution to neighboring rational\n";
    cout << "\t surfaces (0.9). Must have 0.05 < rw < 0.99. If -S is in addition\n";
    cout << "\t selected, then rw is the smallest of the distance between the left\n";
    cout << "\t and right rational surfaces. \n";
//      cout << "-e<eb> For packing near plasma edge set eb > 1.0 (1.0).\n";
    cout << "-i<if> Format (1). Possible values are:\n";
    cout << "\t-1 read binary equilibrium input file (INP1) from CHEASE\n";
    cout << "\t 0 input equilibrium from TRANSP Ufiles (local or MDS+)\n";
    cout << "\t 1 read netCDF equilibrium input file (inp1.cdf) from CHEASE\n";
    cout << "\t 2 read netCDF equilibrium input file (eqdsk.cdf) from JSOLVER\n";
    cout << "\t 3 read G-EQDSK equilibrium input file from EFIT. This can also\n";
    cout << "\t   be an MDSPlus path, for instance:\n";
    cout << "\t   MDS+/REDUCE:BIRCH.PPPL.GOV:8501:EFIT01(103984;t=0.23)\n";
    cout << "\t 4 read G-EQDSK equilibrium input file from EFIT and rerun thru\n";
    cout << "\t   ESC-Q. Can also be an MDSPlus path (see above).\n";
    cout << "\t 5 read Menard netCDF input file in psipqgRZ MKS format\n";
    cout << "\t 6 read Belova netCDF input file format (freeqbe)\n";
    cout << "-f<File> Read equilibrium data from <File>. <File> can be a local \n";
    cout << "\t netCDF file (e.g. inp1.cdf) or, if MDSPlus access has been enabled\n";
    cout << "\t a path such as: MDS+:TRANSPGRID.PPPL.GOV:TRANSP(TFTR.88,37065Y82)\n";
    cout << "-T<time> Set MDSPlus equilibrium time to <time> in [sec] (0). \n";
    cout << "\t Will correct to the min/max value if selected <time> is outside \n";
    cout << "\t interval.\n";
    cout << "-q<newQ> Rescale equilibrium by adjusting the safety factor q=newQ\n";
    cout << "\t at the plasma edge. If negative, then set q on axis to abs(newQ).\n";
    cout << "\t Leave equilibrium unchanged if newQ=0 (default).\n";
    cout << "-s<is> Solver (0). Possible values are:\n";
    cout << "\t 0 use PEST3's built-in linear system solver\n";
    cout << "\t 1 compute lowest delta-W eigenvalue using LAPACK routine\n";
    cout << "-u Uniform radial mesh. Default is mesh packing about rational surfaces.\n";
    cout << "\t In general, automatic mesh packing is desirable. This option is useful\n";
    cout << "\t at high finite element resolution by preventing the clutering of nodes\n";
    cout << "\t about the resonant surface.\n";

    cout << "-a<xa> Extrapolate some profiles over distance xa to axis (0.01)\n";
    cout << "-x<xe> Extrapolate some profiles over distance xe to edge (0.00)\n";
    
    cout << "-P<p1> Use p1 poloidal rays for mapping (129)(direct equilibrium only).\n";
    cout << "-R<rs>  Use rs radial surfaces for mapping (401)(direct equilibrium only).\n";
    cout << "-E<ed>  Separatrix retreat factor (0.999) (direct equilibrium only).\n";

    cout << "-A To Akima spline Jacobian and friends (default is Jacobian derived \n";
    cout << "\t from X and Z).\n";
    
    
    cout << "-h This message.\n"; 
  }        

  int pest3::parseArgs(int narg, char *argc[]){
   int j, flag_loc, k, i;

    execName = string(argc[0]);
    for (i=1; i<narg; i++){
      switch(argc[i][1]){
    case 'm':
       if(isdigit(argc[i][2])) {
	// if a number then must be resonant m
	ms = atoi(&argc[i][2]);
	msin[0] = ms;
	msin_len = 1;
	// msin set later
       } else {
 	// user knows where the resonant surfaces are, and
 	// uses a string format: eg ".**." 
 	ms = 0;
 	j = 2;
 	msin_len = 0;
 	while (argc[i][j] != '-' && argc[i][j] != '\0' && 
 	       msin_len < pest3_msinmax){
	  if(argc[i][j] == '"') continue;
 	  msin[msin_len] = 1;
 	  if(argc[i][j] == '.') {
	    msin[msin_len] = 0;
	    cout << " Rational surface no. " << j-1 << " will be skipped\n";
	  } else {
	    cout << " Rational surface no. " << j-1 << " will be added\n";
	  }
 	  msin_len++;
 	  j++;
 	}
      }
      break;
    case 'n':
      n = atoi(&argc[i][2]);
      break;
    case 'k':
     mm[0] = atoi(&argc[i][2]);
     j = 2;
     k = 0;
     flag_loc = 1;
     while (argc[i][j] != '-' && argc[i][j] != '\0'){
	if(flag_loc==0 && argc[i][j]==' ') flag_loc=1;
	if(flag_loc==1 && argc[i][j]!=' ' && k<pest3_kmax){
	  mm[k] = atoi(&argc[i][j]);
	  flag_loc = 0;
	  k++;
	}
	j++;
     }  
      break;
    case 'l':
      lmax1 = atoi(&argc[i][2]);
      break;
    case 'b':
      b = atof(&argc[i][2]);
      break;
    case 'e':
      betmsh = atof(&argc[i][2]);
      break;
    case 'd':
      dlayb = atof(&argc[i][2]);
      break;
    case 'S':
      lsymhi = 1;
      break;
    case 'i':
      inputFormat = atoi(&argc[i][2]);
      break;
    case 'q':
      newQ = atof(&argc[i][2]);
      break;
    case 'u':
      mesh = 1;
      break;
    case 's':
      solver = atoi(&argc[i][2]);
      break;
    case 'f':
      inputPath = string((const char*) (&argc[i][2]));
      break;
    case 'T':
      equilTime = (double) atof(&argc[i][2]);
      break;
    case 'a':
      x2axis = (double) atof(&argc[i][2]);
      break;
    case 'x':
      x2edge = (double) atof(&argc[i][2]);
      break;
    case 'A':
      akima = 1;
      break;
    case 'P':
      np1 = atoi(&argc[i][2]);
      break;
    case 'R':
      nr = atoi(&argc[i][2]);
      break;
    case 'E':
      separatrix = (double) atof(&argc[i][2]);
      break;
    case 'h':
      usage();
      exit(1);
      return 0;
    default:
      cerr << pest3_inputError << endl;
      usage();
      return (++nErrors);
    }
  }
      return nErrors;
  }

int pest3::checkArgs(void){
  int mm_min = mm[0];
  int i;
  for(i=1; i<pest3_kmax; i++){
    if(mm[i] < mm_min && mm[i]!=0) mm_min = mm[i];
  }

  if(ms < 0){
    cerr << pest3_inputError << " resonant poloidal mode number " << ms << " < 1 \n";
    return (++nErrors);
  } else if(n < 1){
    cerr << pest3_inputError << " toroidal mode number " << n << " < 1 \n";
    return (++nErrors);
  } else if(mm_min < 5){
    cerr << pest3_inputError << " min number of radial nodes " << mm_min << " < 5 \n";
    return (++nErrors);
  } else if(lmax1 < 1){
    cerr << pest3_inputError << " max poloidal Fourier mode " << lmax1 << " < 1 \n";
    return (++nErrors);
  } else if(b < 0.0){
    cerr << pest3_inputError << " negative wall distance " << b << endl;
    return (++nErrors);
  } else if(betmsh < 0.0){
    cerr << pest3_inputError << " negative edge packing exponent " << betmsh << endl;
    return (++nErrors);
  } else if(dlayb < 0.05 || dlayb > 0.99){
    cerr << pest3_inputError << " large solution support " << dlayb << 
      " outside [0.05, 0.99]\n";
    return (++nErrors);
  } else if(inputFormat < -1 || inputFormat > 6){
    cerr << pest3_inputError << " Wrong equilibrium input format " << inputFormat << endl;
    return (++nErrors);
  } else if(x2axis > 0.5 || x2axis < 0.){
    cerr << pest3_inputError << " Extrapolation to axis fraction too large or negative " 
	 << x2axis << " > 0.5\n";
  } else if(x2edge > 0.5 || x2edge < 0.){
    cerr << pest3_inputError << " Extrapolation to edge fraction too large or negative " 
	 << x2edge << " > 0.5\n";
  } else if(np1 < 17){
    cerr << pest3_inputError << " Number of poloidal contouring points too small " 
	 << np1 << " < 17\n";
  } else if(nr < 11){
    cerr << pest3_inputError << " Number of radial contouring points too small " 
	 << nr << " < 11\n";
  } else if(separatrix > 1. || separatrix < 0.1){
    cerr << pest3_inputError << " Separatrix retreat factor should be close to <=1 " 
	 << separatrix << endl;
  } 
  else{
    return nErrors;
  }
}

int pest3::isQsInPlasma(void){

  pstgetqbounds_(&q0, &qa, &qmin, &qmax);
  cout << " Equilibrium q axis/edge-> " <<  q0   << '/' << qa << endl;
  cout << " Equilibrium q min /max -> " << qmin << '/' << qmax << endl;

  if(ms != 0){
    double qs =  (double) ms/ (double) n;
    if( qs < qmin || qs > qmax ){
      cerr << pest3_execError << " No " << ms << "/" << n << " resonance \n";
      ++nErrors;
    }
  } else {
    if( qmax < qmin + 1.0/double(n) ){
      cerr << pest3_execError << " No resonant surface \n";
      ++nErrors; 
    }
  }
    
  return nErrors;
}


int pest3::run(void){

    
  /* Initialize pest3. This step involves reading the equilibrium input
     file in. */

    int i;

    int nmeshMax = mm[0];
    for(i=1; i<pest3_kmax; i++){
      if(mm[i] > nmeshMax) nmeshMax=mm[i];
    }

    int ier;
    long long_len = (long) inputPath_len; // pass string length as last arg
    inputPath_len = (int) inputPath.length();
    pstinit_(&lmax1, &nmeshMax, &inputFormat, 
	     (char *)inputPath.c_str(), &inputPath_len, &equilTime, 
	     &newQ, &akima, &x2axis, &x2edge, &np1, &nr, &separatrix, 
	     &ier,
	     long_len);
    if(ier != 0){
      cerr << pest3_execError << " Failed to initialize PEST3\n";
      cerr << " There can be several reasons for this to occur\n";
      cerr << " 1. Check that the input file/MDSPlus path exists and is valid\n";
      cerr << " 2. Check that the plasma is NOT Mercier unstable\n";
      return ++nErrors;
    }

    double gserror;
    pstgetgsaverageerror_(&gserror);
    cout << " Average Grad Shafranov solution error->" << gserror << endl;
    if(fabs(gserror) > 0.001){
      cout << " *WARNING* inaccurate equilibrium\n";
    }
//      if(fabs(gserror) > 0.01){
//        cout << " *ERROR* inaccurate MHD equilibrium solution \n";
//        return (++nErrors);
//      }

    /* set input pest3 parameters */

    //int one = 1;
    pstsetn_(&n);
    if( isQsInPlasma() != 0 ){
      return (++nErrors);
    }
    
    /* set the msin variable in pest3, which selects the resonant 
       rational surfaces counting from the axis towards the edge.
       For instance, msin = [0, 1, 0, ..] selects the second rational
       surface. For simplicity, the q profile is assumed monotonic. */

    pstsetmsin_(&msin_len, &msin[0]); 

    /* mm, the number of radial nodes, can in general be an array,
       which is useful for performing convergence studies.  */

    int kdim = (int) pest3_kmax;
    pstsetmm_(&kdim, &mm[0]);
    
    pstsetb_(&b);
    pstsetdlayb_(&dlayb);
    pstsetlsymhi_(&lsymhi);
    pstsetbetmsh_(&betmsh);
    pstsetisolver_(&solver);
    pstsetuniformmesh_(&mesh);

    /* run Pest3 */

    pstexec_();

    pstgetmth_(&mth);
    pstgetnosurf_(&nosurf);
    pstgetbetapoloidal_(&betaPoloidal);
    pstgetbetatoroidal_(&betaToroidal);
    pstgettroyong_(&TroyonG);
    pstgetfourpiinductance_(&fourPiInductance);
    pstgettotaltoroidalcurrent_(&totalToroidalCurrent);
    pstgetmajorradius_(&majorRadius);
    pstgetminorradius_(&minorRadius);
    pstgetelongation_(&elongation);
    pstgettriangularity_(&triangularity);
    pstgetb0squarecentre_(&b0SquareCentre);
    pstgetrmagnetic_(&rMagnetic);
    pstgetwlambda_(&wLambda);

    /* layer data */

    pstgetnosing_(&nosing);

    double *twoMuF = new double[pest3_nosingmax];

    pstgetcmatch_(&nosing, &twoMuF[0]);
    pstgetpsisin_(&nosing, &psi0, &psis[0], &psia);
    pstgetxmu_(&nosing, &mu[0]);
    pstgetxsmnus_(&nosing, &lambdas[0]);
    pstgetrsnorm_(&nosing, &rdpdr[0]);
    pstgetdrlay_(&nosing, &dr[0]);
    pstgetqslay_(&nosing, &qs[0]);
    

    /* Delta'. Gamma' etc.. */

    pstgetdelpr_((int *) &pest3_nosingmax, &nosing, 
		 &re_deltap[0][0], &im_deltap[0][0]);
    pstgetgampr_((int *) &pest3_nosingmax, &nosing, 
		 &re_gamtap[0][0], &im_gamtap[0][0]);
    pstgeterror_delpr_((int *) &pest3_nosingmax, &nosing, &err_deltap[0][0]);
    pstgeterror_gampr_((int *) &pest3_nosingmax, &nosing, &err_gamtap[0][0]);


    /* get the log terms */

    double *re_xiLog = new double[pest3_nosingmax];
    double *im_xiLog = new double[pest3_nosingmax];
    pstgetxilogterm_(&nosing, &re_xiLog[0], &im_xiLog[0]);
      

    /* Normalize to psi^(2 mu) etc. */

    for(i = 0; i < nosing; i++){

      if(mu[i] != 0.5) {
	lambdas[i] *= -(1.0-2.*mu[i]); 
      } else {
	lambdas[i] = re_xiLog[i];
      }
      lambdas[i] *= twoMuF[i]*psis[i];

      rdpdr[i] /= (twoMuF[i]*pow(psis[i], 2.0*mu[i]));
      for(int j = 0; j < nosing; j++){
	double fact = sqrt(twoMuF[i]*twoMuF[j]) 
	  * pow(psis[i], mu[i]) 
	  * pow(psis[j], mu[j]);
	re_deltap[i][j] *= fact;
	im_deltap[i][j] *= fact;
	err_deltap[i][j] *= fact;
	re_gamtap[i][j] *= fact;
	im_gamtap[i][j] *= fact;
	err_gamtap[i][j] *= fact;
      }
    }

    delete[] twoMuF;
    delete[] re_xiLog;
    delete[] im_xiLog;
    
    // nodes where ideal instability, if any, was detected
    int dum;
    int zero = 0;
    pstgetineg_(&nIdealInstabilities, &dum, &zero);
    if (wLambda != 0.0){
      cout << " delta-W eigenvalue " << wLambda << endl;
    }
    if (nIdealInstabilities != 0){
      // get the nodes where the instability occurs
      int nsize=nIdealInstabilities;
      int *nodes = new int[nsize];
      pstgetineg_(&nIdealInstabilities, &nodes[0], &nsize);
      cout << " *WARNING*: " << nsize << " ideal instabilities detected at nodes ";
      for(i=0; i<nsize; i++){
	cout << nodes[i] << ' ';
      }
      cout << '\n';
      delete[] nodes;
    }


    /* free memory */

    pstfree_();

    return nErrors;
}
    
#endif /* __pest3__ */
