// Nota sobre el programa: // Intencionalmente se omitiran las tildes. Este programa genera una animación // del bote triple. // #include #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 0.971 //Primera frecuencia revisada #define wf 0.971 //Ultima frecuencia revisada #define dw 1e-4 //Delta-frecuencia #define r0 0.8 // Primer coef de rest revisado #define rf 0.8 //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; // 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 int Contento; // Veces que se se han buscado el bote triple int Graf; // Variable usada para graficar 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 golp double AlturaMaxima; // Maxima altura de la particula en el regimen asintotico double TiempoDeVuelo,TiempoAbsoluto; // Variables del movimiento en un golpe double Tiempito,Altura1,Altura2; // Datos que seran guardados en el archivo de trayectorias double Tsemilla,TiempoFinal; // Tiempo que recibe y entrega el metodo Newton() double nu, t00; // Coeficientes utiles void BarrerFrecuencias(); // Esta funcion barre las frecuencias void EvolucionDelSistema(); // Relaja al sistema 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 void Desarrollo(); // Es el desarrollo principal del programa int Dibujar(); // Diagrama en la pantalla void CerrarDibujo(); // Cierra el grafico y el programa /***********************************************************************/ int main(){ srand48(5); // Deseamos iniciar las variables con cierto grado de aleatoriedad for(r=r0;r<=rf;r+=dr){ // Se barren coeficientes de restitucion 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 BarrerFrecuencias(){ FrecuenciasYaVistas = 0; for(w = w0;w <= wf; w += dw) // Se barren las frecuencias { sprintf(Nombre,"TrayectoriaR%1.3fF%1.4f.dat",r,w); // Archivo de datos Abortar = 0; // Abortar = 0 indica que no hay que Abortar TipodeBote = 0; // Valor por defecto que toma la variable Contento = 0; // Contador de veces que se busca el bote triple PeriodoDeLaBase = 2*pi/w; // Se establece el periodo de la base do {Desarrollo();} while(Contento<100 && TipodeBote!=3); 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 AlturaDeLaParticula = A*sin(w*TiempoAbsoluto) + 0.1; // Altura de la particula vz = A*w*cos(w*TiempoAbsoluto); // Rapidez del suelo V = vz + 20.0*(1+drand48()/2)*g/w; // 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 /* 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(nu4){ pl_erase(); pl_pencolorname ("blue"); pl_fcircle(0,Altura1+0.1,0.1); pl_pencolorname ("black"); pl_fline(-1.5,Altura2,1.5,Altura2); pl_pencolorname ("green"); pl_fline(-1.5,-1,1.5,-1); pl_fline(-1.5,1,1.5,1);} else if (Choques<=3 && Altura1>=AlturaMaxima)AlturaMaxima=Altura1; else if (Choques==4 && Tiempito==0) Dibujar(); }}CerrarDibujo();}} Contento++;} /***********************************************************************/ int Dibujar(){ pl_parampl ("USE_DOUBLE_BUFFERING", "yes"); pl_parampl ("BITMAPSIZE","200x400"); Graf = pl_newpl ("X", stdin, stdout, stderr); pl_selectpl (Graf); if (pl_openpl () < 0) { fprintf (stderr, "No se pudo dibujar\n"); return 1;} pl_fspace (-2,-2,2,AlturaMaxima+1); pl_flinewidth (0.01); pl_filltype(1); pl_pencolorname ("blue"); pl_fillcolorname("red"); } /***********************************************************************/ void CerrarDibujo(){ pl_closepl(); pl_selectpl(0); pl_deletepl(Graf); exit(0); }