/*
 *	fermi.c
 *
 *	simulate the Fermi-Pasta-Ulam problem (FPU)
 *	
 *	1D chain of masses (m=1) connected by linear springs (k=1) with
 *	a small (alpha) nonlinearity (quadratic in the force)
 *
 *	by N Giordano -- July, 2006
 */

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

#define MAX	100
#define PI	3.14159
#define YES	1
#define NO	0

double mode_energy();
double total_energy();
double mode_position();
double mode_velocity();

double x[MAX],x_old[MAX],x_new[MAX];
int n_particles;	/* number of particles in the chain	*/
double alpha;	/* size of nonlinearity - same for all masses	*/
double alpha_prime;	/* normalized alpha			*/

double beta,beta_prime;

double dt;	/* time step	*/
double r2;	/* used for iteration				*/
double t_max;	/* length of simulation	*/
double t_record;

int n_modes;

main(argc,argv)
int argc;
char *argv[];
{
	init();
	run();
}

init()
{
	int i;
	int mode_init;
	double omega_1;
	double amplitude,lambda;
	
	n_particles = 32;	/* there are 32 particles that can move -
				with the end particles connected to rigid
				boundary points				*/
	mode_init = 1;	/* mode that contains all of the initial energy	*/
			/* mode #1 is the lowest mode			*/
	amplitude = 1.0;	/* initial amplitude			*/
	amplitude = 10.0;

	/* string starts from rest					*/
	lambda = 2.0 / (mode_init);
	for(i = 1; i <= n_particles; i++) {
		x[i] = x_old[i] = amplitude * sqrt(2.0 / (n_particles+1)) * 
					sin(2.0 * PI * i / (lambda * (n_particles+1)));
	}
	x[0] = x_old[0] = 0.0;
	x[n_particles+1] = x_old[n_particles+1] = 0.0;


	omega_1 = PI / (n_particles+1);	/* frequency of the lowest mode	*/
	dt = 0.001 * 2.0 * PI / omega_1;
	dt = 1.0;
	dt = 0.3;
	dt = 0.02;
	dt = 0.05;
	r2 = dt * dt;
	t_max = 10000;
	t_max = 100000;
	n_modes = 10; /*5;	/* number of modes to track	*/
	t_record = 20 * dt;

	/* different values of beta and alpha of interest	*/
	beta = 1;	/* non-chaotic	*/
	beta = 3;	/* seems chaotic	*/
	beta = 0;
	beta = 0;
	beta = 3;	/* seems chaotic	*/
	beta = 6;
	beta = 1;	/* non-chaotic	*/
	beta = 0.3;

	beta_prime = beta * dt * dt;

	alpha = 1.0;
	alpha = 0.3;
	alpha = 0.6;
	alpha = 0.25;
	alpha = 0;
	alpha_prime = alpha * dt * dt;

	return;
}

run()
{
	int i;
	FILE *fp[MAX],*fp_total;
	FILE *fp_x;
	FILE *fp_force,*fp_alpha;
	FILE *fp_phase,*fp_phase_2;
	FILE *fp_a_mode[MAX];
	char s[100];
	double t,t_test;
	double f,non_lin_f;

	fp_total = fopen("total_e","w");
	fp_x = fopen("x_t","w");
	fp_force = fopen("force_t","w");
	fp_alpha = fopen("alpha_t","w");
	for(i = 1; i <= n_modes; i++) {
		sprintf(s,"mode.%d",i);
		fp[i] = fopen(s,"w");
		sprintf(s,"x_mode.%d",i);
		fp_a_mode[i] = fopen(s,"w");
	}
	fp_phase = fopen("phase_plot","w");
	fp_phase_2 = fopen("phase_plot_2","w");

	save_profile();

	t_test = t_record;
	t = 0;
	while(t < t_max) {
		for(i = 1; i <= n_particles; i++) {
			x_new[i] = r2 * (x[i+1] + x[i-1] - 2.0*x[i])
			+ alpha_prime*(pow(x[i+1]-x[i],2.0)-pow(x[i]-x[i-1],2.0))
			+ beta_prime*(pow(x[i+1]-x[i],3.0)-pow(x[i]-x[i-1],3.0))
			+ 2.0 * x[i] - x_old[i];
		}

		for(i = 1; i <= n_particles; i++) {
			x_old[i] = x[i];
			x[i] = x_new[i];
		}

		t += dt;
		t_test -= dt;
		if(t_test <= 0) {
			t_test = t_record;
			fprintf(fp_total,"%g\t%g\n",t,total_energy());
			fflush(fp_total);
			fprintf(fp_x,"%g\t%g\n",t,x[n_particles/3]);
			fflush(fp_x);
			for(i = 1; i <= n_modes; i++) {
				fprintf(fp[i],"%g\t%g\n",t,mode_energy(i));
				fflush(fp[i]);
				fprintf(fp_a_mode[i],"%g\t%g\n",
					t,mode_position(i));
				fflush(fp_a_mode[i]);
			}
			f = x[i+1] + x[i-1] - 2.0*x[i];
			fprintf(fp_force,"%g\t%g\n",t,f);
			non_lin_f = alpha*(pow(x[i+1]-x[i],3.0)-pow(x[i]-x[i-1],3.0));
			fprintf(fp_alpha,"%g\t%g\n",t,non_lin_f);
			fprintf(fp_phase,"%g\t%g\n",mode_position(1),mode_velocity(1));
			if(mode_cross(3) == YES) {
				fprintf(fp_phase_2,"%g\t%g\n",
					mode_position(1),mode_velocity(1));
			}
/*
			printf("t = %g\t\t",t);
			save_profile();
*/
		}
	}

	return;
}

