* ANDOYER_PELADO.F : RESOLUCION ANALITICA DEL ANDOYER * IMPLICIT REAL*8 (A-H,K-Z) REAL*8 Y(2),DYDX(2) COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE COMMON /COE/ A,B,EPS,T EXTERNAL DERIVS C ESPECIFICACION DE ARCHIVOS. OPEN (2,FILE='polin.dat',STATUS='UNKNOWN') OPEN (3,FILE='cubica.dat',STATUS='UNKNOWN') OPEN (4,FILE='soluciones.dat',STATUS='UNKNOWN') C INICIALIZACION DE VARIABLES GLOBALES. TWOPI = 8.0D0*DATAN(1.0D0) CERO = 0.0D0 UNO = 1.0D0 DOS = 2.0D0 UNO2 = UNO/DOS UNO3 = UNO/3.0D0 TRES2 = 3.0D0/2.0D0 RAD = TWOPI/360.0 PI = UNO2*TWOPI ERROR = 1.0D-13 C CONDICIONES INICIALES EN J Y TETA DJ = 0.01 TETA = 20.0*cero TETA = TETA*RAD write(*,*)'ondiciones iniciales', dj,teta C COEFICIENTES DEL HAMILTONIANO DE ANDOYER c A = -0.003 c B = 0.04 c EPS = 0.004 A = -0.00001 B = 0.01 EPS = 0.001 T = 0.01 C VALORES DEL KERNEL PARA LA CONDICION INICIAL. F430 = A*DJ + B*DJ*DJ + EPS*T*DSQRT(ABS(DOS*DJ))*DCOS(TETA) WRITE (*,*) 'F430:',F430 C FACTORES DE ESCALA. A1 = (4.0*EPS*T/B)**(0.333333333) A1ES = -3.0*B*A1*A1/4.0D0 H1ES = 0.25*B*(A1**4) WRITE(*,*)'H1ES', H1ES ALFA = A/A1ES C1 = F430/H1ES write(*,*)'A1ES', A1ES WRITE (*,*) 'a1, C1,ALFA:',a1,C1,ALFA C POLINOMIO CARACTERISTICO. DI = DJ/(A1**2) C DI = 0.287381654170384 POL = DOS*DI - (C1 + 3.0*ALFA*DI - 4.0*DI*DI)**2 C WRITE(*,*)'POL EN LA RAIZ',POL C STOP POL1 = (DOS-6.0*C1*ALFA) + DOS*(8.0*C1-9.0*ALFA**2)*DI POL1 = POL1 + 72.0*ALFA*DI**2 - 64.0*DI**3 POL2 = (16.0*C1 - 18*ALFA**2)+ 144.0*ALFA*DI -192*DI**2 C INVARIANTES DE LA ELIPTICA DE WEIERSTRASS G2 = UNO3*(8.0*C1 + 4.5*ALFA**2)**2 - 12.0*ALFA G3 = ((8.0*C1 + 4.5*ALFA**2)**3)/27.0D0 G3 = G3 - DOS*ALFA*(8.0*C1 + 4.5*ALFA**2) + 4.0D0 WRITE (*,*)'INVARIANTES G2,G3:',G2,G3 C CRITERIO DE LIBRACION (LIBRA SI DELTA < 0). DELTA = G2**3 - 27.0*G3*G3 WRITE(*,*)'VALOR DE DELTA:', DELTA C GRAFICO DEL POLINOMIO CARACTERISTICO c DIG = -1 PASO = 0.001 DO I = 0,1000 DIG = PASO*DFLOAT(I) POL = DOS*DIG - (C1 + 3.0*ALFA*DIG - 4.0*DIG*DIG)**2 WRITE(2,*) DIG, POL END DO C RESOLUCION NUMERICA DEL POLINOMIO CARACTERISTICO POL. DIN1= 0.15 11 POL = DOS*DIN1 - (C1 + 3.0*ALFA*DIN1 - 4.0*DIN1**2)**2 POL1 = (DOS-6.0*C1*ALFA) + DOS*(8.0*C1-9.0*ALFA**2)*DIN1 POL1 = POL1 + 72.0*ALFA*DIN1**2 - 64.0*DIN1**3 DIN2 = DIN1 - POL/POL1 POLI2 = DOS*DIN2 - (C1 + 3.0*ALFA*DIN2 - 4.0*DIN2**2)**2 if(Dabs(POLI2).gt.1.0E-10) then DIN1 = DIN2 go to 11 end if WRITE(*,*)'raiz del polinomio',DIN2 DI0 = DIN1 C GRAFICO DE LA ECUACION CUBICA PASO = 0.005 DO J = 0,1000 EE = -2.0+ PASO*DFLOAT(J) ECU = 4.0*ee**3-g2*ee-g3 WRITE(3,*) EE, ECU END DO C RESOLUCION DE LA RESOLVENTE CUBICA 4E^3 - G2E -G3 =0 EE_1=2.0 9 EE_2= EE_1 - (4.0*EE_1**3-g2*EE_1-g3)/(12.*EE_1**2-g2) if(abs(4.0*EE_2**3-g2*EE_2-g3).gt.1.0d-9) then EE_1 = EE_2 go to 9 end if WRITE(*,*) 'raiz de la cubica resolvente', EE_2 EE1 = EE_2 IF (DELTA.LT.CERO) THEN ETA = DSQRT(3.0*EE1*EE1 - G2/4.0) KA = DSQRT(UNO2 - 3.0*EE1/ETA/4.0) WRITE (*,*) 'ETA,KA',ETA,KA ELSE WRITE (*,*) 'CIRCULA!!' END IF C AMPLITUD DE LIBRACION. AP = DOS - DOS*(C1+3.0*ALFA*DI0-4.0*DI0*DI0)*(3.0*ALFA-8.0*DI0) AP = AP/4.0D0/ETA BP = -DOS*(3.0*ALFA-8.0*DI0)*(3.0*ALFA-8.0*DI0) + * DOS*(C1+3.0*ALFA*DI0-4.0*DI0*DI0)*8.0 BP = (BP + 24.0*(ETA-EE1))/24.0/ETA QQ = AP/(UNO - BP) WRITE (*,*) 'AP,BP',AP,BP WRITE (*,*) 'QQ,QQ [deg]',QQ,QQ/RAD C CENTRO DE LIBRACION. DIC = DI0 + UNO2*QQ write(*,*)'DIC', DIC C INTEGRAL ELIPTICA K(k). KA2 = KA*KA KELIP = UNO + 0.25*KA2 + 9.0*KA2*KA2/64.0 + 25.0*(KA2**3)/256.0 KELIP = KELIP + 1225.0*(KA2**4)/16384.0 + 3969.0*(KA2**5)/65536.0 KELIP = UNO2*PI*KELIP C PERIODO DE LIBRACION (ESCALEADO). PER = 4.0*KELIP/DOS/DSQRT(ETA) C PERIODO DE LIBRACION (EN TIEMPO FISICO). T0 = PER*(4.0/A1/A1/B) WRITE (*,*) 'PER,T0',PER,T0 stop CCCC LAZO SOBRE EL TIEMPO. CALCULO DE LA ORBITA ASOCIADA. c TETA0 = TETA Y(1) = DJ Y(2) = TETA C DELT = T0/100.0 SIGMAX = CERO DO 10 IN=0,1000 TIEMPO = CERO + DFLOAT(IN)*DELT CALL DERIVS (Y,DYDX) CALL RK4 (Y,DYDX,2,TIEMPO,DELT,Y,DERIVS) c CALL RK4 (Y,DYDX,2,TIEMPO,pass,Y,DERIVS) DJ = Y(1) TETA = Y(2) TETA = DMOD(TETA,TWOPI) IF (TETA.LT.CERO) TETA = TETA + TWOPI IF (TETA.GT.PI) TETA = TETA - TWOPI KNUM = DSQRT(ABS(DOS*DJ))*DCOS(TETA) HNUM = DSQRT(ABS(DOS*DJ))*DSIN(TETA) C SOLUCION ANALITICA. W1 = TWOPI*TIEMPO/T0 COU = BP/AP DIA = DI0 + UNO2*QQ -UNO2*QQ*DCOS(W1) DIA = DIA - UNO2*COU *Q2 + UNO2*COU *Q2*DCOS(DOS*W1) STETA = QQ*DSIN(W1)+ (UNO/8.0/DI0- UNO2*COU)*Q2*DSIN(DOS*W1) STETA = UNO2*TWOPI/PER/DSQRT(DOS*DIC)*STETA DIA1 = DI0 + UNO2*QQ -UNO2*QQ*DCOS(W1) STETA1 = QQ*DSIN(W1) STETA1 = UNO2*TWOPI/PER/DSQRT(DOS*DIC)*STETA1 C LA SOLUTION ANALITICA. DJA = DIA*A1*A1 T1CA = DASIN(STETA) IF (T1CA.LT.CERO) T1CA = T1CA + TWOPI KA = DSQRT(DOS*DJA)*DCOS(T1CA) HA = DSQRT(DOS*DJA)*DSIN(T1CA) DJA1 = DIA1*A1*A1 T1CA1 = DASIN(STETA1) IF (T1CA1.LT.CERO) T1CA1 = T1CA1 + TWOPI KA1 = DSQRT(DOS*DJA1)*DCOS(T1CA1) HA1 = DSQRT(DOS*DJA1)*DSIN(T1CA1) WRITE (4,100) KA1, HA1, KNUM, HNUM WRITE (*,*) KA1, HA1, KNUM, HNUM 10 CONTINUE C 100 FORMAT (1P9E15.6) C END C********************************************************************* C SUBRUTINAS SUBROUTINE DERIVS (Y,DYDX) C IMPLICIT REAL*8 (A-H,K-Z) REAL*8 Y(2),DYDX(2) COMMON /CTE/ TWOPI,CERO,UNO,DOS,UNO2,TRES2,RAD,ERROR,BASE COMMON /COE/ A,B,EPS,T C DJ = Y(1) TETA = Y(2) C HAMILTONIANO. F = A* DJ + B*DJ**2 + EPS*T*DSQRT(DOS*DJ)*DCOS(TETA) C ECUACIONES DE MOVIMIENTO. DYDX(1) = - EPS*T*DSQRT(abs(DOS*DJ))*DSIN(TETA) DYDX(2) = A + DOS*B*DJ + EPS*T*DCOS(TETA)/DSQRT(abs(DOS*DJ)) C DYDX(2) = DYDX(2) RETURN C END C******************************************************************* SUBROUTINE RK4 (Y,DYDX,IN,X,H,YOUT,DERIVS) C IMPLICIT REAL*8 (A-H,K-Z) PARAMETER (INMAX=10) REAL*8 Y(IN),DYDX(IN),YOUT(IN),YT(INMAX),DYT(INMAX),DYM(INMAX) C HH = 0.5D0*H H6 = H/6.0D0 XH = X + HH DO 11 I = 1,IN YT(I) = Y(I) + HH*DYDX(I) 11 CONTINUE CALL DERIVS (YT,DYT) DO 12 I = 1,IN YT(I) = Y(I) + HH*DYT(I) 12 CONTINUE CALL DERIVS (YT,DYM) DO 13 I = 1,IN YT(I) = Y(I) + H*DYM(I) DYM(I) = DYT(I) + DYM(I) 13 CONTINUE CALL DERIVS (YT,DYT) DO 14 I = 1,IN YOUT(I) = Y(I) + H6*(DYDX(I) + DYT(I) + 2.0*DYM(I)) 14 CONTINUE C RETURN C END