// Nota sobre el programa: // Intencionalmente se omitiran las tildes. Este programa genera un archivo de datos // para graficar las trayectorias de vuelo // #include #include #include #define pi M_PI // Numero pi #define Relajacion 5000 // La particula choca "Relajacion" 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.970 //Primera frecuencia revisada #define wf 0.971 //Ultima frecuencia revisada #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 ha encontrado el bote triple 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 w, PeriodoDeLaBase; // Caracteristicas del movimiento de la base cuando la amplitud es fija double dw; //Delta-frecuencia double AlturaDeLaParticula,V,VANT,z,vz; // Variables del movimiento en un golp 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 Escribir(); // Esta funcion crea o abre el archivo donde se guardaran los datos 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 GuardarDatos(); // Almacena datos del movimiento periodico void Desarrollo(); // Es el desarrollo principal del 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 dw = 1e-4; // Primer valor para dw sprintf(Nombre,"ApareceR%1.3f.dat",r); // Archivo de datos Escribir(); 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) fprintf(archivo," # Nota sobre el programa que genero este archivo:\n"); fprintf(archivo," # La primera columna guarda el coef. de rest. y la segunda la frecuencia\n"); fprintf(archivo," # a la que aparece el bote triple\n"); fclose(archivo); } /***********************************************************************/ void BarrerFrecuencias(){ for(w = w0;w <= wf; w += dw) // Se barren las frecuencias { Abortar = 0; // Abortar = 0 indica que no hay que Abortar TipodeBote = 0; // Valor inicial PeriodoDeLaBase = 2*pi/w; // Se establece el periodo de la base Desarrollo(); if(TipodeBote==3 && Contento>30) break; // La instruccion anterior imprime en pantalla el progreso del programa por cada 25 frecuencias revisadas }} /***********************************************************************/ // 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 '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) break; actualizar(); nu = (TiempoAbsoluto-t00)/PeriodoDeLaBase; if(nu