// Nota sobre el programa: // Esta aplicacion busca obtener VER como aparece el bote triple // En particular se revisa la aparicion del bote triple // // Intencionalmente se omitiran las tildes. // #include #include #include #define pi M_PI // Numero pi #define BotedeInteres 3 // 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 w0 1.00 // Primera frecuencia revisada #define wf 1.00 // Ultima frecuencia revisada #define dw 1.00 // Delta-frecuencia #define r0 0.80 // Primer coef de rest revisado #define rf 0.80 // Ultimo coef de rest revisado #define dr 0.01 // Delta-Coef_de_roce #define Tolerancia 1e-5 // Tolerancia #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 short int Estable; // Estable = 1 significa que el bote es estable. long int ContadorDeBotes, MiniBotes; // Botes que han ocurrido y Cantidad de botes cortos long int Incremento; // Se usa para guardar solo ciertos resultados 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 double Alfa; // double AlturaDeLaParticula,V,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,V0; // Coeficiente util y C.I. double dF1dx,dF2dx,dF1dy,dF2dy; // Variables que se usan para determinar la estabilidad void CrearArchivos(); // Crea los archivos de escritura void EvolucionDelSistema(); // Hace que la particula colisione hasta relajarse void Iniciar(); // Da valores iniciales a las alturas, tiempos y velocidades double DiferenciaDeAltura(double tiempo); // Funcion que entrega la diferencia de altura entre la bola 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 Diferencial(double Phi,double Vi); // Encuentra una aproximacion del diferencial // de la funcion que lleva de un punto a otro punto igual 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(){ for(r=r0;r<=rf;r+=dr){ // Se barren coeficientes de restitucion for(w = w0;w <= wf; w += dw) {Abortar=0;PeriodoDeLaBase=2*pi/w; printf("\nPara el coeficiente r = %f se analiza la frecuencia W = %f\n",r,w); CrearArchivos(); for(V0=0.1;V0<=6;V0+=0.01) { for(t00=0;t00<=2*pi/w;t00+=0.001){Incremento=0;Iniciar();GuardarDatos();}} }} printf("\n\n Se termino la ejecucion del programa\n\n"); return 0;} /***********************************************************************/ void CrearArchivos(){ sprintf(Nombre,"R%1.2fF%1.4fMayorEstable.dat",r,w) ; archivo=fopen(Nombre,"wt");fclose(archivo); sprintf(Nombre,"R%1.2fF%1.4fMayorInestable.dat",r,w); archivo=fopen(Nombre,"wt");fclose(archivo); sprintf(Nombre,"R%1.2fF%1.4fMenorEstable.dat",r,w) ; archivo=fopen(Nombre,"wt");fclose(archivo); sprintf(Nombre,"R%1.2fF%1.4fMenorInestable.dat",r,w); archivo=fopen(Nombre,"wt");fclose(archivo); } /***********************************************************************/ 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 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() { TiempoAbsoluto = t00; vz = A*w*cos(w*TiempoAbsoluto); // Rapidez del suelo AlturaDeLaParticula = A*sin(w*TiempoAbsoluto) + 1e-6;// Altura de la particula V = V0; } /***********************************************************************/ 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 Diferencial(double phi,double vi) {double dtdx,dtdy,tsiguiente,tanterior; // Se quiere calcular la derivada de la funcion del tiempo c/r al tiempo absoluto TiempoAbsoluto=(phi+epsilon)/w; AlturaDeLaParticula=A*sin(phi); V=vi; CalcularTiempoDeVuelo(); tsiguiente=TiempoDeVuelo; TiempoAbsoluto=(phi-epsilon)/w; CalcularTiempoDeVuelo(); tanterior=TiempoDeVuelo; dtdx=(tsiguiente-tanterior)/(2*epsilon); /*derivada temporal c/r al tiempo absoluto*/ // Se quiere calcular la derivada de la funcion del tiempo c/r a la velocidad TiempoAbsoluto=phi/w; AlturaDeLaParticula=A*sin(phi); V=vi+epsilon; CalcularTiempoDeVuelo(); tsiguiente=TiempoDeVuelo; V=vi-epsilon; CalcularTiempoDeVuelo(); tanterior=TiempoDeVuelo; dtdy=(tsiguiente-tanterior)/(2*epsilon); /*derivada temporal c/r a la velocidad*/ dF1dx=1+w*dtdx; dF1dy=w*dtdy; dF2dx=-(1+r)*A*w*sin(phi)*(1+w*dtdx)+r*g*dtdx; dF2dy=-(1+r)*A*w*sin(phi)*w*dtdy-r + r*g*dtdy; // Ahora conviene que TiempoDeVuelo sea el valor que realmente le corresponde TiempoAbsoluto=phi/w; AlturaDeLaParticula=A*sin(phi); V=vi; CalcularTiempoDeVuelo(); } /***********************************************************************/ void Estabilidad(double Phi,double Vi){ double dF1dx2,dF2dx2,dF1dy2,dF2dy2,b,c,det,P11,P12,P21,P22; int ContadorEst; // Primero se calcula el diferencial para el argumento de esta funcion Diferencial(Phi,Vi); dF1dx2=dF1dx; dF2dx2=dF2dx; dF1dy2=dF1dy; dF2dy2=dF2dy; Phi=Phi+TiempoDeVuelo*w; Vi= (1+r)*A*w*cos(Phi) - r*(Vi - TiempoDeVuelo*g); for(ContadorEst =2;ContadorEst<=BotedeInteres;++ContadorEst){ Diferencial(Phi,Vi); P11=dF1dx*dF1dx2+dF1dy*dF2dx2; P12=dF1dx*dF1dy2+dF1dy*dF2dy2; P21=dF2dx*dF1dx2+dF2dy*dF2dx2; P22=dF2dx*dF1dy2+dF2dy*dF2dy2; dF1dx2=P11; dF2dx2=P21; dF1dy2=P12; dF2dy2=P22; Phi=Phi+TiempoDeVuelo*w; Vi= (1+r)*A*w*cos(Phi) - r*(Vi - TiempoDeVuelo*g); } b = -(dF1dx2+dF2dy2); c = dF1dx2*dF2dy2 - dF1dy2*dF2dx2; det= b*b-4*c; if (det<=0 && fabs(-b/2)<1) Estable = 1; else if (det>0 && fabs((-b+sqrt(det))/2)<1 && fabs((b+sqrt(det))/2)<1) Estable = 1; else Estable = 0; } /***********************************************************************/ void GuardarDatos(){ int Bote; double T11, V1; for(Bote=1;Bote<=BotedeInteres;Bote++){ CalcularTiempoDeVuelo(); actualizar();if(Abortar) break;} if(Abortar==0) {T11 = TiempoAbsoluto; V1 = V; Estabilidad(w*t00,V0); if(V0V1 && fmod(w*t00,2*pi)>fmod(w*T11,2*pi)) { if(Estable == 1 && Incremento%1 == 0){ sprintf(Nombre,"R%1.2fF%1.4fMenorEstable.dat",r,w); archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.3f %8.3f\n",V0,t00); fclose(archivo);} else if(Estable == 0 && Incremento%5 == 0){ sprintf(Nombre,"R%1.2fF%1.4fMenorInestable.dat",r,w); archivo=fopen(Nombre,"at"); fprintf(archivo,"%8.3f %8.3f\n",V0,t00); fclose(archivo);} } Incremento++; }}