/*************************    Program 2.9    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define f(x) (e2/(x*x)-a0*exp(-x/r0)/r0)
#define g(x) (e2/x-a0*exp(-x/r0))
double e2=14.4,a0=1090,r0=0.33;

main()
/* Main program to calculate the bond length of NaCl.
   Copyright (c) Tao Pang 1997. */
{
int istep;
double dl=1e-6;
double xi,di,xf,df;
void m_secant();

xi = 2;
di = 0.1;
m_secant (dl,xi,di,&xf,&df,&istep);
printf("%4d %16.8lf %16.8lf\n", istep,xf,df);
}

void m_secant (dl,xi,di,xf,df,istep)
/* Subroutine for the root of f(x) = dg(x)/dx = 0 with the
   secant method with the search toward the maximum of g(x).
   Copyright (c) Tao Pang 1997. */
int *istep;
double dl,xi,di;
double *xf,*df;
{
double d,g0,g1,g2,x0,x1,x2;

x0 = xi;
g0 = g(x0);
x1 = x0+di;
g1 = g(x1);
*df = di;
*istep = 0;
if (g1 < g0) x1 = x0-*df;
while (fabs(*df) > dl)
  {
  d = f(x1)-f(x0);
  *df = -(x1-x0)*f(x1)/d;
  x2 = x1+*df;
  g2 = g(x2);
  if (g2 < g1) x2 = x1-*df;
  x0 = x1;
  x1 = x2;
  g1 = g2;
  *istep = *istep+1;
  }
  *xf = x2;
}