PROGRAM ORBITAL
      IMPLICIT REAL*8 (A-H,O-Z)
C...
      COMMON/POT/POT
      CHARACTER*8 POT
C...
      READ(*,*)POT
      READ(*,*)E,AL,AM1,AM2,TAU,RI
      IF(POT.EQ.'mgz     ')POT='MGZ     '
      IF(POT.EQ.'armonico')POT='ARMONICO'
      IF(POT.EQ.'gravit  ')POT='GRAVIT  '
      G=9.8D0
      DELTA=0.01D0
      NI=1
      NF=TAU/DELTA+1
      IF(POT.EQ.'MGZ     ')THEN
      AMU=AM1
      AMU1=AM1+AM2
      AM=1.D0
      CT=AM2*G
      ENDIF
      IF(POT.EQ.'ARMONICO')THEN
      AMU=AM1*AM2/(AM1+AM2)
      AMU1=AMU
      AM=AMU/AM2
      CT=0.5D0
      ENDIF
      IF(POT.EQ.'GRAVIT  ')THEN
      AMU=AM1*AM2/(AM1+AM2)
      AMU1=AMU
      AM=AMU/AM2
      CT=-AM1*AM2
      ENDIF
      TI=0.D0-DELTA
      R=RI
      THETA=0.D0
      SIG=1.D0
      RAD=E-AL**2/(2.D0*AMU*RI**2)-CT*F(RI)
      IF(RAD.LT.0.D0)GOTO 3
      DO I=NI,NF
      TI=TI+DELTA
      RAD=E-AL**2/(2.D0*AMU*R**2)-CT*F(R)
      IF(RAD.LT.0.D0.OR.R.LT.0.D0)THEN
      SIG=-1.D0*SIG
   2  R=R+SIG*0.00001D0
      RAD=E-AL**2/(2.D0*AMU*R**2)-CT*F(R)
      IF(RAD.LT.0.D0.OR.R.LT.0.D0)GOTO 2
      ENDIF
      R=R+SIG*(DSQRT(2.D0/AMU1)*DSQRT(RAD))*DELTA
      THETA=THETA+(AL/(AMU*R**2))*DELTA
      X1=AM*R*DCOS(THETA)
      Y1=AM*R*DSIN(THETA)
      X2=-(AM2/(AM1+AM2))*R*DCOS(THETA)
      Y2=-(AM2/(AM1+AM2))*R*DSIN(THETA)
      IF(POT.EQ.'MGZ     ')THEN      
      WRITE(*,1)TI,X1,Y1      
      ELSE
      WRITE(*,5)TI,X1,Y1,X2,Y2
      ENDIF      
      ENDDO
      GOTO 4
   3  WRITE(*,*)' **COMBINACION DE PARAMETROS INVALIDA**'
      WRITE(*,*)'   RADICANDO=',RAD
      WRITE(*,*)'          **PRUEBE OTRA VEZ!**'
   4  CONTINUE
   1  FORMAT(3(2X,F12.7))
   5  FORMAT(5(2X,F12.7))
      END 
C...
      FUNCTION F(R)
      IMPLICIT REAL*8 (A-H,O-Z)
C...
      COMMON/POT/POT
      CHARACTER*8 POT
C...
      IF(POT.EQ.'MGZ     ')F=R
      IF(POT.EQ.'ARMONICO')F=R**2
      IF(POT.EQ.'GRAVIT  ')F=1.D0/R
      RETURN
      END