/* 	Nota sobre el programa:

* Esta aplicacion busca obtener informacion para el sistema Bouncing-Ball una vez que ha pasado una cierta cantidad de botes.
*
*
* Se guardan datos como las rapideces terminales, las posiciones de colision, 
* tiempo de vuelo, tipo de bote y estabilidad (si el movimiento es periodico), entre otros.
*
*
* El programa es una modificacion al programa  "1D2.c",
* hecho por el profesor Patricio Cordero S. del Departamento de Fisica de la FCFM U. de Chile
*
* La informacion para cada coeficiente de restitucion "r" se encuentra en un archivo de nombre "ResultadosCoefr.dat"
* Aqui hay informacion bastate completa.
* Los datos separados por ramas estan en archivos llamados "RamaRCoefr.dat", donde R es el numero de la rama.
* Estos ultimos archivos se abren y usan, asi que si hay textos con ese nombre los datos nuevos y antiguos seran combinados.
*
*
* Intencionalmente se omitiran las tildes.
*/

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

#define    pi             M_PI		  /*  Numero pi						     				 */
#define    BotesIniciales 5000		  /*  La particula choca "BotesIniciales" veces antes de que se guarden sus datos 	 */
#define    Revision	  10		  /*  Veces que uno revisa que el tipo de bote sea coincidente			 	 */
#define    g              1.0     	  /*  Aceleracion de gravedad              		  		     	    	 */
#define    A              1.0	          /*  Amplitud del movimiento del suelo				         	         */
#define    PeriodoMaximo  100	          /*  Periodo maximo a revisar			   	    				 */
#define    w0             0.9	          /*  Primera frecuencia revisada            	  	    				 */	
#define    wf             1.1	          /*  Ultima frecuencia revisada          		    				 */
#define    dw             1e-4		  /*  Delta-frecuencia                    		    				 */
#define    r0             0.80	          /*  Primer coef de rest revisado            	  	    				 */	
#define    rf             0.80	          /*  Ultimo coef de rest revisado          		    				 */
#define    dr             0.10		  /*  Delta-Coef_de_roce                    		    				 */
#define    nu0            1e-2	          /*  Controla que no haya botes muy pequeños   		    			 */
#define    Tolerancia     1e-4	          /*  Dos valores cuya distancia sea menor a 'Tolerancia' seran vistos como iguales	 */
#define    epsilon        1e-5		  /*  Numero que sirve para aproximar derivadas		  			         */


FILE       *archivo;			  /* Respresentacion virtual del archivo de datos	     */

short int  abortar,hayraiz; 			// Indican si un bote o secuencia de botes son validos o no

int  	   TipodeBote; 				// Muestra el tipo de bote, ejemplo: TipodeBote = 1 => bote simple

long int   ContadorDeBotes, MiniBotes;		// La 1º variable indica los botes que han ocurrido y la 2º la cantidad de botes cortos

long int   FrecuenciasYaVistas;			// Frecuencias que se han ya revisado

char	   Nombre[30];				// Servira para guardar archivos de distintos nombres

double 	   r;					// Coeficiente de restitucion de la particula

double     w, PeriodoDeLaBase;			// Caracteristicas del movimiento de la base cuando la amplitud es fija

double     AlturaDeLaParticula,V,VANT,z,vz,TiempoDeVuelo,TiempoAbsoluto,energia; // Variables del movimiento en un golpe

double	   Tsemilla,TiempoFinal; 		// Tiempo que recibe y entrega el metodo Newton()

double	   nu, t00;

double	   dF1dx,dF2dx,dF1dy,dF2dy, Est;        // Variables que se usan para determinar la estabilidad del movimiento    

void Escribir();				// Esta funcion crea o abre el archivo donde se guardaran los datos

void BarrerFrecuencias();			// Esta funcion barre las frecuencias

void EvolucionDelSistema();			// Hace que la particula colisione una cierta cantidad de veces para luego analizar el movimiento

void Iniciar();					// Da valores iniciales a las alturas, tiempos y velocidades

double DiferenciaDeAltura(double tiempo);	// Funcion que entrega la diferencia de altura entre la particula y la base en un tiempo dado

