// Nota sobre el programa: // Intencionalmente se omitiran las tildes. Este programa genera un archivo de datos // con el modulo de la parte real de los autovalores de la funcion // #include #include #include #define pi M_PI // Numero pi #define BotesIniciales 1000 // 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.98 //Primera frecuencia revisada #define wf 1.097 //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 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 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 double dF1dx,dF2dx,dF1dy,dF2dy,det; // Variables que se usan para determinar la estabilidad double Real1,Real2; // Parte real de los Autovalores buscados double Imag; // Parte imaginaria de Autovalores buscados void Escribir1(); // Esta funcion crea o abre el archivo donde se guardaran los datos void Escribir2(); // 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 Diferencial(double Phi,double Vi); // Encuentra una aproximacion del diferencial de la funcion void Estabilidad(double phi,double vi); // Da la estabilidad de la funcion mencionada anteriormente 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 sprintf(Nombre,"RealR%1.3f.dat",r); // Archivo de datos Escribir1(); sprintf(Nombre,"ImagR%1.3f.dat",r); // Archivo de datos Escribir2(); 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 Escribir1(){ archivo=fopen(Nombre,"wt"); // Se crea el archivo (o si se abre se borra el contenido) fprintf(archivo," # Este programa genera un archivo de datos\n"); fprintf(archivo," # con la parte real de los autovalores de la funcion\n"); fprintf(archivo," # y tambien su modulo\n"); fprintf(archivo," # w R1 R2 fasb(R1) fasb(R2)\n\n"); fclose(archivo); } /***********************************************************************/ void Escribir2(){ archivo=fopen(Nombre,"wt"); // Se crea el archivo (o si se abre se borra el contenido) fprintf(archivo," # Este programa genera un archivo de datos\n"); fprintf(archivo," # con el modulo de la parte imaginaria de los autovalores de la funcion\n"); fprintf(archivo," # y tambien su modulo\n"); fprintf(archivo," # w Im\n\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 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<30 && 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(nu0 && fabs((-b+sqrt(det))/2)<1 && fabs((b+sqrt(det))/2)<1) {Real1=(-b+sqrt(det))/2; Real2=(b+sqrt(det))/2;} } /***********************************************************************/ void Desarrollo(){ if(Abortar==0){ double AlturaAux; int RevBotes, TipodeBoteRespaldo,ContadorDeBotes2; int Choques; RevBotes = 1; EvolucionDelSistema(); DeterminarTipoDeBote(); TipodeBoteRespaldo=TipodeBote; while ( RevBotes <= Revision){ DeterminarTipoDeBote(); if (TipodeBoteRespaldo!=TipodeBote) {TipodeBote=0; break;} RevBotes++;} // La idea de las iteraciones anteriores es verificar si el tipo de bote es 'único' TiempoAbsoluto = fmod(TiempoAbsoluto,PeriodoDeLaBase); nu=nu/TipodeBote; if (fabs(nu-ceil(nu))<0.4)nu=ceil(nu); else nu=(int)nu; // nu se puede ver como la familia vertical if(TipodeBote==3 && nu==1.0){ Estabilidad(w*TiempoAbsoluto,V); GuardarDatos(); }} Contento++;} /***********************************************************************/ void GuardarDatos(){ sprintf(Nombre,"RealR%1.3f.dat",r); archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.5f %8.5f %8.5f %8.5f %8.5f \n",w,Real1,Real2,fabs(Real1),fabs(Real2)); fclose(archivo); if(det<=0) { sprintf(Nombre,"ImagR%1.3f.dat",r); archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.5f %8.5f\n",w,Imag); fclose(archivo);} printf ("\nFalta el %2.2f por ciento\n",100*(wf-w)/(wf-w0)); }