/* 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 #include #include #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 1 /* Primera frecuencia revisada */ #define wf 1 /* Ultima frecuencia revisada */ #define dw 1e-4 /* Delta-frecuencia */ #define r0 0.15 /* Primer coef de rest revisado */ #define rf 0.90 /* Ultimo coef de rest revisado */ #define dr 0.0005 /* 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 CoeficientesYaVistos; // 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; void Escribir(); // Esta funcion crea o abre el archivo donde se guardaran los datos void BarrerCoeficientesDeRestitucion(); // 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 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(w=w0;w<=wf;w+=dw){ // Se barren las frecuencias sprintf(Nombre,"Frecuencia%1.2f.dat",w); // Cuando no importan las ramas todo se guarda en "ResultadosCoefr.dat" Escribir(); //Se escribe un mensaje inicial BarrerCoeficientesDeRestitucion(); // 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," # La primera columna guarda el coef de restitucion\n\n"); fprintf(archivo," # La segunda columna guarda la velocidad de despegue de la particula\n"); fclose(archivo); } /***********************************************************************/ void BarrerCoeficientesDeRestitucion(){ CoeficientesYaVistos = 0; for(r = r0;r <= rf; r += dr) // 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(CoeficientesYaVistos%25==0) printf("\nSe revisa el coeficiente r= %f para la frecuencia W = %f\n",r,w); // La instruccion anterior imprime en pantalla el progreso del programa por cada 25 frecuencias revisadas CoeficientesYaVistos++; // 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(nuPeriodoMaximo ) && 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(){ archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.5f %8.5f\n",r,V); fclose(archivo); CalcularTiempoDeVuelo(); actualizar(); // Estos comando llevan el sistema al siguiente bote; se asume abortar=0 }