Programa baseado no método de Leapfrog para a simulação de uma corda real

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

void main()
{
	int i, n, j;

	// Corda C2
	//double l = 1.90, T = 750, M = 35e-3, s = M/l, L = 100, f = 16.0e3, e = 7.5e-6;

	// Corda C4
	//double l = 0.62, T = 670, M = 3.93e-3, s = M/l, L = 50, f = 32.0e3, e = 3.82e-5;

	// Corda C7
	//double l = 0.09, T = 750, M = 0.467e-3, s = M/l, L = 16, f = 96.0e3, e = 8.67e-4;

	//double *u, *g, *h, t, c = sqrt(T/s), dx = l/L, dt = 1.0/f, r = c*dt/dx;

	// Geral
	double l = 2;
	double *u, *g, *h, t, r, c = 300, dx = 0.01, dt = dx/(4*c), e = 7.5e-6, L = l/dx;
	char choice[20], file[20];
	FILE *arq;

	double N = L, tmax = 0.04;
	// Alocação de memória dinâmica dos ponteiros em forma de vetor

	u = malloc(L * sizeof(double));
	g = malloc(L * sizeof(double));
	h = malloc(L * sizeof(double));

	if(l != 1.9 && l != 0.62 && l != 0.09)
		if(e == 0)
			r = 1.0;
		else
			r = 1.0/4.0;


	///////////////////////////////////////
	// Condição de inicial: Estacionária //
	///////////////////////////////////////
	/*for(i = 0; i < L; i++)
	{
		u[i] = sin(M_PI*i/10.0);
		g[i] = sin(M_PI*i/10.0);
		h[i] = sin(M_PI*i/10.0);
	}

	u[0] = 0;
	g[0] = 0;
	u[L-1] = 0;
	g[L-1] = 0;*/


	///////////////////////////////////////////
	// Condição de inicial: Pacote gaussiano //
	///////////////////////////////////////////

	if(L == 16)
		for(i = 0; i < L; i++)
		{
			u[i] = exp(-pow(i-L/2,2)*0.5);
			g[i] = exp(-pow(i-L/2,2)*0.5);
			h[i] = exp(-pow(i-L/2,2)*0.5);
		}
	else
		for(i = 0; i < L; i++)
		{
			u[i] = exp(-pow(i-L/2,2)/(5*5));
			g[i] = exp(-pow(i-L/2,2)/(5*5));
			h[i] = exp(-pow(i-L/2,2)/(5*5));
		}

	//for(i=1;i<L-1;i++)
	//	{
	//		h[i] = r*r*u[i+1] + 2*u[i]*(1-r*r) + r*r*u[i-1] - h[i];
	//	}


	// Rotina onda em corda ideal
	/*for(t = 0; t <= tmax; t+=dt)
	{
		for(i=1;i<L-1;i++)
		{
			g[i] = r*r*u[i+1] + 2*u[i]*(1-r*r) + r*r*u[i-1] - h[i];
		}

		for(i=0;i<L;i++)
		{
			h[i] = u[i];
			u[i] = g[i];
		}
		printf("\nset grid\nset title 'Tempo = %.5lf'\nset yrange [-1:1]\npl \'-' w lp pt 7 ps 0.7\n", t);
		for (i = 0; i < L; i++)
		{
			printf("%lf\t%lf\n", i*dx, u[i]);
		}
		printf("e\npause 0.1\n");
	}*/
	printf("set term gif animate\nset o 'Giordano.gif'\n");
	// Rotina onda em corda real
	for(t = 0; t <= tmax; t+=dt)
	{
		for(i=1;i<L-1;i++)
		{
			if(i == 1)
				g[i] = (2-2*r*r-6*e*r*r*N*N)*u[i] - h[i] + r*r*(1+4*e*N*N)*(u[i+1]+u[i-1]) - e*r*r*N*N*(u[i+2]+(-u[i]));
			else if (i == L-2)
				g[i] = (2-2*r*r-6*e*r*r*N*N)*u[i] - h[i] + r*r*(1+4*e*N*N)*(u[i+1]+u[i-1]) - e*r*r*N*N*((-u[i])+u[i-2]);
			else
				g[i] = (2-2*r*r-6*e*r*r*N*N)*u[i] - h[i] + r*r*(1+4*e*N*N)*(u[i+1]+u[i-1]) - e*r*r*N*N*(u[i+2]+u[i-2]);
		}

		for(i=0;i<L;i++)
		{
			h[i] = u[i];
			u[i] = g[i];
		}
		
		if(t > tmax/10.0)
		{
			printf("\nset grid\nset title 'Tempo = %.5lf'\nset xrange [0:%lf]\nset yrange [-1:1]\npl \'-' t 'Giordano' w l lw 2 \n", t, l-dx);
			for (i = 0; i < L; i++)
			{
				printf("%lf\t%lf\n", i*dx, u[i]);
			}
			printf("e\npause 0.000001\n");
		}
	}

}