/*
 *	program to simulate waves on a piano string
 *
 *	the string is perfectly flexible, and is excited by a blow from
 *	a hammer
 *
 *	need to simulate both the string and the hammer
 *
 *	one end of the string is held fixed (i.e., perfectly rigid)
 *	while the other is connected to a flexible bridge
 *	it is assumed that for the bridge the F = g1 v
 *
 *	add the effect of air damping with F(drag) ~ v^2 as in my book
 *
 *	also add the longitudinal vibration
 *
 *	N. Giordano     2-98
 *
 *	to compile	cc -o nonlinear nonlinear.c -lm
 *
 *	Warning: we do not include any mechanism for energy loss from
 *	the longitudinal mode, so its amplitude will continue to grow with
 *	time.  To be more realistic we should include some loss mechanisms.
 */

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

#define		YES	1
#define		NO	0

/* this program follows the same pattern as the one which treats only
the transverse vibrations - so only the main differences are documented here */
double y_new[2000],y_old[2000],y_current[2000]; /* transverse displacement   */
double w_new[2000],w_old[2000],w_current[2000]; /* longitudinal displacement */
double y_hammer,v_hammer,K,p,m_hammer,mu,dm,tension,length;
int hammer_contact;
double hammer_profile[40];
int hammer_size;
double t;
double x_save,x_hammer;
int n_max,n_save,n_hammer;
double dx,dt,c,r2;
double young,area;		/* Young's modulus and cross-sectional area */
double dt_lon,c_lon,r2_lon;	/* parameters for the longitudinal wave
					equation	*/
int n_lon;
double imped_y,drag_coeff,radius;	/* bridge impedance, and
						air drag coefficient	*/
FILE *f_save_1,*f_save_2,*f_save_3;

main()
{
	int i;
	int i_loop = 30;
	long int i_end = 10024;

	i_loop = 1;

	init();
	while(1) {
		for(i = 0; i < i_loop; i++) propagate();
		if(t > .02) break;
		save_y();
	}
	fclose(f_save_1);
	fclose(f_save_2);
	fclose(f_save_3);
}

init()
{
	int i;
	                   	/* all units are mks	*/
	m_hammer = 3e-3;	/* hammer mass		*/
	mu = 6e-3;		/* mass per unit length of string	*/
	tension = 670;		/* tension		*/
	length = 0.84;		/* string length for G3	*/
	length = 0.62;		/* string length for C4 (middle C)	*/
	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;
	young = 2.0e11;		/* Young's modulus	*/
	
	n_max = 100;
	c = sqrt(tension / mu);
	x_hammer = length * 6.0 / 7.0;	/* hammer strike point	*/
	t = 0.0;

	imped_y = 1e6;  /*1000;		/* mks	*/
	radius = 5e-4;		/* string radius	*/
	area = 3.14 * radius * radius;

	dx = length / n_max;
	dt = dx / c;

	n_max = length / dx;
	c_lon = sqrt(young * area / mu);
	dt_lon = dx / c_lon;
	n_lon = dt / dt_lon;
	dt_lon = dt / n_lon;
	c_lon = dx / dt_lon; 

	dm = mu * dx;
	n_hammer = x_hammer * n_max / length;

	/* an empirical guess at the shape of the hammer
		for the value of dx used here, this hammer interacts with the
		string over a distance of a little less than 1 cm

		hammer_profile[0] is fraction of the hammer force which acts
			at the center of the hammer, etc. for hammer_profile[1]

		that the sum of these numbers should be unity
		(the values for indices !=0 are used twice)
	*/
	hammer_profile[0] = 0.18;
	hammer_profile[1] = 0.17;
	hammer_profile[2] = 0.13;
	hammer_profile[3] = 0.075;
	hammer_profile[4] = 0.035;
	hammer_size = 4;
	
	drag_coeff = 1.3 / (3.14 * radius * 8e3);
	drag_coeff = 0.0;
	x_save = 0.90 * length;
	n_save = x_save * n_max / length;
	f_save_1 = fopen("y_non.new","w"); /* transverse displacement */
	f_save_2 = fopen("b_non.new","w"); /* transverse force on the bridge  */
	f_save_3 = fopen("w_non.new","w"); /* longitudinal force on the bridge*/

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

propagate()
{
	int i,j;
	double force,velocity;
	double dc2;

	r2 = dt * dt * c * c / (dx * dx);
	r2_lon = dt_lon * dt_lon * c_lon * c_lon / (dx * dx);
	dc2 = dt_lon * dt_lon * (c_lon * c_lon - c * c) / (2.0*dx*dx*dx);

	y_new[0] = y_current[0] + 
		dt * tension * (y_current[1] - y_current[0]) / (dx * imped_y);

/* first update the transverse displacement	*/
	for(i = 1; i <= n_max - 1; i++) {
		velocity = y_current[i] - y_old[i];
		y_new[i] = 2 * (1-r2) * y_current[i] - y_old[i]
			+ r2 * (y_current[i+1] + y_current[i-1])
			- drag_coeff * fabs(velocity) * velocity;
	}

/* 
 *	next update the longitudinal displacement n_lon times	
 *	because the longitudinal speed is much faster than the transverse speed
 *	the longitudinal time step must be much smaller to maintain 
 *	numerical stability	
 */
	for(j = 0; j < n_lon; j++) {
		for(i = 1; i <= n_max - 1; i++) {
			w_new[i] = 2 * (1-r2_lon) * w_current[i] - w_old[i]
			+ r2_lon * (w_current[i+1] + w_current[i-1]);
			w_new[i] += dc2 * (y_current[i+1]-y_current[i-1])
			* (y_current[i+1] + y_current[i-1] - 2*y_current[i]);
		}
		for(i = 0; i <= n_max - 1; i++) {
			w_old[i] = w_current[i];
			w_current[i] = w_new[i];
		}
	}

		/* hammer is extremely narrow - width = dx	*/
		/* 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] += hammer_profile[0] * dt * dt * force / dm;  
		for(i = 1; i <= hammer_size; i++) {
			y_new[n_hammer+i] += 
				hammer_profile[i] * dt * dt * force / dm;  
			y_new[n_hammer-i] += 
				hammer_profile[i] * 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 = 0; i <= n_max - 1; i++) {
		y_old[i] = y_current[i];
		y_current[i] = y_new[i];
	}

	t += dt;

	return;
}

save_y()
{
	fprintf(f_save_1,"%g\t%g\n",t,y_current[n_save]); 
						/* transverse displacement */
	fflush(f_save_1);

	fprintf(f_save_2,"%g\t%g\n",t,(y_current[1]-y_current[0]) * tension/dx);
				/* transverse force on the bridge	*/
	fflush(f_save_2);
	
	fprintf(f_save_3,"%g\t%g\n",t,w_current[1]-w_current[0]);
				/* longitudinal displacement at the end 
				also proportional to the longitudinal force
				on the bridge				*/
	fflush(f_save_3);
	
	return;
}