/*****************************************************************************
 * Metodo de Euler Explicito
 * 
 * Eduardo Oda - jan/2008
 * 
 * Esta eh uma implementacao do Metodo de Euler:
 * 	| y_0
 * 	| y_{k+1}=y_k+h*f(x_k,y_k)
 * 	
 * para integrar um problema de Cauchy:
 *  	| dy/dx = f(x,y)
 *  	| y(x_0) = y_0
 *  	| x \in [x_0,x_m]
 *  	| y \in R^n
 * 
 * O programa recebe como parametros o intervalo [x_0,x_m], a condicao inicial
 * y_0, a quantidade de passos m e a dimensao do espaco n. O tamanho do passo
 * eh calculado como (x_m-x_0)/m. A funcao f deve ser fornecida diretamente no
 * codigo. Para isso, edite a funcao
 * 
 * 	double f(double x, double y){...}
 * 
 * A saida sao tres arquivos: eulern.out, eulern.plot e eulern.rfase, gravados
 * no mesmo diretorio onde o comando foi chamado. O arquivo euler.plot
 * sobrescreve arquivos existentes, mas os arquivos eulern.out e eulern.rfase
 * acrescentam os dados ao final do arquivo deixando uma linha em branco.
 * 
 * O arquivo eulern.out eh uma tabela com n+1 colunas, a primeira eh o ponto
 * x_k e as seguintes sao aproximacoes de y_i(x_k), ou seja y_k[i].
 * 
 * O arquivo eulern.rfase tem os dados para plotar um retrato de fase do
 * problema se n==2.
 *
 * O arquivo eulern.plot eh um script do gnuplot que plota um retrato de fase
 * se n=2. Para executar utilize o comando:
 * 	$ gnuplot eulern.plot
 *
 * Este arquivo ja é gravado como um cabecalho de shell script, assim, se vc
 * der permissao de execucao para ele, ele roda esse comando automaticamente.
 *
 * Para compilar este programa utilize o comando:
 * 	$ gcc -Wall -ansi -lm eulern.c -o eulern
 * 
 * E para rodar o programa:
 * 	$ eulern m n x_0 x_m y_0[1] ... y_0[n] 
 *
 * Este programa foi escrito no padrao ANSI, portanto deve rodar em qualquer
 * maquina com configuracoes usuais, mas foi testado em:
 * 
 * - Slackware Linux 12
 * 	$ uname -a
 * 	Linux honir 2.6.19.1 #7 Fri Oct 19 00:22:45 BRST 2007 i686 Mobile AMD
 * 	Sempron(tm) Processor 2800+ AuthenticAMD GNU/Linux
 *
 * Pre-requisitos:
 * 	- gcc 4.1.2 ou superior
 * 	- gnuplot 4.2 ou superior
 * 
 *****************************************************************************/

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

/*****************************************************************************
 * Declaração da Função Derivada
 *****************************************************************************/
double f(double x, double y[], int i){
	if(i==0)
		return(y[1]);
	else if(i==1)
		return(-y[0]);
	return(1);
}

int main(int argc,char *argv[]){
	int m, n, i, j;
	double x_0,x_m, h, x_k;
	FILE *out, *rfase, *plot;

	/* 
	 * Converte os parametros
	 */
	if(argc<6){
		printf("Modo de uso:\n\t$ eulern m n x_0 x_m y_0[1] ... y_0[n] \n\n");
		return(1);
	}

	m=atoi(argv[1]);
	n=atoi(argv[2]);
	double y_0[n], y_k[n], y_aux[n];
	x_0=atof(argv[3]);
	x_m=atof(argv[4]);
	for(i=0;i<n;i++)
		y_0[i]=atof(argv[i+5]);

	/* 
	 * Calcula do Passo de Integração 
	 */
	h=(x_m-x_0)/m;

	/* 
	 * Inicializa os variaveis e executa o metodo, mandado os resultados
	 * diretamente para a saida
	 */
	out=fopen("eulern.out","a");
	rfase=fopen("eulern.rfase","a");
	fprintf(out,"# x_k\t\t");
	for(i=0;i<n;i++)
		fprintf(out,"y_k[%d]\t\t",i);
	fprintf(out,"\n");

	x_k=x_0;
	fprintf(out,"%lf\t",x_k);
	for(i=0;i<n;i++){
		y_k[i]=y_0[i];
		fprintf(out,"\t%lf", y_k[i]);
	}
	fprintf(out,"\n");

	for(j=0;j<m;j++){
		fprintf(out,"%lf\t",x_k+h);
		for(i=0;i<n;i++){
			y_aux[i]+=h*f(x_k,y_k,i);
			fprintf(out,"\t%lf", y_aux[i]);
			fprintf(rfase,"\t%lf", y_aux[i]);
		}
		for(i=0;i<n;i++){
			y_k[i]=y_aux[i];
		}
		fprintf(out,"\n");
		fprintf(rfase,"\n");
		x_k+=h;
	}
	fprintf(out,"\n");
	fprintf(rfase,"\n");
	fclose(out);
	fclose(rfase);

	if(n==2 || n==3){
		/*
		 * Cria o script do gnuplot
		 */
		plot=fopen("eulern.plot","w");
		fprintf(plot,"#!/usr/bin/gnuplot -persist\n\n");
		if(n==2)
			fprintf(plot,"plot \"eulern.rfase\" title \"Retrato de fase\" with lines");
		else
			fprintf(plot,"splot \"eulern.rfase\" title \"Retrato de fase\" with lines");
		fclose(plot);

		/* 
		 * Roda o gnuplot caso receba o parametro -plot
		 */
		system("gnuplot -persist eulern.plot");
	}
	return(0);
}