/* Nota sobre el programa: * Esta aplicacion busca obtener informacion para el sistema Bouncing-Ball una vez que ha pasado una cierta cantidad de botes. * * 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 2000 // La particula choca "BotesIniciales" veces antes de que se guarden sus datos #define Revision 1 // 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.6 // Primera frecuencia revisada #define wf 1.8 // Ultima frecuencia revisada #define dw 1e-4 // Delta-frecuencia #define r0 0.15 // Primer coef de rest revisado #define rf 0.15 // 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; // Botes que han ocurrido y 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; // 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 nu, t00; // Coeficientes utiles void Escribir(); // Esta funcion crea o abre el archivo donde se guardaran los datos void BarrerFrecuencias(); // Esta funcion barre las frecuencias void EvolucionDelSistema(); // 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 bas 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 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,"GraficoCoef%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"); 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(nunu0) && (TipodeBote != 0)) {ContadorDeBotes2=1; while(ContadorDeBotes2<=TipodeBote*2 ) {if (abortar) break; GuardarDatosYAvanzar();ContadorDeBotes2++;} }} /***********************************************************************/ void GuardarDatosYAvanzar(){ archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.5f %8.5f\n",w,V); fclose(archivo); CalcularTiempoDeVuelo(); actualizar(); // Estos comando llevan el sistema al siguiente bote; se asume abortar=0 }