/*
 *	program to simulate waves on a guitar string
 *
 *	use a plucked profile for the initial excitation
 *	
 *	treat a perfectly flexible string
 * 
 *	N. Giordano	2-98
 *
 *	must be compiled with the math library
 *	cc -o guitar guitar.c -lm
 */

#include <stdio.h>
#include <math.h>

/*	these arrays contain the current, old, and new (latest) 
	string profiles		*/
double y_new[2000],y_old[2000],y_current[2000];
double t;	/* time	*/
double x_save,x_pluck;	/* location to pluck the string 	*/
int n_max,n_save;
double dx,dt,c,r2,length;	/* variables discussed in CIP article	*/
FILE *f_save,*f_save_2;

main()
{
	int i;
	int i_loop = 30;	/* loop 30 times through the system before 
					saving the results	*/
	int i_end = 4096;	/* save this many values - convenient for
					FFT analysis		*/

	init();
	while(1) {
		for(i = 0; i < i_loop; i++) propagate();
		save_y();
		if(--i_end == 0) break;
	}
	fclose(f_save);
	fclose(f_save_2);
}

init()
{
	int i,n_pluck;

				/* parameters appropriate for G string (!)
					of a guitar		*/
	length = 0.65;		/* string is 65 cm long	*/
	n_max = 1000;		/* number of spatial steps	*/
	c = 200.0;		/* wave speed for T=120 N and mu=3e-3 kg/m */
	x_pluck = 0.3;
	x_save = 0.10;
	f_save = fopen("y_t","w");  /* save the position of the string here */
	f_save_2 = fopen("f_t","w");/* save the force on the bridge here    */

	dx = length / n_max;
	dt = dx / c;
	n_pluck = x_pluck * n_max;
	n_save = x_save * n_max;
	r2 = dt * dt * c * c / (dx * dx);
	t = 0.0;

				/* some initialization	*/
				/* first the ends	*/
	y_current[0] = y_old[0] = y_current[n_max] = y_old[n_max] = 0.0;
				/* now the initial plucked profile	*/
	for(i = 1; i <= n_pluck; i++) {
			y_current[i] = y_old[i] = ((double)i) / n_pluck;
	}
	for(i = n_pluck+1; i<= n_max-1; i++) { 
		y_current[i] = y_old[i] = 
		-((double)n_max)/(n_pluck-n_max) + ((double)i)/(n_pluck-n_max);
	}

	return;
}

/* make one pass through the system -- update by one time step	*/
propagate()
{
	int i;

	for(i = 1; i <= n_max - 1; i++) 
		y_new[i] = r2*(y_current[i+1]+y_current[i-1])
			+2 * (1.0 - r2) * y_current[i] - y_old[i];
	
	for(i = 1; i <= n_max - 1; i++) {
		y_old[i] = y_current[i];
		y_current[i] = y_new[i];
	}
	t += dt;

	return;
}

/* save the results as discussed above	*/
save_y()
{
	fprintf(f_save,"%g\t%g\n",t,y_current[n_save]);
	fflush(f_save);
	
	fprintf(f_save_2,"%g\t%g\n",t,(y_current[1]-y_current[0])/dx);
	fflush(f_save_2);
	
	return;
}