double
total_energy()
{
	int i;
	double v,dx;
	double sum;

	sum = 0.0;
	for(i = 1; i <= n_particles+1; i++) {
		v = (x[i] - x_old[i]) / dt;
		dx = (x[i] + x_old[i] - x[i-1] - x_old[i-1]) / 2.0;
		sum += 0.5 * v * v + 0.5 * dx * dx;
	}

	return(sum);
}

double 
mode_energy(k)
int k;
{
	int i;
	double a,a_old,sum;

	a_old = a = 0.0;
	for(i = 1; i <= n_particles; i++) {
		a += x[i] * sin(i * k * PI / (n_particles+1));
		a_old += x_old[i] * sin(i * k * PI / (n_particles+1));
	}
/* not needed I think ??	*/
	a *= sqrt(2.0 / (n_particles+1));
	a_old *= sqrt(2.0 / (n_particles+1));

	sum = 0.5 * (a - a_old) * (a - a_old) / (dt * dt) + 
		0.5 * (a + a_old) * (a + a_old) * 
			sin(PI*k/(2*(n_particles+1))) * sin(PI*k/(2*(n_particles+1)));
	return(sum); 
}

save_profile(n)
int n;
{
	int i;
	FILE *fp;
	
	fp = fopen("prof.junk","w");
	for(i = 0; i <= n_particles+1; i++) {
		fprintf(fp,"%d\t%g\n",i,x[i]);
	}
	fclose(fp);
	
	return;
}

double 
mode_position(k)
int k;
{
	int i;
	double a,a_old,sum;

	a = a_old = 0.0;
	for(i = 1; i <= n_particles; i++) {
		a += x[i] * sin(i * k * PI / (n_particles+1));
		a_old += x_old[i] * sin(i * k * PI / (n_particles+1));
	}
	a *= sqrt(2.0 / (n_particles+1));
	a_old *= sqrt(2.0 / (n_particles+1));

	return((a + a_old) / 2.0); 
}

double 
mode_velocity(k)
int k;
{
	int i;
	double a,a_old,sum;

	a = a_old = 0.0;
	for(i = 1; i <= n_particles; i++) {
		a += x[i] * sin(i * k * PI / (n_particles+1));
		a_old += x_old[i] * sin(i * k * PI / (n_particles+1));
	}
	a *= sqrt(2.0 / (n_particles+1));
	a_old *= sqrt(2.0 / (n_particles+1));

	return((a - a_old) / dt); 
}

mode_cross(k)
int k;
{
	int i;
	double a,a_old;
	double sum;
	
	a = a_old = 0.0;
	for(i = 1; i <= n_particles; i++) {
		a += x[i] * sin(i * k * PI / (n_particles+1));
		a_old += x_old[i] * sin(i * k * PI / (n_particles+1));
	}

	if(a * a_old > 0) {
		return(NO);
	}
	else {
		return(YES);
	}
}