void Newton(); 					// Este metodo entrega el tiempo "exacto" en que ocurre la colision

void CalcularTiempoDeVuelo();			// Calcula el tiempo de vuelo usando entre otras cosas el metodo Newton()

void actualizar();				// Actualiza alturas, tiempos y velocidades despues de cada choque

void DeterminarTipoDeBote();			// Esta funcion dice que tipo de bote hay en el caso, hace TipodeBote = 0 si el mov. no es periodico

void Diferencial(double Phi,double Vi);		// Encuentra una aproximacion del diferencial de la funcion que lleva de un punto a otro punto igual

void Estabilidad(double phi,double vi);		// Da la estabilidad de la funcion mencionada anteriormente

void GuardarDatosYAvanzar();			// Almacena datos del movimiento periodico
	
void Desarrollo();				// Es el desarrollo principal del programa


/***********************************************************************/

 int main(){

 srand48(100); // Deseamos iniciar las variables con cierto grado de aleatoriedad
 for(r=r0;r<=rf;r+=dr){ // Se barren coeficientes de restitucion
 sprintf(Nombre,"ResultadosCoef%1.2f.dat",r); // Cuando no importan las ramas todo se guarda en "ResultadosCoefr.dat"
 Escribir(); //Se escribe un mensaje inicial
 BarrerFrecuencias(); // Se llama a la funcion que barre las frecuencias
 } printf("\n\n Se termino la ejecucion del programa\n\n");
 return 0; // Just to be polite! Jeje
 }

/***********************************************************************/
 void Escribir(){

 archivo=fopen(Nombre,"wt"); // Se crea el archivo (o si se abre se borra el contenido)
 
 // A continuacion se guardara en el archivo un pequeño mensaje aclaratorio 
//,energia,TipodeBote,Est);
 
 fprintf(archivo," # 		Nota sobre el programa que genero este archivo:\n # ");
 fprintf(archivo,"Esta aplicacion busca obtener informacion para el sistema Bouncing-Ball ");
 fprintf(archivo,"una vez que ha pasado una cierta cantidad de botes.\n # ");
 fprintf(archivo,"El programa es una modificacion al programa  1D2.c, ");
 fprintf(archivo,"hecho por el profesor Patricio Cordero S. del Departamento de Fisica de la FCFM U. de Chile\n \n");
 fprintf(archivo," # La primera columna guarda la frecuencia de oscilaciones de la base\n\n");
 fprintf(archivo," # La segunda columna guarda la velocidad de despegue de la particula 'ahora'\n");
 fprintf(archivo," # La tercera columna guarda la velocidad de despegue anterior\n\n");
 fprintf(archivo," # La cuarta columna guarda la posicion en que ocurre el choque\n");
 fprintf(archivo," # La quinta columna guarda el tiempo de vuelo entre colisiones\n");
 fprintf(archivo," # La sexta  columna guarda la energia por unidad de masa\n");
 fprintf(archivo," # La septima  columna guarda el tipo de bote\n");
 fprintf(archivo," # La octava columna guarda una medida de la estabilidad de movimiento\n\n");
 fclose(archivo); 
 }
/***********************************************************************/

 void BarrerFrecuencias(){
 FrecuenciasYaVistas    = 0;
 for(w = w0;w <= wf; w += dw)   // Se barren las frecuencias
 {  abortar    = 0;		// abortar = 0 indica que no hay que abortar
    PeriodoDeLaBase  = 2*pi/w;  // Se establece el periodo de la base
    Desarrollo();		
    if(FrecuenciasYaVistas%25==0) printf("\nPara el coeficiente r= %f se analiza la frecuencia W = %f\n",r,w); 
    // La instruccion anterior imprime en pantalla el progreso del programa por cada 25 frecuencias revisadas
    FrecuenciasYaVistas++;	// Se agrega 1 al contador de frecuencias
 }}
/***********************************************************************/

/*  Esta es la funcion que dado un tiempo "tiempo" entrega la diferencia de altura entre la particula y el piso  */
double DiferenciaDeAltura(double tiempo)
{   double TiempoAuxiliar,AlturaDeLaBase;
    TiempoAuxiliar = TiempoAbsoluto + tiempo;
    AlturaDeLaBase   = A*sin(w*TiempoAuxiliar);
    return(AlturaDeLaParticula + V*tiempo - 0.5*g*tiempo*tiempo - AlturaDeLaBase);
}

