/* 	Nota sobre el programa:

* Esta aplicacion busca obtener informacion para el sistema Bouncing-Ball una vez que ha pasado una cierta cantidad de botes.
*
* Intencionalmente se omitiran las tildes.
*/

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

#define    pi             M_PI		  /*  Numero pi						     				 */
#define    BotesEvol 	  30000	  	  /*  La particula choca "BotesEvol" veces antes de que se guarden sus datos 	 */
#define    Revision	  500		  /*  Veces que uno revisa que el tipo de bote sea coincidente			 	 */
#define    Condiciones	  200000	  /*  Condiciones iniciales revisadas		 	 */
#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.5	          /*  Primera frecuencia revisada            	  	    				 */	
#define    wf             1.1	          /*  Ultima frecuencia revisada          		    				 */
#define    dw             1e-4		  /*  Delta-frecuencia                    		    				 */
#define    r              0.95	          /*  Primer coef de rest revisado            	  	    				 */	
#define    Cuociente0     1e-2	          /*  Controla que no haya botes muy pequeños   		    			 */
#define    Tolerancia     1e-6		  /*  Dos valores cuya distancia sea menor a 'Tolerancia' seran vistos como iguales	 */


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

FILE       *archivo;			  	// Respresentacion virtual del archivo de datos				         
int  	   TipodeBote, FrecuenciasYaVistas;   	// Muestra el tipo de bote, ejemplo: TipodeBote = 1 => bote simple, y las frecuencias ya vistas
int	   RamaAnalizada;			// Rama en la que se analizan las frecuencias
long int   ContadorDeBotes;			// Indica los botes que han ocurrido 
char	   Nombre[30];				// Servira para guardar archivos de distintos nombres
double     w, PeriodoDeLaBase;			// Caracteristicas del movimiento de la base cuando la amplitud es fija
double     Probabilidad, CasosTotales;		// Variables para calcular la probabilidad
double     V,Fase; 				// Variables del movimiento en un golpe
double	   Cuociente;				//Cuociente que sirve para comparar cosas

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

void Escribir();				// Esta funcion crea o abre el archivo donde se guardaran los datos
void Avanzar();					// Actuliza posicion y velocidad post-choque
void DeterminarTipoDeBote();			// Esta funcion dice que tipo de bote hay en el caso, hace TipodeBote = 0 si el mov. no es periodico
void GuardarDatos();				// Almacena datos del movimiento periodico
void Desarrollo();				// Es el desarrollo principal del programa

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

 int main(){
 int Cuentas;
 srand48(115);
 printf("\nSe inicia el programa\n\n",r,w);
 for(RamaAnalizada=1; RamaAnalizada<=7;RamaAnalizada++){
 FrecuenciasYaVistas    = 0;
 sprintf(Nombre,"Rama%dProbabilidadCoef%1.2f.dat",RamaAnalizada,r); // Nombre del archivo
 Escribir(); //Se escribe un mensaje inicial
 for(w = w0;w <= wf; w += dw)   // Se barren las frecuencias
 {  Probabilidad = 0;           // La probilidad de encontrar el punto triple es a priori es cero 
    CasosTotales = 0;           // Para esta frecuencia aun no se analizan casos
    PeriodoDeLaBase  = 2*pi/w;  // Se establece el periodo de la base


    if(FrecuenciasYaVistas%25==0) {
    printf("** Para el coeficiente r = %f se analiza la rama # %d en la frecuencia W = %f **\n\n",r,RamaAnalizada,w);}
    // La instruccion anterior imprime en pantalla el progreso del programa por cada 25 frecuencias revisadas

    // Se recorren las condiciones iniciales deseadas
    for(Cuentas=0;Cuentas<=Condiciones;Cuentas++)
    { if(Cuentas%50==0) 
    {printf("Se ha revisado un %f por ciento de las c.i. para esta frecuencia\n",100.0*Cuentas/Condiciones);}
    Fase = 2*pi*drand48(); V = 20*drand48(); Desarrollo(); CasosTotales++;}

    //Se asigna una probabilidad a cada frecuencia:
    Probabilidad=Probabilidad/CasosTotales;
    GuardarDatos();
    FrecuenciasYaVistas++;}}	// Se agrega 1 al contador de 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 
 
 fprintf(archivo," # 		Nota sobre el programa que genero este archivo:\n # ");
 fprintf(archivo,"Esta aplicacion busca obtener la probabilidad de hallar un cierto tipo de mov. ");
 fprintf(archivo,"una vez que ha pasado una cierta cantidad de botes.\n # ");
 fprintf(archivo," # La primera columna guarda el coeficiente de restitucion");
 fprintf(archivo," # La segunda columna guarda la frecuencia de oscilaciones de la base\n\n");
 fprintf(archivo," # La tercera columna guarda la Probabilidad de estar en una determinada rama\n");
 fprintf(archivo," #  r           w          P\n");
 fclose(archivo);
 }

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

void Avanzar()
 {
 Fase += 2*w*V/g; 
 V = (1+r)*A*w*cos(Fase)+r*V;
 ContadorDeBotes++;
 }
/***********************************************************************/
 void DeterminarTipoDeBote(){
 double Vauxiliar,Fase0;
 Fase0       =  Fase;
 Vauxiliar =  V;
 TipodeBote=1.0;
 while(TipodeBote<=PeriodoMaximo) 		//Se determina el tipo de bote que hay
  { Avanzar();
    Cuociente = (Fase-Fase0)/(PeriodoDeLaBase*w);
    if(Cuociente<Cuociente0 || V<0 || fabs(V)<1.3) {TipodeBote=0; break;} 	
    if(fabs(V - Vauxiliar) < Tolerancia) break;
    ++TipodeBote;}
 }

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

void Desarrollo()
{int Rama;
 int RevBotes, TipodeBoteRespaldo;
 RevBotes = 1;
 ContadorDeBotes = 0;

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

 if(TipodeBote==1)	 Rama=Cuociente/TipodeBote;
 else 			 Rama=0;	 
 if (fabs(Rama-ceil(Rama))<0.4) Rama = ceil(Rama); 
 else Rama = (int)Rama;

 if(Rama==RamaAnalizada) Probabilidad++;
 
 }
/***********************************************************************/

void GuardarDatos(){
 archivo=fopen(Nombre,"at"); 
 fprintf(archivo,"%8.5f    %8.5f   %8.5f\n",r,w,Probabilidad);
 fclose(archivo);
 }