// Nota sobre el programa: // // Esta aplicacion busca obtener informacion para el sistema Bouncing-Ball // una vez que ha pasado una cierta cantidad de botes. // // En particular, esta version asigna una probabilidad a un determinado tipo de bote. // // Se omitiran las tildes #include #include #include #define pi M_PI // Constante pi #define Relajacion 5000 // Botes de relajacion #define BotedeInteres 3 // Tipo de bote que se busca caracterizar #define RamaAnalizada 1.0 // Rama que interesa #define Revision 50 // Veces que se verifica el tipo de bote #define CondicionesTotales 100000 // Condiciones iniciales vistas #define g 1.0 // Aceleracion gravitatoria #define A 1.0 // Amplitud de la base #define PeriodoMaximo 10 // Tipo de bote maximo para considerar periodico el mov. #define w 1.00 // Ultima frecuencia de simulacion #define r0 0.753 // Primer coef. de restitucion a simular #define rf 0.803 // ultimo coef. de restitucion a simular #define dr 0.002 // Paso numerico entre coeficientes #define Cuociente0 1e-2 // Cuociente para discriminar entre botes cortos y normales #define Tolerancia 1e-4 // Dos cantidades que difieran en "Tolerancia" seran consideradas iguales 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 long int ContadorDeBotes, MiniBotes; // Botes que han ocurrido y Cantidad de botes cortos char Nombre[30]; // Servira para guardar archivos de distintos nombres double r; // Coeficiente de restitucion de la particula double PeriodoDeLaBase; // Frecuencia y Periodo de la base double Probabilidad; // Probabilidad de hayar el bote de interes double AlturaDeLaParticula,V,vz; // Variables del movimiento en un golpe double TiempoDeVuelo,TiempoAbsoluto; // Variables del movimiento en un golpe double Tsemilla,TiempoFinal; // Tiempo que recibe y entrega el metodo Newton() double Cuociente, t00; // void Escribir(); // Esta funcion crea o abre el archivo donde se guardaran los datos void EvolucionDelSistema(); // Hace que la particula pase por una etapa de relajacion void Iniciar(); // Da valores iniciales a las alturas, tiempos y velocidades double DiferenciaDeAltura(double tiempo); // Da la diferencia de altura entre la particula y la base 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 void GuardarDatos(); // Almacena los datos de interes void Desarrollo(); // Es el desarrollo principal del programa /***********************************************************************/ int main(){ srand48(120); int Condicion; // Contador de condiciones iniciales printf("\nSe inicia el programa\n\n",r,w); sprintf(Nombre,"ProbabilidadFreq%1.2f.dat",w); // Nombre del archivo Escribir(); //Se escribe un mensaje inicial for(r=r0;r<=rf;r+=dr) // Se barren coeficientes de restitucion { Probabilidad = 0; // La probilidad de encontrar el punto triple es a priori es cero PeriodoDeLaBase = 2*pi/w; // Se establece el periodo de la base for(Condicion=0;Condicion<=CondicionesTotales;Condicion++) Desarrollo(); Probabilidad=Probabilidad/CondicionesTotales; GuardarDatos(); }printf("\n\n Se termino la ejecucion del programa\n\n"); return 0; // Se informa que el programa termino sin problemas su ejecucion } /***********************************************************************/ void Escribir(){ Archivo=fopen(Nombre,"wt"); // Se crea el archivo (o si se abre se borra el contenido) fprintf(Archivo," # Este archivo contiene la probabilidad de hayar un bote determinado.\n"); fprintf(Archivo," # r P\n"); fclose(Archivo); } /***********************************************************************/ /* 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*drand48()/w; // Iniciamos TiempoAbsoluto entre 0 y 2pi/w AlturaDeLaParticula = A*sin(w*TiempoAbsoluto) + 0.01;// Altura de la particula V = 20*drand48(); // Rapidez de la particula } /***********************************************************************/ void CalcularTiempoDeVuelo() { double TiempoAuxiliar,DeltaTiempo,AlturaParticulaAUX,AlturaBaseAUX; 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 // Estas teraciones 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 'Relajacion' botes contra el piso ContadorDeBotes = 0; MiniBotes = 0; Iniciar(); while(ContadorDeBotes<=Relajacion) { 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) {TipodeBote=0; break;} actualizar(); Cuociente = (TiempoAbsoluto-t00)/PeriodoDeLaBase; if(Cuociente