/***********************************************************************/

/*Esta funcion encuentra el instante preciso en que ocurre la colision*/

// El funcionamiento de este algoritmo esta bien detallado en el apunte del profesor y por ende no se mencionara aca

void Newton()
{ double t0, t1, f0, f1, Af1;
  int mucho;
    hayraiz = 0;
    t0      = Tsemilla;
    f0      = DiferenciaDeAltura(t0);
    t1      = t0 + 1e-3/w;
    f1      = DiferenciaDeAltura(t1);
    mucho   = 0;
    do
    {if(f1!=f0) TiempoFinal = (t0*f1-t1*f0)/(f1-f0);
     t0     = t1;
     f0     = f1;
     t1     = TiempoFinal;
     f1     = DiferenciaDeAltura(t1);
     Af1    = fabs(f1);
     mucho++;
    }while(  (Af1>1e-8) && (mucho<100)  ) ;
    if(mucho<100) hayraiz=1;
}

/***********************************************************************/

void Iniciar()
{   abortar  			= 0; 				// Indicamos que no queremos abortar
    TiempoAbsoluto     		= (2*pi/w)*drand48();		// Iniciamos TiempoAbsoluto entre 0 y 2pi/w
    vz       			= A*w*cos(w*TiempoAbsoluto);	// Rapidez del suelo
    AlturaDeLaParticula      	= A*sin(w*TiempoAbsoluto) + 0.1;// Altura de la particula
    V       			= vz + 20.0*drand48()*g/w;	// Rapidez de la particula
}

/***********************************************************************/

void CalcularTiempoDeVuelo()
{ double TiempoAuxiliar,DeltaTiempo,AlturaParticulaAUX,AlturaBaseAUX;
    abortar= 0;			  // No queremos abortar
    DeltaTiempo    = 0.01/w;      // Difinimos un pequeño salto temporal
    TiempoAuxiliar   = 0.0;	  // Tiempo auxiliar que vale 0 cuando ocurre el choque anterior
    do

    {do                    /* iteraciones que dejan a la particula bajo el suelo */
     {TiempoAuxiliar   = TiempoAuxiliar + DeltaTiempo;			// Se actualiza el tiempo
      AlturaBaseAUX   = A*sin(w*(TiempoAbsoluto+TiempoAuxiliar));	// Se actualiza la altura de la base
      AlturaParticulaAUX   = AlturaDeLaParticula + V*TiempoAuxiliar - 0.5*g*TiempoAuxiliar*TiempoAuxiliar; 
     }while(AlturaParticulaAUX>AlturaBaseAUX);

     Tsemilla  = TiempoAuxiliar;
     Newton();       		// Se calcula el instante exacto de colision
     TiempoDeVuelo = TiempoFinal;
     DeltaTiempo = 0.5*DeltaTiempo;
     if(TiempoDeVuelo<1e-5)  MiniBotes++;
     else MiniBotes=0;             /* MiniBotes = veces que bote es muy breve*/
     if(  (MiniBotes>20) || (hayraiz==0) )
     { abortar=1;
       printf("Abortando, se han contado %4ld > 20 botes cortos o bien no hay raiz\n",MiniBotes);
       break;
     }
    }while(TiempoDeVuelo<0 && abortar==0);
}

/***********************************************************************/

void actualizar()
{TiempoAbsoluto  = TiempoAbsoluto + TiempoDeVuelo; // Se actualiza el tiempo
 vz    = A*w*cos(w*TiempoAbsoluto);		   // Se actualiza la rapidez del suelo
 AlturaDeLaParticula   =A*sin(w*TiempoAbsoluto);   // Se cambia la pos. de la particula
 V    = (1+r)*vz - r*(V - TiempoDeVuelo*g);	   // Nueva velocidad de despegue
 ContadorDeBotes++;				   // Se aumenta el contador de botes
}
/***********************************************************************/

 void EvolucionDelSistema(){ // Ocurren 'BotesIniciales' botes contra el piso
 ContadorDeBotes = 0;
 MiniBotes       = 0;
 Iniciar();
 while(ContadorDeBotes<=BotesIniciales)
 { CalcularTiempoDeVuelo();
 if(abortar) break;
 actualizar();
 }}

