#include <iostream>
#include "nrutil.h"
using namespace std;
//     Here we define various functions called by the main program

double E_function(double *x);
void dE_function(double *x, double *g);
void dfpmin(double p[], int n, double gtol, int *iter, double *fret,
	    double(*func)(double []), void (*dfunc)(double [], double []));
//   Main function begins here
int main()
{
     int n, iter;
     double gtol, fret;
     double alpha;
     n = 1;
     cout << "Read in guess for alpha" << endl;
     cin >> alpha;
//   reserve space in memory for vectors containing the variational
//   parameters
     double *p = new double [2];
     gtol = 1.0e-5;
//   now call dfmin and compute the minimum
     p[1] = alpha;
     dfpmin(p, n, gtol, &iter, &fret,&E_function,&dE_function);
     cout << "Value of energy minimum = " << fret << endl;
     cout << "Number of iterations = " << iter << endl;
     cout << "Value of alpha at minimu = " << p[1] << endl;
      delete [] p;
      return 0;
}  // end of main program
//  this function defines the Energy function
double E_function(double x[])
{
  double value = x[1]*x[1]*0.5+1.0/(8*x[1]*x[1]);
  return value;
} // end of function to evaluate

//  this function defines the derivative of the energy 
void dE_function(double x[], double g[])
{
  g[1] = x[1]-1.0/(4*x[1]*x[1]*x[1]);
} // end of function to evaluate