* 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 (4,FILE='ando_numerico.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.1 C DJ = 0.1 TETA = 30.0*cero TETA = TETA*RAD A = 0.01 B = -0.005 EPS = 0.1 T = -0.1 C VALORES DEL KERNEL PARA LA CONDICION INICIAL. F430 = A*DJ + B*DJ*DJ +EPS*T*DSQRT(DOS*DJ)*DCOS(TETA) WRITE (*,*) 'F430:',F430 CCCC LAZO SOBRE EL TIEMPO. CALCULO DE LA ORBITA ASOCIADA. Y(1) = DJ Y(2) = TETA C DELT = 1. SIGMAX = CERO DO 10 IN=0,10000 TIEMPO = CERO + DFLOAT(IN)*DELT CALL DERIVS (Y,DYDX) CALL RK4 (Y,DYDX,2,TIEMPO,DELT,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 SALIDA DE RESULTADOS. WRITE (4,*)tiempo, KNUM, HNUM C WRITE (4,*)tiempo, TETA/rad, DJ C 10 CONTINUE 20 CONTINUE C 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 ECUACIONES DE MOVIMIENTO. DYDX(1) = EPS*T*DSQRT(DOS*DJ)*DSIN(TETA) DYDX(2) = A + DOS*B*DJ + EPS*T*DCOS(TETA)/DSQRT(DOS*DJ) 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