/***********************************************************************/
 void DeterminarTipoDeBote(){
 double Vauxiliar;
 t00       =  TiempoAbsoluto;
 Vauxiliar =  V;
 TipodeBote=1.0;
 while(TipodeBote<=PeriodoMaximo) 		//Se determina el tipo de bote que hay
  { CalcularTiempoDeVuelo();
    if(abortar) break;
    actualizar();
    nu = (TiempoAbsoluto-t00)/PeriodoDeLaBase;
    if(nu<nu0) {TipodeBote=0; break;} 		// Los botes muy cortos no son de gran interes
    if(fabs(V - Vauxiliar) < Tolerancia) break;
    ++TipodeBote;
 }}

/***********************************************************************/
void Diferencial(double phi,double vi)
 {double dtdx,dtdy,tsiguiente,tanterior;

 // Se quiere calcular la derivada de la funcion del tiempo c/r al tiempo absoluto
 TiempoAbsoluto=(phi+epsilon)/w;
 AlturaDeLaParticula=A*sin(phi);
 V=vi;
 CalcularTiempoDeVuelo();
 tsiguiente=TiempoDeVuelo;
 TiempoAbsoluto=(phi-epsilon)/w;
 CalcularTiempoDeVuelo();
 tanterior=TiempoDeVuelo;
 dtdx=(tsiguiente-tanterior)/(2*epsilon); /*derivada temporal c/r al tiempo absoluto*/
 
 // Se quiere calcular la derivada de la funcion del tiempo c/r a la velocidad
 TiempoAbsoluto=phi/w;
 AlturaDeLaParticula=A*sin(phi);
 V=vi+epsilon;
 CalcularTiempoDeVuelo();
 tsiguiente=TiempoDeVuelo;
 V=vi-epsilon;
 CalcularTiempoDeVuelo();
 tanterior=TiempoDeVuelo;
 dtdy=(tsiguiente-tanterior)/(2*epsilon); /*derivada temporal c/r a la velocidad*/
 dF1dx=1+w*dtdx;
 dF1dy=w*dtdy;
 dF2dx=-(1+r)*A*w*sin(phi)*(1+w*dtdx)+r*g*dtdx;
 dF2dy=-(1+r)*A*w*sin(phi)*w*dtdy-r + r*g*dtdy; 
 // Ahora conviene que TiempoDeVuelo sea el valor que realmente le corresponde
 TiempoAbsoluto=phi/w;
 AlturaDeLaParticula=A*sin(phi);
 V=vi;
 CalcularTiempoDeVuelo();
 }
/***********************************************************************/
 void Estabilidad(double Phi,double Vi){
 double dF1dx2,dF2dx2,dF1dy2,dF2dy2,b,c,det,P11,P12,P21,P22; 
 int ContadorEst;
 // Primero se calcula el diferencial para el argumento de esta funcion
 Diferencial(Phi,Vi);
 dF1dx2=dF1dx; 
 dF2dx2=dF2dx;
 dF1dy2=dF1dy;
 dF2dy2=dF2dy; 
 Phi=Phi+TiempoDeVuelo*w;
 Vi= (1+r)*A*w*cos(Phi) - r*(Vi - TiempoDeVuelo*g);
 for(ContadorEst =2;ContadorEst<=TipodeBote;++ContadorEst){
 Diferencial(Phi,Vi);
 P11=dF1dx*dF1dx2+dF1dy*dF2dx2;
 P12=dF1dx*dF1dy2+dF1dy*dF2dy2;
 P21=dF2dx*dF1dx2+dF2dy*dF2dx2;
 P22=dF2dx*dF1dy2+dF2dy*dF2dy2;
 dF1dx2=P11;
 dF2dx2=P21;
 dF1dy2=P12;
 dF2dy2=P22;  
 Phi=Phi+TiempoDeVuelo*w;
 Vi= (1+r)*A*w*cos(Phi) - r*(Vi - TiempoDeVuelo*g);
 }
 b = -(dF1dx2+dF2dy2);
 c = dF1dx2*dF2dy2 - dF1dy2*dF2dx2;
 det= b*b-4*c; 
 if (det<=0 && fabs(-b/2)<1) Est= (fabs(b/2));
 else if (det>0  && fabs((-b+sqrt(det))/2)<1 && fabs((b+sqrt(det))/2)<1) Est= ((fabs((-b+sqrt(det))/2) + fabs((b+sqrt(det))/2))/2); 
 else Est= 1.5;
 }
