/*
 *	sound.c
 *
 *	N Giordano	July 2006
 *
 *	solve the coupled acoustic equations for v and p in a 1D
 *	geometry
 *
 *	drive one end as a piston and use a real acoustic impedance
 *	Z = p/v at the other
 *
 *	use staggered grids in space and time
 *		see book for details
 *
 *	for space:

 0   1/2   1   3/2   .....		N-1   N-1/2   N
 p    v    p    v                        p     v      p
wall                                                 wall

for the arrays containing v:  v[1] = v[1/2] etc.

	for time:
just update first p, then v, then p, etc.

 *
 *	init() -- initialize variables
		c = speed of sound
		rho = density of air
		length = length of pipe
		n_wall = number of spatial units along pipe
		Z = acoustic impedance of wall

	drive the piston at a given frequency

	The program is current configured to do a calculation at a single
	frequency, and record the pressure at the wall and the air velocity
	as a function of time

	However, it contains routines to scan the frequency and look at
	resonant behavior when the frequency passes through the resonant
	frequency of the pipe.
 */

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

#define MAX	3001
#define PI	3.14159

double v_speaker();	/* velocity at the speaker wall	*/
double propagate();

double v[MAX],v_old[MAX];
double p[MAX],p_old[MAX];

double length,dx,dt;
int n_wall;
double rho,c;
double Z;

double t_end;
double t_start_meas;

double f_start,f_end,df;
double f_current;

/* drive the end at i = 1/2 with this velocity 	*/
double v_0,frequency;
double t_speaker_off;	/* turn off the speaker	now	*/
double rand_max;

main()
{
	double pressure;
	FILE *fp_v_max;

	init();
	propagate();
/*
	fp_v_max = fopen("v_2_freq","w");
	f_current = f_start;
	while(f_current <= f_end) {
		pressure = propagate();
		fprintf(fp_v_max,"%g\t%g\n",f_current,pressure);
		f_current += df;
	}
	fclose(fp_v_max);
*/
}

init()
{
	int i;

	length = 1.0;
	rho = 1.3;
	c = 330;
	n_wall = 100;
	dx = length / n_wall;
	dt = dx / c;

	Z = 100;	/* does not work	*/
	Z = 6e3;	/*3e4;	/* 1e6;	*/
	v_0 = 1e-3;
	frequency = c / length;
	rand_max = pow(2.0,31.0) - 1.0;

	t_end = 200 / frequency;
	t_start_meas = t_end / 2;
	t_speaker_off = t_end;

	for(i = 0; i <= n_wall; i++) {
		p[i] = p_old[i] = 0.0;
		v[i] = v_old[i] = 0.0;
	}

	f_start = frequency / 4;
	f_end = 1.8 * frequency;
	df = (f_end - f_start) / 60;

	f_current = frequency / 4;

	srand(32);

	return;
}

double 
v_speaker(t)
double t;
{
	double val;
	
	if(t < t_speaker_off) {
		val = v_0 * sin(2 * PI * f_current * t);	

/*	want to add a white noise generator eventually	*/
/*	val = v_0 * (0.5 - rand() / rand_max);		
		printf("%g\t%g\n",t,val);		*/

	}
	else {
		val = 0.0;
	}

	return(val);
}

double
propagate()
{
	int i,n_v;
	double t,p_tmp,v_tmp,a1,a2;
	FILE *fp_p_wall,*fp_v;
	double p_max,v_max_2,tmp;

	fp_p_wall = fopen("p_wall", "w"); /* record p at the far wall */
	fp_v = fopen("v_signal", "w");
	n_v = n_wall / 2.5;	/* record v here	*/

	a1 = 1.0 - Z * dt / (rho * dx);
	a2 = 1.0 + Z * dt / (rho * dx);
	p_tmp = rho * c * c * dt / dx;
	v_tmp = dt / (rho * dx);
	t = 0;

	p_max = v_max_2 = 0.0;
	while(t < t_end) {
		fprintf(fp_p_wall,"%g\t%g\n",t,p[n_wall-1]);
		fprintf(fp_v,"%g\t%g\n",t,v[n_v]);

		if(p[n_wall-1] > p_max) p_max = p[n_wall-1];

/* update p first - only need to do the interior	*/
		for(i = 1; i < n_wall; i++) {
			p[i] += - p_tmp * (v[i] - v[i-1]);
		}

/* next update v - start with the interior	*/
		for(i = 1; i < n_wall - 1; i++) {
			v[i] += - v_tmp * (p[i+1] - p[i]);
		}
		v[0] = v_speaker(t);
		v[n_wall-1] += (a1 / a2) * v[n_wall-1]
			+ p[n_wall-1] * 2 * dt / (rho*dx*a2);

		if(t > t_start_meas) {
			tmp = 0.0;
			for(i = 1; i < n_wall - 1; i++) {
				tmp += v[i] * v[i];
			}
			if(tmp > v_max_2) v_max_2 = tmp;
		}

		t += dt;
	}

	fclose(fp_p_wall);
	fclose(fp_v);

	return(v_max_2);
}