/*
 *	program to simulate waves on a piano string
 *
 *	the string is perfectly flexible, and is excited by a blow from
 *	a hammer, which is modeled by the force law F=Kz^p
 *
 *	need to simulate both the string and the hammer
 *
 *	the ends of the string are held fixed (i.e., perfectly rigid)
 *
 *	N. Giordano    2-98
 *
 *	to compile:	cc -o hammer hammer.c -lm
 */

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

#define		YES	1
#define		NO	0

			/* old, current, and new string profiles	*/
double y_new[2000],y_old[2000],y_current[2000];
double y_hammer,v_hammer,K,p,m_hammer,mu,dm;  /* notation follows CIP article*/
int hammer_contact;	/* flag which indicates if hammer is in contact with
				the string or not	*/
double t;
double x_save,x_hammer;
int n_max,n_save,n_hammer;
double dx,dt,c,r2;
FILE *f_save;

main()
{
	int i;
	int i_loop = 30;	/* make i_loop updates of the system before
					storing new values	*/
	int i_end = 1024;

	i_loop = 30;

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

init()
{
	int i;
	                   	/* all units are mks	*/
	m_hammer = 3e-3;	/* hammer mass		*/
	mu = 3e-3;		/* mass per unit length of string	*/
	p = 2.5;		/* force exponent	*/
	K = 5e9;		/* force coefficient	*/
	v_hammer = 2.0;		/* velocity and position of hammer	*/
	y_hammer = 0.0;
	hammer_contact = YES;
	
	n_max = 1000;
	c = 300.0;
	x_hammer = 1.0 / 7.0;	/* hammer strike point	*/
	t = 0.0;

	dx = 1.0 / n_max;	/* string is 1 m long	*/
	dt = dx / c;
	dm = mu * dx;
	n_hammer = x_hammer * n_max;
	
	x_save = 0.10;
	n_save = x_save * n_max;
	f_save = fopen("y_t","w");

	for(i = 0; i <= n_max; i++) y_new[i] = y_current[i] = y_old[i] = 0.0;
	
	return;
}

propagate()
{
	int i;
	double force;

	r2 = dt * dt * c * c / (dx * dx);

		/* assume a perfectly flexible string -- no stiffness	*/
	for(i = 1; i <= n_max - 1; i++) 
		y_new[i] = 2 * (1-r2) * y_current[i] - y_old[i]
			+ r2 * (y_current[i+1] + y_current[i-1]);

		/* hammer force is distributed over spatial units n_hammer 
			and n_hammer +- 1	*/
		/* force = K |y_hammer - y_string|^p		*/

	if(hammer_contact == YES) {
		force = K * pow(y_hammer - y_current[n_hammer],p);
		y_new[n_hammer-1] += 0.25 * dt * dt * force / dm;  
		y_new[n_hammer] += 0.5 * dt * dt * force / dm;  
		y_new[n_hammer+1] += 0.25 * dt * dt * force / dm;  
		v_hammer = v_hammer - dt * force / m_hammer;
		y_hammer = y_hammer + v_hammer * dt;
		if(y_hammer < y_new[n_hammer]) hammer_contact = NO;
	}

	for(i = 1; i <= n_max - 1; i++) {
		y_old[i] = y_current[i];
		y_current[i] = y_new[i];
	}

	t += dt;

	return;
}

/* 
 *	this saves only the string position at location n_save
 *
 *	could also save the entire string profile -- y_current[n]
 *	if desired
 */
save_y()
{
	fprintf(f_save,"%g\t%g\n",t,y_current[n_save]);
	fflush(f_save);
	
	return;
}