/* This is a driver program for root-finding using the 
    Newton-Raphson algorithm as implmented in the Numerical
    Recipes library.  Although general for any function, 
    we use it here to find a root of the Van der Waals 
    equation of state discussed in class.

    G.C. Rutledge
    October, 2002
*/

#include "nr.h"                 /* prototypes zbrac and rtsafe */
#include "nrutil.h"             /* prototypes nrerror */
#include <stdio.h>
#define GAS_CONSTANT 0.08206    /* liter*atm/mol/K */
#define A 5.464                 /* liter^2*atm/mole^2, for H2O */
#define B 0.03049               /* liter/mole, for H2O */
#define SWAP(a,b) {float temp;temp=(a);(a)=(b);(b)=temp;}     
                                /* macro to swap float values */
float func_vdw(float);          /* prototype for func_vdw() */
float deriv_func_vdw(float);    /* prototype for deriv_func_vdw() */
void func(float,float *,float *); /* prototype for function call */

static float pressure,temperature;  /* these are needed in function calls*/

int main(void)
{
  float xlo,xhi,xacc;

  /*Communication with user...    */
  printf("input Pressure and Temperature, in atm and K:\n");
  scanf("%f%f",&pressure,&temperature);

  /*From here on, the routine is general */
  printf("Input an initial guess for an interval to bracket the root:\n");
  scanf("%f%f",&xlo,&xhi);
  printf("Input desired accuracy in f(x)\n");
  scanf("%f",&xacc);
  if(xlo>xhi) SWAP(xlo,xhi);

  /* use NR library routine to ensure that bracket is good */
  if(zbrac(func_vdw,&xlo,&xhi)!=1) 
    nrerror("unacceptably large bracket in zbrac");

  /* use NR library routine to locate root */
  printf("The root of the equation is %f\n",rtsafe(func,xlo,xhi,xacc));
  return (0);
}

/* This function calls func and deriv_func to satisfy syntax of rtsafe */
void func(float v,float *f,float *df)
{
  *f  = func_vdw(v);
  *df = deriv_func_vdw(v);
  printf("%f %f\n",*f,*df);
  return;
}

/* This function returns the residual for the Van der Waals
   equation of state */
float func_vdw(float v)
{
  if(v==0.0) nrerror("volume cannot be zero");
  return ( (pressure+A/v/v)*(v-B)-GAS_CONSTANT*temperature );
}
  
/* This function returns the derivative of the residual with
   respect to volume for the Van der Waals equation of state */
float deriv_func_vdw(float v)
{
  float v2;
  if(v==0.0) nrerror("volume cannot be zero");
  v2=v*v;
  return (pressure+A/v2-2.*A*B/v2/v);
}