/***********************************************************************/
void Desarrollo()
{ 
 int RevBotes, TipodeBoteRespaldo,ContadorDeBotes2;
 RevBotes = 1;
 EvolucionDelSistema();
 DeterminarTipoDeBote(); TipodeBoteRespaldo=TipodeBote;

 while ( RevBotes <= Revision){
 DeterminarTipoDeBote();
 if (TipodeBoteRespaldo!=TipodeBote) {TipodeBote=0; break;}
 RevBotes++;} // La idea de las iteraciones anteriores es verificar si el tipo de bote es 'único' 
 VANT=V;  CalcularTiempoDeVuelo(); actualizar();

 // Ahora debemos guardar la inf para el movimiento que no mostro un periodo menor o igual al valor PeriodoMaximo:
 if(( TipodeBote==0 || TipodeBote>PeriodoMaximo ) && nu>=nu0)
 {TipodeBote=0; ContadorDeBotes2=1; 
 while(ContadorDeBotes2<=PeriodoMaximo) {if(abortar) break; GuardarDatosYAvanzar(); ContadorDeBotes2++;}} 

 //Caso en que el periodo existe y es menor a PeriodoMaximo
 if((nu>nu0)  && (TipodeBote != 0)) 
 {ContadorDeBotes2=1;  nu=nu/TipodeBote; if (fabs(nu-ceil(nu))<0.4)nu=ceil(nu); else nu=(int)nu;
 // nu se puede ver como la 'rama del diagrama de bifurcaciones'
 while(ContadorDeBotes2<=TipodeBote*2 ) {if (abortar) break; GuardarDatosYAvanzar();ContadorDeBotes2++;}
 }}
/***********************************************************************/

void GuardarDatosYAvanzar(){
 Est=1.5;
 energia=g*AlturaDeLaParticula+0.5*V*V;
 if (TipodeBote !=0 ){ // Se mide la estabilidad del movimiento
 Estabilidad(TiempoAbsoluto*w,V);				// Se calcula la estabilidad
 sprintf(Nombre,"Rama%dCoef%1.2f.dat",(int)nu,r); 
 archivo=fopen(Nombre,"at");
 fprintf(archivo,"%8.5f	8.5f	%8.3f	%8.3f	%8.3f	%8.3f	%8.3f	%d\n",w,V,VANT,AlturaDeLaParticula,TiempoDeVuelo,energia,TipodeBote,Est);
 fclose(archivo);
 sprintf(Nombre,"ResultadosCoef%1.2f.dat",r); // Cuando no importan las ramas todo se guarda en "ResultadosCoefr.dat"
 archivo=fopen(Nombre,"at"); 
 fprintf(archivo,"%8.5f	%8.5f	%8.3f	%8.3f	%8.3f	%8.3f	%d	%8.3f\n",w,V,VANT,AlturaDeLaParticula,TiempoDeVuelo,energia,TipodeBote,Est);
 fclose(archivo);
 }
 else {
 sprintf(Nombre,"ResultadosCoef%1.2f.dat",r); // Cuando no importan las ramas todo se guarda en "ResultadosCoefr.dat"
 archivo=fopen(Nombre,"at"); 
 fprintf(archivo,"%8.5f	%8.5f	%8.3f	%8.3f	%8.3f	%8.3f	%d	%8.3f\n",w,V,VANT,AlturaDeLaParticula,TiempoDeVuelo,energia,TipodeBote,1.5);
 fclose(archivo);}
 VANT=V;    // Se guarda el valor de la rapidez 
 CalcularTiempoDeVuelo(); actualizar(); // Estos comando llevan el sistema al siguiente bote; se asume abortar=0
 }