/*
 *	guitar.c
 *
 *	simulate a guitar string - model a perfectly flexible loss-less string
 *
 *	the output is the force on the bridge, which is just the component
 *		of the tension perpendicular to the top plate
 *	saved in the file "bridge_f"
 *
 *	N Giordano	4-3-04
 *
 *	just run it
 *
 *	init() sets initial values of the various parameters, including 
		n_string = number of discrete string elements
		length = length of string
		radius = radius of string
		frequency = desired resonant frequency
		beta = plucking point in fraction of string length
		amplitude = amplitude of pluck

		t_plot -- same the string profile every t_plot seconds
		t_end = time to run simulation

	init_string() initializes the string - gives it a plucked profile
		pluck it n_pluck units from the end

	can save the results to a file every t_plot seconds if desired

	also included is a routine to smooth out the kink at the plucking
		point, but this is not currently used
 */

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

#define	MAX	3001

#define PI	3.14159

#define	ON	0
#define OFF	1

double y_new[MAX],y[MAX],y_old[MAX];
double t,dt,dx;
double r2;	/* r2 = c^2 dt^2 / dx^2	*/
int n_string;	/* the string runs from n=0 to n_string - these are the ends*/
		/* bridge is at n = 0 so pluck near here	*/
double length,c,tension,density,radius;
double rho,mu;

double t_end;	/* length (in time) of the simulation	*/
double frequency;
double beta;	/* plucking point - fraction of sting length	*/
int n_pluck,n_smooth;
double amplitude;

double t_plot;	/* how often to plot	*/

char y_string_file[100];

main()
{

	init();
	save_y(0);
	pluck();
}

pluck()
{
	int i;
	double t_next_plot;
	int i_plot,test_flag;
	FILE *fp_force;

	fp_force = fopen("bridge_f","w");

	t_next_plot = t_plot;
	i_plot = 1;

	while(t < t_end) {
		/*	first propagate the string in the usual way	*/
		for(i = 1; i < n_string; i++) {
			y_new[i] = 2 * y[i] - y_old[i] 
				+ r2 * (y[i+1] + y[i-1] - 2.0 * y[i]);
		}

		/* now update current and old profiles		*/
		for(i = 0; i <= n_string; i ++) {
			y_old[i] = y[i];
			y[i] = y_new[i];
		}

		fprintf(fp_force,"%g\t%g\n",t,tension * (y[1] - y[0]));

		t += dt;
		t_next_plot -= dt;
		if(t_next_plot < 0.0) {
			save_y(i_plot);		
			++i_plot;
			t_next_plot = t_plot;
		}
	}

	fclose(fp_force);

	return;
}

init()
{
	double period;

	t = 0.0;
	n_string = 100; /* 1000; /* 1000; /* 300; /*100;	*/

	length = 0.65;
	rho = 8000;
	radius = 0.24 / 2.0;	/* seems to small	*/
	tension = 149;
	r2 = 1.0;
	frequency = 247;
	beta = 1.0 / 20.0;
	beta = 1.0 / 5.0;
	amplitude = 0.005;

	c = 2 * length * frequency;
	mu = PI * radius * radius * rho;
	tension = c * c * mu;

	sprintf(y_string_file,"y_string.");

	dx = length / n_string;
	dt = sqrt(r2 * dx * dx / (c * c));
	period = 1.0 / frequency;
	n_pluck = beta * n_string;

	t_plot = 0.05 / frequency;
	t_end = 1 * period;

	init_string();	/* initialize with a simple pluck		*/

	return;
}

/* initialize the string profile with a pluck at the bowing point	*/

init_string()
{
	int i;
	double amp;	/* amplitude of the pluck	*/
	double slope;

	amp = amplitude;

	slope = amp / (n_pluck * dx);
	for(i = 0; i <= n_pluck; i++) {
		y[i] = y_old[i] = i * dx * slope;
	}

	slope = amp / ((n_string - n_pluck) * dx);
	for(i = n_pluck; i <= n_string; i++) {
		y[i] = y_old[i] = amp - (i - n_pluck) * dx * slope;
	}

	return;
}

save_y(n)
int n;
{
	FILE *fp;
	char name[100];
	int i;
	double offset;

	offset = amplitude * n / 10;

	sprintf(name,"%s%d",y_string_file,n);
	fp = fopen(name,"w");
	for(i = 0; i <= n_string; i++) {
		fprintf(fp,"%g\t%g\n",i*dx,y[i] - offset);
	}
	fclose(fp);

	return;
}

/*	smooth out the string  near the pluck
 *	try to remove singularities associated
 *	with bow kinks
 *
 *	average y over +- n_smooth adjacent elements
 */
smooth_string()
{
	int i,j,n;
	double val;

	val = 0.0;
	n = 0;
	for(j = n_smooth; j < n_string - n_smooth; j++) {
		for(i = -n_smooth; i <= n_smooth; i++) {
			val += y[j + i];
			++n;
		}	
		y[j] = val / n;
	}

	return;
}