/* 	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		
#define    g              1.0     	
#define    A              1.0	     
#define    TipodeBote	  3	     
#define    w              0.9983	    
#define    r              0.80	      	  	    					
#define    nu0            1e-2
#define    Tolerancia     0.0005	
#define    epsilon        1e-5	
//Aca se ponen los valores que se desean probar

	   
short int  abortar,hayraiz; 			
long int   MiniBotes;		
char	   Nombre[30];				
double     PeriodoDeLaBase,TiempoDePrueba,VelocidadDePrueba;			
double     AlturaDeLaParticula,V,vz,TiempoDeVuelo,TiempoAbsoluto; 
double	   Tsemilla,TiempoFinal; 		
double	   nu, t00;
double	   dF1dx,dF2dx,dF1dy,dF2dy, Est;        
void EvolucionDelSistema();
double DiferenciaDeAltura(double tiempo);
void Newton(); 				
void CalcularTiempoDeVuelo();		
void actualizar();				
void Diferencial(double Phi,double Vi);		
void Estabilidad(double phi,double vi);		

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

 int main(){
 printf("\nSe usan coeficiente r = %f y la frecuencia w = %f. Solo se analiza el bote triple. \n\n", r, w);
 printf("Ingresar Tiempo a probar:\n");
 scanf("%lf",&TiempoDePrueba);
 printf("Se ingreso el tiempo T = %f\n\n",TiempoDePrueba);
 printf("Ingresar la Velocidad a probar:\n");
 scanf("%lf",&VelocidadDePrueba);
 printf("Se ingreso la Velocidad V = %f\n\n",VelocidadDePrueba);
 AlturaDeLaParticula   =A*sin(w*TiempoAbsoluto);
 Estabilidad(TiempoDePrueba*w,VelocidadDePrueba); 
 if(Est) printf("\nLos valores introducidos dan un bote estable\n\n");
 else    printf("\nLos valores introducidos dan un bote inestable\n\n");
 return 0; // Just to be polite! Jeje
 }


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

/*  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 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
}

/***********************************************************************/
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=1;						//estable
 else if (det>0  && fabs((-b+sqrt(det))/2)<1 && fabs((b+sqrt(det))/2)<1) Est=1;	//estable
 else Est= 0;									//inestable
 }
/***********************************************************************/