/* 
************************************************************************
*  rk4.c: 4th order Runge-Kutta solution for harmonic oscillator       *
*                      *
* From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
   by RH Landau, MJ Paez, and CC BORDEIANU 
   Copyright Princeton University Press, Princeton, 2008.
   Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
   MJ Paez, Univ Antioquia, 2008; & CC BORDEIANU, Univ Bucharest, 2008
   Support by National Science Foundation                              
*
************************************************************************
*/
       
#include <stdio.h>

#define N 2                                   /* number of equations */
#define dist 0.01                              /* stepsize */
#define MIN 0.0                               /* minimum x */
#define MAX 10.0                              /* maximum x */

void runge2(double x, double y[], double step);
double f(double x, double y[], int i);


main()
{
	double x, y[N];
	int j;
   
	FILE *outputp, *outputv;
	outputp = fopen("rk2p.dat","w");
	outputv = fopen("rk2v.dat","w");
   
	y[0] = 1.0;                                /* initial position  */
	y[1] = 0.0;                                /* initial velocity  */
   
	fprintf(outputp, "%f\t%f\n", x, y[0]);
	fprintf(outputv, "%f\t%f\n", x, y[1]);

	for(x = MIN; x <= MAX ; x += dist)
	{
		runge2(x, y, dist);
		fprintf(outputp, "%f\t%f\n", x, y[0]);   /* position vs. time */
		fprintf(outputv, "%f\t%f\n", x, y[1]);   /* velocity vs. time */
	}

	fclose(outputv);
	fclose(outputp);
}
/*-----------------------end of main program--------------------------*/

/* Runge-Kutta subroutine */
void runge2(double x, double y[], double step)
{
	double h=step/2.0,                         /* the midpoint */
		t1[N], t2[N],
		k1[N], k2[N];          /* for Runge-Kutta  */
	int i;
 
	for (i=0; i<N; i++) t1[i] = y[i]+0.5*(k1[i]=step*f(x, y, i));
	for (i=0; i<N; i++) k2[i]=step*f(x+h, t1, i);

	for (i=0; i<N; i++) y[i] += k2[i];
}
/*--------------------------------------------------------------------*/

/* definition of equations - this is the harmonic oscillator */
double  f(double x, double y[], int i)
{
	double k = 4,
	       m = (1/3.14159)*(1/3.14159);
	if (i == 0) return(y[1]);               /* RHS of first equation */
	if (i == 1) return(k/m*(-y[0]));              /* RHS of second equation */

}