/*************************    Program 2.4    ****************************/
/*                                                                      */
/************************************************************************/
/* 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 NMAX 101

main()
/* Main program for derivatives of f(x) = sin(x).  F1: f';
   F2: f";  D1: error in f'; and D2: error in f". 
   Copyright (c) Tao Pang 1997. */
{
int i,n;
double pi,h;
double x[NMAX],f[NMAX],f1[NMAX],d1[NMAX],f2[NMAX],d2[NMAX];
void three();

pi = 4*atan(1);
n = NMAX;
h  = pi/(2*(n-1));
for (i=0; i < n; ++i)
  {
  x[i]  = h*i;
  f[i]  = sin(x[i]);
  f1[i] = 0;
  f2[i] = 0;
  }
three (n,h,f,f1,f2);
for (i=0; i < n; ++i)
  {
  d1[i] = f1[i]-cos(x[i]);
  d2[i] = f2[i]+sin(x[i]);
  printf("%10.6lf %10.6lf %10.6lf %10.6lf %10.6lf\n",
          x[i],f1[i],d1[i],f2[i],d2[i]);
  }
}

void three(n,h,f,f1,f2)
/* Subroutine for 1st and 2nd order derivatives with the three-
   point formulas. Extrapolations are made at the boundaries.
   FI: input f(x); H: interval; F1: f'; and F2: f".
   Copyright (c) Tao Pang 1997. */
int n;
double h;
double f[],f1[],f2[];
{
int i;
/* f' and f" from three-point formulas */

for (i = 1; i < n-1; ++i)
  {
  f1[i] = (f[i+1]-f[i-1])/(2*h);
  f2[i] = (f[i+1]-2*f[i]+f[i-1])/(h*h);
  }
/* Linear extrapolation for the boundary points */

f1[0]   = 2*f1[1]-f1[2];
f1[n-1] = 2*f1[n-2]-f1[n-3];
f2[0]   = 2*f2[1]-f2[2];
f2[n-1] = 2*f2[n-2]-f2[n-3];
}