C **********************************************************************
C "METHODS.FOR" - PROGRAM COMPUTING THEORETICAL METEOR SHOWER RADIANTS
C                 (VERSION 2)
C **********************************************************************
C
C STATUS OF USING:
C ====================
C
C   1. THE USE OF THE PROGRAM AS A SUBJECT OF COMMERCE IS NOT ALLOWED.
C   2. FOR OTHER PURPOSES, THE AUTHORS PERMIT TO COPY AND USE THE
C      PROGRAM OR ITS PART(S) FREELY;
C   3. EACH USING OF THE PROGRAM OR ANY OF ITS PARTS (SUBROUTINES)
C      SHOULD BE QUOTED IN THE OFFICIAL PRESENTATION; THE PROGRAM WAS
C      INTRODUCED IN THE "ASTRONOMY AND ASTROPHYSICS", PLEASE SEE:
C      NESLUSAN L., SVOREN J., PORUBCAN V.: 1998, "A COMPUTER PROGRAM
C      FOR CALCULATION OF A THEORETICAL METEOR-STREAM RADIANT", ASTRON.
C      ASTROPHYS. 331, 411-413.  
C
C **********************************************************************
C
C BRIEF DESCRIPTION:
C ==================
C
C      THE PROGRAM COMPUTES A THEORETICAL RADIANT OF A METEOR SHOWER,
C      IF THE ORBIT OF AN ASSUMED PARENT BODY IS WELL-KNOWN
C       
C      THE PROGRAM INCLUDES SIX DIFFERENT METHODS:
C      - (Q) Q-ADJUSTMENT (HASEGAWA; 1990)
C      - (B) VARIATION OF PERIHELION DISTANCE AND EXCENTRICITY (SVOREN
C        ET AL; 1993)
C      - (W) ROTATION OF THE LINE OF APSIDES (STEEL AND BAGGALEY; 1985)
C      - (A) ROTATION AROUND THE LINE OF APSIDES (SVOREN ET AL; 1993)
C      - (H) OMEGA-ADJUSTMENT (HASEGAWA; 1990)
C      - (P) PORTER'S METHOD (PORTER; 1952)
C     (FOR MORE DETAILED DESCRIPTION AND REFERENCES, PLEASE, SEE IN THE
C     ABOVE MENTIONED PAPERS. REFERENCES ARE ALSO GIVEN IN BRIEF INDI-
C     VIDUAL DESCRIPTIONS OF SUBROUTINES BELOW.)
C
C ---------------------------------------------------------------------
C MAIN PROGRAM
C
C   INPUT:
C   ORBITAL ELEMENTS OF THE ORBIT OF THE PARENT BODY
C      EA(1) - PERIHELION DISTANCE [AU]
C      EA(2) - ECCENTRICITY
C      EA(3) - ARGUMENT OF PERIHELION [DEGREES]
C      EA(4) - LONGITUDE OF NODE [DEGREES]
C      EA(5) - INCLINATION TO ECLIPTIC [DEGREES]
C      YEAR, MONTH, DAY - TIME OF PERIHELION PASSAGE OF THE PARENT BODY
C      (CONVERTED TO: EA(6) - THE TIME IN JULIAN CENTURIES FROM
C                             JD = 2451545.0)
C
C   OUTPUT:
C   THE QUANTITIES DESCRIBING THEORETICAL RADIANTS AT BOTH ARCS (PRE-
C   PERIHELION AND POST-PERIHELION) OF THE ORBIT OF THE PARENT BODY;
C   THE OUTPUT QUANTITIES AT THE POST-PERIHELION ARC ARE:
C      R1(1) - SOLAR LONGITUDE OF THE MAXIMUM OF POTENTIAL SHOWER 
C              [DEGRRES]
C      R1(2) - RIGHT ASCENSION OF THE PREDICTED RADIANT [DEGREES]
C      R1(3) - DECLINATION OF THE RADIANT [DEGREES]
C      R1(4) - PREDICTED GEOCENTRIC VELOCITY OF METEORS [KM/S]
C      R1(5) - PREDICTED HELIOCENTRIC VELOCITY OF METEORS [KM/S]
C      R1(6) - PREDICTED TIME OF THE SHOWER MAXIMUM [JULIAN CENTURIES
C              FROM 2000.0]
C      WM1(1) - MOMENT OF MINIMUM DISTANCE BETWEEN THE EARTH AND THE
C               ORBIT OF THE PARENT BODY [JULIAN CENTURIES FROM 2000.0]
C      WM1(2) - MINIMUM DISTANCE BETWEEN POST-PERIHELION ARC OF THE
C               PARENT BODY ORBIT AND EARTH'S ORBIT [AU]
C      WM1(3) - DISTANCE OF THE PARENT BODY FROM THE EARTH AT THE MOMENT
C               OF MAXIMUM [AU]
C      WM1(4) - TIME INTERVAL BETWEEN THE PASSAGES OF THE PARENT BODY
C               AND EARTH ACROSS THE NEAREST POINTS OF THE POST-PERI-
C               HELION ARC OF THE PARENT BODY ORBIT AND EARTH'S ORBIT
C               [DAYS]
C      WM1(5) - TRUE ANOMALY OF THE PARENT BODY AT THE MOMENT OF ITS
C               NEAREST POST-PERIHELION APPROACH TO THE EARTH'S ORBIT
C               [DEGREES]
C      WM1(6) - TRUE ANOMALY OF THE EARTH AT THE MOMENT OF ITS NEAREST
C               APPROACH TO THE POST-PERIHELION ARC OF THE PARENT BODY
C               ORBIT [DEGREES]
C   ORBITAL ELEMENTS OF MODIFIED ORBIT CROSSING THE EARTH'S ORBIT
C      EB1(1) - PERIHELION DISTANCE [AU]
C      EB1(2) - ECCENTRICITY
C      EB1(3) - ARGUMENT OF PERIHELION [DEGREES]
C      EB1(4) - LONGITUDE OF NODE [DEGREES]
C      EB1(5) - INCLINATION TO ECLIPTIC [DEGREES]
C      DD1 - D-CRITERION CHARACTERIZING THE DIFFERENCE BETWEEN THE ORBIT
C            OF THE PARENT BODY (INPUT ORBIT) AND MODIFIED ORBIT
C            CROSSING THE EARTH'S ORBIT (OUTPUT ORBIT)
C   THE ANALOGOUS QUANTITIES CONCERNING THE PRE-PERIHELION ARC ARE
C   STORED IN THE FIELDS OF VALUES "R2(I)", "WM2(I)", AND "EB2(I)",
C   WHERE I=1, 2,..., 6 (5, RESPECTIVELY), AND IN VALUE "DD2".
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*1 MET,YN
      CHARACTER*4 MON
      DIMENSION EA(6),EB1(6),EB2(6),R1(6),R2(6),RY1(6,6),RY2(6,6)
      DIMENSION EY1(6,6),EY2(6,6),DDY1(6),DDY2(6),MET(6),MON(12)
      DIMENSION WM1(6),WM2(6)
C
      MON(1)='JAN.'
      MON(2)='FEB.'
      MON(3)='MAR.'
      MON(4)='APR.'
      MON(5)='MAY '
      MON(6)='JUNE'
      MON(7)='JULY'
      MON(8)='AUG.'
      MON(9)='SEP.'
      MON(10)='OCT.'
      MON(11)='NOV.'
      MON(12)='DEC.'
C
      MET(1)='Q'
      MET(2)='B'
      MET(3)='W'
      MET(4)='A'
      MET(5)='H'
      MET(6)='P'
C
C      EA(1)=0.981656D0
C      EA(2)=0.904433D0
C      EA(3)=1.725606D2
C      EA(4)=2.351162D2
C      EA(5)=1.627073D2
C      IROK=1965
C      IMES=4
C      DEN=3.00010D1
C      GOTO 918
C
  5   CONTINUE
      WRITE(*,7)
  7   FORMAT(//)
  8   CONTINUE
      WRITE(*,10)
 10   FORMAT(' PERIHELION DISTANCE IN AU _._______')
      READ(*,20) EA(1) 
 20   FORMAT(F10.7)
      IF(EA(1).GT.1.0D-8.AND.EA(1).LT.1.0D1) GOTO 25
      WRITE(*,22)
 22   FORMAT(' ONLY A VALUE IN THE RANGE FROM 0.0000001 TO 9.9999999 IS 
     *ACCEPTABLE.')
      GOTO 8
C
 25   CONTINUE
      WRITE(*,30)
 30   FORMAT(' ECCENTRICITY _._______')
      READ(*,20) EA(2) 
      IF(EA(2).GE.0.0D0.AND.EA(2).LT.1.0D1) GOTO 35
      WRITE(*,32)
 32   FORMAT(' ONLY A VALUE IN THE RANGE FROM 0.0 TO 9.9999999 IS ACCEPT
     *ABLE.')
      GOTO 25
C
 35   CONTINUE
      WRITE(*,40)
 40   FORMAT(' ARGUMENT OF PERIHELION IN DEGREES ___._____')
      READ(*,50) EA(3) 
 50   FORMAT(F10.5)
      IF(EA(3).GE.0.0D0.AND.EA(3).LT.0.36D3) GOTO 55
      WRITE(*,52)
 52   FORMAT(' ONLY A VALUE IN THE RANGE FROM 0.0 TO 359.99999 DEGREES I 
     *S ACCEPTABLE.')
      GOTO 35
C	
 55   CONTINUE
      WRITE(*,60)
 60   FORMAT(' NODE IN DEGREES ___._____')
      READ(*,50) EA(4)
      IF(EA(4).GE.0.0D0.AND.EA(4).LT.0.36D3) GOTO 65
      WRITE(*,52)
      GOTO 55
C
 65   CONTINUE
      WRITE(*,70)
 70   FORMAT(' INCLINATION IN DEGREES ___._____')
      READ(*,50) EA(5)
      IF(EA(5).GE.0.0D0.AND.EA(5).LT.0.18D3) GOTO 75
      WRITE(*,72)
 72   FORMAT(' ONLY A VALUE IN THE RANGE FROM 0.0 TO 179.99999 DEGREES I
     *S ACCEPTABLE.')
      GOTO 65
C
 75   CONTINUE
      WRITE(*,900)
 900  FORMAT(' DATE OF PERIHELION PASSAGE - YEAR _____')
      READ(*,905) IROK
 905  FORMAT(I5)
      IF(IROK.GE.-9999.AND.IROK.LE.9999) GOTO 908
      WRITE(*,907)
 907  FORMAT(' ONLY A VALUE IN THE RANGE FROM -9999 TO 9999 IS ACCEPTABL
     *E.')
      GOTO 75
 908  CONTINUE
      WRITE(*,910)
 910  FORMAT('    - MONTH __')
      READ(*,905) IMES
      IF(IMES.GE.1.AND.IMES.LE.12) GOTO 913
      WRITE(*,912)
 912  FORMAT(' ONLY A VALUE IN THE RANGE FROM 1 TO 12 IS ACCEPTABLE.')
      GOTO 908
 913  CONTINUE
      WRITE(*,915)
 915  FORMAT('    - DAY __._____')
      READ(*,50) DEN
      IF(DEN.GE.0.0D0.AND.DEN.LT.3.2D1) GOTO 918
      WRITE(*,917)
 917  FORMAT(' ONLY A VALUE IN THE RANGE FROM 0.0 TO 31.99999 IS ACCEPTA
     *BLE.')
      GOTO 913
 918  CONTINUE
      CALL JULIAN(IROK,IMES,DEN,T19,T20)
      EA(6)=T20
 919  CONTINUE
      WRITE(*,920)
 920  FORMAT(/,' YEAR OF INVESTIGATION OF POTENTIAL SHOWER _____')
      READ(*,905) IROKI
      IF(IROK.GE.-9999.AND.IROK.LE.9999) GOTO 925
      WRITE(*,907)
      GOTO 919
C
 925  DEN=1.0D0
      IMES=1
      CALL JULIAN(IROKI,IMES,DEN,T19,TMIN)
      IROKI=IROKI+1
      CALL JULIAN(IROKI,IMES,DEN,T19,TMAX)
      IROKI=IROKI-1
C
      WRITE(*,930)
 930  FORMAT(//,' WAIT FOR A MINUTE, PLEASE, THE COMPUTATION IS EXECUTED
     *.')
      CALL APPROACH(TMIN,TMAX,EA,WM1,WM2)
C
 210  CONTINUE
      CALL QMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 215 J=1,6
      RY1(1,J)=R1(J)
      RY2(1,J)=R2(J)
      EY1(1,J)=EB1(J)
      EY2(1,J)=EB2(J)
 215  CONTINUE
      DDY1(1)=DD1
      DDY2(1)=DD2
C
 220  CONTINUE
      CALL BMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 225 J=1,6
      RY1(2,J)=R1(J)
      RY2(2,J)=R2(J)
      EY1(2,J)=EB1(J)
      EY2(2,J)=EB2(J)
 225  CONTINUE
      DDY1(2)=DD1
      DDY2(2)=DD2
C
 230  CONTINUE
      CALL WMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 235 J=1,6
      RY1(3,J)=R1(J)
      RY2(3,J)=R2(J)
      EY1(3,J)=EB1(J)
      EY2(3,J)=EB2(J)
 235  CONTINUE
      DDY1(3)=DD1
      DDY2(3)=DD2
C
 240  CONTINUE
      CALL AMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 245 J=1,6
      RY1(4,J)=R1(J)
      RY2(4,J)=R2(J)
      EY1(4,J)=EB1(J)
      EY2(4,J)=EB2(J)
 245  CONTINUE
      DDY1(4)=DD1
      DDY2(4)=DD2
C
 250  CONTINUE
      CALL HMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 255 J=1,6
      RY1(5,J)=R1(J)
      RY2(5,J)=R2(J)
      EY1(5,J)=EB1(J)
      EY2(5,J)=EB2(J)
 255  CONTINUE
      DDY1(5)=DD1
      DDY2(5)=DD2
C
 260  CONTINUE
      CALL PMETH(WM1,WM2,EA,R1,EB1,DD1,R2,EB2,DD2)
      DO 265 J=1,6
      RY1(6,J)=R1(J)
      RY2(6,J)=R2(J)
      EY1(6,J)=EB1(J)
      EY2(6,J)=EB2(J)
 265  CONTINUE
      DDY1(6)=DD1
      DDY2(6)=DD2
C
      DD1=DDY1(1)
      I1=1
      DO 270 I=2,6
      IF(DDY1(I).LT.DD1) I1=I
      IF(DDY1(I).LT.DD1) DD1=DDY1(I)
 270  CONTINUE
      DD2=DDY2(1)
      I2=1
      DO 280 I=2,6
      IF(DDY2(I).LT.DD2) I2=I
      IF(DDY2(I).LT.DD2) DD2=DDY2(I)
 280  CONTINUE
C
      IF(DD1.LT.9.9999D1.OR.DD2.LT.9.9999D1) GOTO 80
      WRITE(*,290)
 290  FORMAT(//,' NO METHOD IS APPLICABLE FOR THIS TYPE OF PARENT BODY O
     *RBIT.',/)
 292  CONTINUE
      WRITE(*,295)
 295  FORMAT(/,' OPTIONS: <D/Y> - REPEAT OF COMPUTING WITH ANOTHER INPUT
     * DATA/YEAR OF INVEST.;')
      WRITE(*,175)
      READ(*,185) YN
      IF(YN.EQ.'d'.OR.YN.EQ.'D') GOTO 190
      IF(YN.EQ.'y'.OR.YN.EQ.'Y') GOTO 919
      IF(YN.EQ.'q'.OR.YN.EQ.'Q') STOP
      GOTO 292
C
  80  CONTINUE
      WRITE(*,90)
  90  FORMAT(///,' ')
      WRITE(*,93) IROKI
  93  FORMAT('              EQUINOX: 2000.0;       DATA FOR YEAR:',I5)
      WRITE(*,110)
 110  FORMAT(' ---------------------------------------------------------
     *-----------')
      WRITE(*,100)
 100  FORMAT(' METH.  ALPHA  DELTA    VG      VH        L     DATE-MAX. 
     *   D-DISC.')
      DO 130 I=1,6
      IF(DDY2(I).GE.9.9999D1) GOTO 132
      TJD=RY2(I,6)
      CALL INVJUL(TJD,IRK,IMS,DNI)
      WRITE(*,120) MET(I),RY2(I,2),RY2(I,3),RY2(I,4),RY2(I,5),RY2(I,1),M
     *ON(IMS),DNI,DDY2(I)
      GOTO 130
 132  CONTINUE
      WRITE(*,122) MET(I)
 130  CONTINUE
 120  FORMAT(2X,'-',A1,4X,F5.1,2X,F5.1,2X,F6.2,2X,F6.2,4X,F5.1,3X,A4,1X,
     *F4.1,3X,F6.3)
 122  FORMAT(2X,'-',A1,'    * THIS METHOD IS NOT APPLICABLE IN THAT CASE
     * *')
      DO 135 I=1,6
      IF(DDY1(I).GE.9.9999D1) GOTO 137
      TJD=RY1(I,6)
      CALL INVJUL(TJD,IRK,IMS,DNI)
      WRITE(*,125) MET(I),RY1(I,2),RY1(I,3),RY1(I,4),RY1(I,5),RY1(I,1),M
     *ON(IMS),DNI,DDY1(I)
      GOTO 135
 137  CONTINUE
      WRITE(*,127) MET(I)
 135  CONTINUE
 125  FORMAT(2X,A1,'+',4X,F5.1,2X,F5.1,2X,F6.2,2X,F6.2,4X,F5.1,3X,A4,1X,
     *F4.1,3X,F6.3)
 127  FORMAT(2X,A1,'+    * THIS METHOD IS NOT APPLICABLE IN THAT CASE *'
     *)
      WRITE(*,110)
      WRITE(*,140)
 140  FORMAT(' FIRST/SECOND (-/+) SET OF DATA CONCERNS THE PRE-/POST-PER
     *IHELION ARC')
      WRITE(*,150) MET(I2),DD2
 150  FORMAT(' THE BEST METHOD - PRE-PERIHELION ARC:  ',A1,'     (D =',F
     *6.3,')')
      WRITE(*,155) MET(I1),DD1
 155  FORMAT('                 - POST-PERIHELION ARC: ',A1,'     (D =',F
     *6.3,')')
C
 180  CONTINUE
      WRITE(*,160)
 160  FORMAT(' OPTIONS: <E> - THE ELEMENTS OF FITTED ORBIT;')
      WRITE(*,165)
 165  FORMAT('          <D/Y> - REPEAT OF COMPUTING WITH ANOTHER INPUT D
     *ATA/YEAR OF INVEST.;')
      WRITE(*,170)
 170  FORMAT('          <C> - COMMENT OF DENOTATION AND UNITS OF QUANTIT
     *IES DISPLAYED;')
      WRITE(*,175)
 175  FORMAT('          <Q> - QUIT THE PROGRAM.')
      READ(*,185) YN
 185  FORMAT(A1)
      IF(YN.EQ.'e'.OR.YN.EQ.'E') GOTO 300
      IF(YN.EQ.'d'.OR.YN.EQ.'D') GOTO 190
      IF(YN.EQ.'y'.OR.YN.EQ.'Y') GOTO 919
      IF(YN.EQ.'c'.OR.YN.EQ.'C') GOTO 500
      IF(YN.EQ.'q'.OR.YN.EQ.'Q') STOP
      GOTO 180
C
 190  CONTINUE
      WRITE(*,195)
 195  FORMAT(//////////////////////////////,' ')
      GOTO 5
C
 300  CONTINUE
      WRITE(*,90)
      WRITE(*,303) IROKI
 303  FORMAT('        EQUINOX: 2000.0;        DATA FOR YEAR:',I5)
      WRITE(*,310)
 310  FORMAT(' ---------------------------------------------------------
     *-')
      WRITE(*,305)
 305  FORMAT(' METH.     Q       E      PERI.    NODE   INCL.     D-DISC
     *.')
      DO 330 I=1,6
      IF(DDY2(I).GE.9.9999D1) GOTO 332
      WRITE(*,320) MET(I),EY2(I,1),EY2(I,2),EY2(I,3),EY2(I,4),EY2(I,5),D
     *DY2(I)
      GOTO 330
 332  CONTINUE
      WRITE(*,122) MET(I)
 330  CONTINUE
 320  FORMAT(2X,'-',A1,3X,F7.4,2X,F6.4,4X,F5.1,3X,F5.1,3X,F5.1,4X,F6.3)
      DO 335 I=1,6
      IF(DDY1(I).GE.9.9999D1) GOTO 337
      WRITE(*,325) MET(I),EY1(I,1),EY1(I,2),EY1(I,3),EY1(I,4),EY1(I,5),D
     *DY1(I)
      GOTO 335
 337  CONTINUE
      WRITE(*,127) MET(I)
 335  CONTINUE
 325  FORMAT(2X,A1,'+',3X,F7.4,2X,F6.4,4X,F5.1,3X,F5.1,3X,F5.1,4X,F6.3)
      WRITE(*,310)
      WRITE(*,350)
 350  FORMAT(' ADDITIONAL INFORMATION:')
      IF(WM2(5).GT.0.36D3) GOTO 354
      WRITE(*,353) WM2(3),WM2(4),WM2(2)
 353  FORMAT('    (-) DIST. =',F8.3,' AU;   DT =',F8.1,' DAYS;   MIN. DI
     *ST. =',F7.4,' AU')
      GOTO 356
 354  WRITE(*,355)
 355  FORMAT('    (-)   * NO APPROACH OF PRE-PERIHELION ARC TO THE EARTH
     *`S ORBIT *')
 356  CONTINUE
      IF(WM1(5).LT.0.0D0) GOTO 358
      WRITE(*,357) WM1(3),WM1(4),WM1(2)
 357  FORMAT('    (+) DIST. =',F8.3,' AU;   DT =',F8.1,' DAYS;   MIN. DI
     *ST. =',F7.4,' AU')
      GOTO 380
 358  WRITE(*,359)
 359  FORMAT('    (+)   * NO APPROACH OF POST-PERIHELION ARC TO THE EART
     *H`S ORBIT *')
C
 380  CONTINUE
      WRITE(*,360)
 360  FORMAT(' OPTIONS: <G> - BACK TO THE GEOCENTRIC DATA;')
      WRITE(*,165)
      WRITE(*,170)
      WRITE(*,175)
      READ(*,185) YN
      IF(YN.EQ.'g'.OR.YN.EQ.'G') GOTO 80
      IF(YN.EQ.'d'.OR.YN.EQ.'D') GOTO 190
      IF(YN.EQ.'y'.OR.YN.EQ.'Y') GOTO 919
      IF(YN.EQ.'c'.OR.YN.EQ.'C') GOTO 500
      IF(YN.EQ.'q'.OR.YN.EQ.'Q') STOP
      GOTO 380
C
 500  CONTINUE
      WRITE(*,90)
      WRITE(*,510)
 510  FORMAT('       THE QUANTITIES CONCERNING THE POTENTIAL METEOR SHOW
     *ER DISPLAYED:')
      WRITE(*,515)
 515  FORMAT(' ---------------------------------------------------------
     *--------------------')
      WRITE(*,520)
 520  FORMAT(' ALPHA, DELTA - EQUATORIAL COORDINATES OF RADIANT [DEGREES
     *]')
      WRITE(*,530)
 530  FORMAT(' VG, VH - GEOCENTRIC AND HELIOCENTRIC VELOCITIES [KM/S]')
      WRITE(*,540)
 540  FORMAT(' L - SOLAR LONGITUDE IN THE MOMENT OF PREDICTED MAXIMUM [D
     *EGREES]')
      WRITE(*,545)
 545  FORMAT(' DATE-MAX. - DATE OF THE PREDICTED MAXIMUM [MONTH DAY]')
      WRITE(*,550)
 550  FORMAT(' THE ELEMENTS OF MEAN ORBIT: Q - PERIHELION DISTANCE [AU],
     * E - ECCENTRICITY,')
      WRITE(*,560)
 560  FORMAT('    PERI. - ARGUMENT OF PERIHELION [DEGREES], NODE - LONGI
     *TUDE OF NODE [DEGR.],')
      WRITE(*,570)
 570  FORMAT('    INCL. - INCLINATION TO ECLIPTIC [DEGREES]')
      WRITE(*,580)
 580  FORMAT(' D-DISC. - THE VALUE OF D-DISCRIMINANT; D < 0.1 (ROUGHLY) 
     *- GOOD PREDICTION,')
      WRITE(*,585)
 585  FORMAT('           D > 0.5 (ROUGHLY) - UNREAL PREDICTION')
      WRITE(*,590)
 590  FORMAT(' DIST. - DISTANCE BETWEEN THE EARTH AND PARENT BODY AT THE
     * MOMENT OF')
      WRITE(*,595)
 595  FORMAT('         THE PREDICTED MAXIMUM [AU]') 
      WRITE(*,600)
 600  FORMAT(' DT - TIME ELAPSED FROM THE PASSAGE OF THE PARENT BODY THR
     *OUGH THE POINT OF ITS')
      WRITE(*,605)
 605  FORMAT('      ORBIT NEAREST TO THE EARTH`S ORBIT AT THE INVESTIGAT
     *ED SHOWER MAXIMUM')
      WRITE(*,607)
 607  FORMAT('      (IF NEGATIVE, THE BODY WILL PASS THE POINT AFTER THE
     * MAXIMUM) [DAYS]')
      WRITE(*,610)
 610  FORMAT(' MIN. DIST. - THE MINIMUM DISTANCE BETWEEN THE GIVEN ARC O
     *F THE PARENT BODY')
      WRITE(*,615)
 615  FORMAT('              AND EARTH`S ORBITS [AU]')
C
 640  CONTINUE
      WRITE(*,620)
 620  FORMAT(/,' OPTIONS: <G> - BACK TO THE GEOCENTRIC DATA;')
      WRITE(*,630)
 630  FORMAT('          <E> - BACK TO THE ELEMENTS OF THE MEAN ORBIT;')
      WRITE(*,165)
      WRITE(*,175)
      READ(*,185) YN
      IF(YN.EQ.'g'.OR.YN.EQ.'G') GOTO 80
      IF(YN.EQ.'e'.OR.YN.EQ.'E') GOTO 300
      IF(YN.EQ.'d'.OR.YN.EQ.'D') GOTO 190
      IF(YN.EQ.'y'.OR.YN.EQ.'Y') GOTO 919
      IF(YN.EQ.'q'.OR.YN.EQ.'Q') STOP
      GOTO 640
C
      END 
C ---------------------------------------------------------------------
      SUBROUTINE QMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (Q) Q-ADJUSTMENT
C     (HASEGAWA I.: 1990, PUBL. ASTRON. SOC. JAPAN 42, 175.)
C
C     "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C     POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0. NEXT
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2,
C     3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS
C     THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D-
C     DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC
C     OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI2=2.0D0*PI
      PI180=PI/0.18D3
C
      DO 50 J=2,5
      EB1(J)=EA(J)
      EB2(J)=EA(J)
 50   CONTINUE
      DT0=(TMAX-TMIN)/3.6525D2
      DT05=DT0/2.0D0
      DD1=0.999D99
      DD2=0.999D99
C
      IF(EA(3).GT.0.18D3) GOTO 100
      FB=0.18D3-EA(3)
      IF(FB.GT.1.799999999D2.AND.EA(2).GT.0.9999999D0) GOTO 200
      VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
      GOTO 110
 100  FB=0.36D3-EA(3)
      VO0=EA(4)
 110  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 120  CONTINUE
      DO 130 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 130
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=T20
      DWM=DABS(WL-VO0)
 130  CONTINUE
      IF(DT.LT.DT05) GOTO 150
      TA=R1(6)-DT
      TB=R1(6)+DT
      DT=0.1D-1*DT
      GOTO 120
C
 150  QWV=FB*PI180
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      EB1(1)=RZM*(1.0D0+EA(2)*DCOS(QWV))/(1.0D0+EA(2))
      IF(EB1(1).LT.0.0D0) GOTO 200
      DD1=DABS(EA(1)-EB1(1))
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
C
 200  CONTINUE
      IF(EA(3).GT.0.18D3) GOTO 205
      FB=0.36D3-EA(3)
      IF(FB.LT.1.800000001D2.AND.EA(2).GT.0.9999999D0) GOTO 300
      VO0=EA(4)
      GOTO 210
 205  FB=0.54D3-EA(3)
      VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
 210  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 220  CONTINUE
      DO 230 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 230
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=T20
      DWM=DABS(WL-VO0)
 230  CONTINUE
      IF(DT.LT.DT05) GOTO 250
      TA=R2(6)-DT
      TB=R2(6)+DT
      DT=0.1D-1*DT
      GOTO 220
C
 250  QWV=FB*PI180
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      EB2(1)=RZM*(1.0D0+EA(2)*DCOS(QWV))/(1.0D0+EA(2))
      IF(EB2(1).LT.0.0D0) GOTO 300
      DD2=DABS(EA(1)-EB2(1))
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
C
 300  CONTINUE
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE BMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (B) VARIATION OF PERIHELION DISTANCE AND EXCENTRICITY
C     (SVOREN J., NESLUSAN L., PORUBCAN V.: 1993, CONTRIB. ASTRON. OBS.
C     SKALNATE PLESO 23, 23.)
C
C     "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C     POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0. NEXT
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2,
C     3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS
C     THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D-
C     DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC
C     OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI2=2.0D0*PI
      PI180=PI/0.18D3
C
      DO 50 J=3,5
      EB1(J)=EA(J)
      EB2(J)=EA(J)
 50   CONTINUE
      DT0=(TMAX-TMIN)/3.6525D2
      DT05=DT0/2.0D0
      DD1=0.999D99
      DD2=0.999D99
C
      IF(EA(3).GT.0.18D3) GOTO 100
      FB=0.18D3-EA(3)
      IF(FB.GT.1.799999999D2.AND.EA(2).GT.0.9999999D0) GOTO 200
      VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
      GOTO 110
 100  FB=0.36D3-EA(3)
      VO0=EA(4)
 110  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 120  CONTINUE
      DO 130 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 130
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=T20
      DWM=DABS(WL-VO0)
 130  CONTINUE
      IF(DT.LT.DT05) GOTO 150
      TA=R1(6)-DT
      TB=R1(6)+DT
      DT=0.1D-1*DT
      GOTO 120
C
 150  DCMIN=0.999D99
      CF=DCOS(FB*PI180)
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      DO 160 EC=0.0D0,5.0D0,0.1D0
      DE=EC-EA(2)
      QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC)
      DQ=EA(1)-QX
      DDX=DSQRT(DQ*DQ+DE*DE)
      IF(DDX.GE.DCMIN) GOTO 160
      DCMIN=DDX
      EB1(2)=EC
      EB1(1)=QX
 160  CONTINUE
C
      ED=EB1(2)-0.2D0
      IF(ED.LT.0.0D0) ED=0.0D0
      EU=EB1(2)+0.2D0
      DO 170 EC=ED,EU,0.1D-3
      DE=EC-EA(2)
      QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC)
      DQ=EA(1)-QX
      DDX=DSQRT(DQ*DQ+DE*DE)
      IF(DDX.GE.DCMIN) GOTO 170
      DCMIN=DDX
      EB1(2)=EC
      EB1(1)=QX
 170  CONTINUE
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
      DD1=DCMIN
C
 200  CONTINUE
      IF(EA(3).GT.0.18D3) GOTO 205
      FB=0.36D3-EA(3)
      IF(FB.LT.1.800000001D2.AND.EA(2).GT.0.9999999D0) GOTO 300
      VO0=EA(4)
      GOTO 210
 205  FB=0.54D3-EA(3)
      VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
 210  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 220  CONTINUE
      DO 230 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 230
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=T20
      DWM=DABS(WL-VO0)
 230  CONTINUE
      IF(DT.LT.DT05) GOTO 250
      TA=R2(6)-DT
      TB=R2(6)+DT
      DT=0.1D-1*DT
      GOTO 220
C
 250  DCMIN=0.999D99
      CF=DCOS(FB*PI180)
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      DO 260 EC=0.0D0,5.0D0,0.1D0
      DE=EC-EA(2)
      QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC)
      DQ=EA(1)-QX
      DDX=DSQRT(DQ*DQ+DE*DE)
      IF(DDX.GE.DCMIN) GOTO 260
      DCMIN=DDX
      EB2(2)=EC
      EB2(1)=QX
 260  CONTINUE
C
      ED=EB2(2)-0.2D0
      IF(ED.LT.0.0D0) ED=0.0D0
      EU=EB2(2)+0.2D0
      DO 270 EC=ED,EU,0.1D-3
      DE=EC-EA(2)
      QX=RZM*(1.0D0+EC*CF)/(1.0D0+EC)
      DQ=EA(1)-QX
      DDX=DSQRT(DQ*DQ+DE*DE)
      IF(DDX.GE.DCMIN) GOTO 270
      DCMIN=DDX
      EB2(2)=EC
      EB2(1)=QX
 270  CONTINUE
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
      DD2=DCMIN
C
 300  CONTINUE
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE WMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (W) ROTATION OF THE LINE OF APSIDES
C     (STEEL D.I., BAGGALEY J.W.: 1985, MON. NOT. R. ASTRON. SOC. 212,
C     817.)
C
C     "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C     POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0. NEXT
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2,
C     3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS
C     THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D-
C     DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC
C     OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI180=PI/0.18D3
      PI05=PI/2.0D0
      PI2=2.0D0*PI
C
      DO 10 J=1,5
      EB1(J)=EA(J)
      EB2(J)=EA(J)
 10   CONTINUE      
      IF(EA(2).GT.0.999D0) GOTO 40
      IF(EA(1).GT.1.1D0) GOTO 20
      IF(EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)).LT.0.9D0) GOTO 20
      IF(EA(2).LT.1.0D-12) GOTO 20
      GOTO 40
 20   DD1=0.999D99
      DD2=0.999D99
      RETURN
C
 40   DT0=(TMAX-TMIN)/3.6525D2
      DT05=DT0/2.0D0
C
      IF(EA(3).GT.0.18D3) GOTO 100
      VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
      GOTO 110
 100  VO0=EA(4)
 110  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 120  CONTINUE
      DO 130 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 130
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=T20
      DWM=DABS(WL-VO0)
 130  CONTINUE
      IF(DT.LT.DT05) GOTO 150
      TA=R1(6)-DT
      TB=R1(6)+DT
      DT=0.1D-1*DT
      GOTO 120
C
 150  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      CFB=(EA(1)*(1.0D0+EA(2))/RZM-1.0D0)/EA(2)
      DD1=0.999D99
      IF(DABS(CFB).GT.1.000001D0) GOTO 200
      IF(CFB.LT.-1.0D0) CFB=-1.0D0
      IF(CFB.GT.1.0D0) CFB=1.0D0
      IF(DABS(CFB).LT.1.0D-12) GOTO 160
      SFB=DSQRT(1.0D0-CFB*CFB)
      FB=DATAN(SFB/CFB)
      IF(FB.LT.0.0D0) FB=FB+PI
      GOTO 170
 160  FB=PI05
 170  FB=FB/PI180
      EB1(3)=0.18D3-FB
      IF(EA(3).GT.0.18D3) EB1(3)=0.36D3-FB
      CALL DCRIT(EA,EB1,DD1)
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
C
 200  CONTINUE
      IF(EA(3).GT.0.18D3) GOTO 205
      VO0=EA(4)
      GOTO 210
 205  VO0=EA(4)+0.18D3
      IF(VO0.GE.0.36D3) VO0=VO0-0.36D3
 210  VO0=VO0*PI180
      DWM=2.0D0*PI2
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 220  CONTINUE
      DO 230 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 230
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=T20
      DWM=DABS(WL-VO0)
 230  CONTINUE
      IF(DT.LT.DT05) GOTO 250
      TA=R2(6)-DT
      TB=R2(6)+DT
      DT=0.1D-1*DT
      GOTO 220
C
 250  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZM=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      CFB=(EA(1)*(1.0D0+EA(2))/RZM-1.0D0)/EA(2)
      DD2=0.999D99
      IF(DABS(CFB).GT.1.000001D0) GOTO 300
      IF(CFB.LT.-1.0D0) CFB=-1.0D0
      IF(CFB.GT.1.0D0) CFB=1.0D0
      IF(DABS(CFB).LT.1.0D-12) GOTO 260
      SFB=DSQRT(1.0D0-CFB*CFB)
      FB=DATAN(SFB/CFB)
      IF(FB.LT.0.0D0) FB=FB+PI
      GOTO 270
 260  FB=PI05
 270  FB=0.36D3-FB/PI180
      EB2(3)=0.36D3-FB
      IF(EA(3).GT.0.18D3) EB2(3)=0.54D3-FB
      CALL DCRIT(EA,EB2,DD2)
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
C
 300  CONTINUE
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE AMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (A) ROTATION AROUND THE LINE OF APSIDES
C     (SVOREN J., NESLUSAN L., PORUBCAN V.: 1993, CONTRIB. ASTRON. OBS.
C     SKALNATE PLESO 23, 23.)
C
C     "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C     POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0. NEXT
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2,
C     3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS
C     THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D-
C     DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC
C     OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI05=PI/2.0D0
      PI2=2.0D0*PI
      PI180=PI/0.18D3
C
      DD1=0.999D99
      DD2=0.999D99
      IF(EA(2).LT.1.0D-12) GOTO 5
      IF(EA(2).GT.0.999D0) GOTO 10
      IF(EA(1).LT.1.1D0) GOTO 10
      IF(EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2)).GT.0.9D0) GOTO 10
   5  CONTINUE
      RETURN
C
  10  EB1(1)=EA(1)
      EB1(2)=EA(2)
      EB2(1)=EA(1)
      EB2(2)=EA(2)
      CSO=DCOS(PI180*EA(3))
      SSO=DSIN(PI180*EA(3))
      CHO=DCOS(PI180*EA(4))
      SHO=DSIN(PI180*EA(4))
      CSK=DCOS(PI180*EA(5))
      SSK=DSIN(PI180*EA(5))
      XQ=CSO*CHO-SSO*CSK*SHO
      YQ=CSO*SHO+SSO*CSK*CHO
      ZQ=SSO*SSK
      RQ=DSQRT(XQ*XQ+YQ*YQ+ZQ*ZQ)
      IF(DABS(RQ-1.0D0).GT.1.0D-12) WRITE(*,11)
 11   FORMAT(/,'  ERROR IN SUBROUTINE `AMETH` (RQ.NE.1).')
      IF(DABS(RQ-1.0D0).GT.1.0D-12) STOP
C
      DT0=(TMAX-TMIN)/3.6525D2
      DT05=DT0/2.0D0
C
      DO 13 T20=TMIN,TMAX,DT05
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      RZ0=WEL(1)*(1.0D0-WEL(2)*WEL(2))/(1.0D0+WEL(2)*DCOS(WEL(7)))
      CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2)
      IF(DABS(CF0).LE.1.0D0) GOTO 15
 13   CONTINUE
      GOTO 5
C
 15   CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 17
      IF(DABS(CF0).GT.1.000000001D0) GOTO 200
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 20
  17  F0=PI05
  20  EB1(3)=PI-F0
      IF(EA(3).GT.0.18D3) EB1(3)=PI2-F0
      CSOB=DCOS(EB1(3))
      SSOB=DSIN(EB1(3))
      EB1(3)=EB1(3)/PI180
      QKOB=SSOB*SSOB+CSOB*CSOB
      IF(DABS(QKOB-1.0D0).GT.1.0D-12) WRITE(*,30)
  30  FORMAT(/,'  ERROR IN SUBROUTINE `AMETH` (QKOB.NE.1).')
      IF(DABS(QKOB-1.0D0).GT.1.0D-12) STOP
      IF(DABS(SSOB).LT.1.0D-9.AND.DABS(SSO*SSK).GT.1.0D-9) GOTO 200
      IF(DABS(SSOB).LT.1.0D-9) GOTO 40
      SIB=SSO*SSK/SSOB
      IF(DABS(SIB).GT.1.0000000001D0) GOTO 200
      GOTO 60
  40  SIB=SSK
  60  CONTINUE
      IF(SIB.GT.1.0D0) SIB=1.0D0
      IF(SIB.LT.-1.0D0) SIB=-1.0D0
      CIB=DSQRT(1.0D0-SIB*SIB)
      JIB=1
  65  CONTINUE
      IF(DABS(CIB).LT.1.0D-12) EB1(5)=PI05
      IF(DABS(CIB).LT.1.0D-12) GOTO 70
      EB1(5)=DATAN(SIB/CIB)
      IF(CIB.LT.0.0D0) EB1(5)=PI+EB1(5)
  70  CONTINUE
      EB1(5)=EB1(5)/PI180
      WEVO=CSOB*CSOB+SSOB*SSOB*CIB*CIB
C      IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) DD1=0.999D99
      IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) GOTO 200
      QQS=(YQ*CSOB-XQ*SSOB*CIB)/WEVO
      QQC=(XQ*CSOB+YQ*SSOB*CIB)/WEVO
      QKONTR=QQS*QQS+QQC*QQC
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,80)
  80  FORMAT(/,'  ERROR IN SUBROUTINE `AMETH` ((SIN**2(EB(4)) + COS**(EB
     *(4))).GT.1).')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP     
      IF(DABS(QQC).LT.1.0D-12) EB1(4)=PI05
      IF(DABS(QQC).LT.1.0D-12) GOTO 90
      EB1(4)=DATAN(QQS/QQC)
      EB1(4)=DABS(EB1(4))
      IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) EB1(4)=PI2-EB1(4)
      IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) EB1(4)=PI-EB1(4)
      IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) EB1(4)=PI+EB1(4)
  90  CONTINUE
      EB1(4)=EB1(4)/PI180
      CALL DCRIT(EA,EB1,DD1)
      IF(JIB.EQ.2) GOTO 110
      SOBO=EB1(3)
      HOBO=EB1(4)
      SKBO=EB1(5)
      DDO=DD1
 100  JIB=2
      CIB=-CIB
      GOTO 65
 110  CONTINUE
      TA=TMIN
      TB=TMAX
      DT=DT0
      IF(DD1.LT.DDO) GOTO 120
      EB1(3)=SOBO
      EB1(4)=HOBO
      EB1(5)=SKBO
      DD1=DDO
C
 120  DWM=2.0D0*PI2
      VO0=EB1(4)*PI180+PI
      IF(EB1(3).GT.0.18D3) VO0=EB1(4)*PI180
      IF(VO0.GE.PI2) VO0=VO0-PI2
      DO 130 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 130
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=T20
      DWM=DABS(WL-VO0)
 130  CONTINUE
      IF(DT.LT.DT05) GOTO 150
      TA=R1(6)-DT
      TB=R1(6)+DT
      DT=0.1D-1*DT
      GOTO 120
C
 150  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))      
      IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 160
      RZ0=RZ
      CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2)
      IF(DABS(CF0).LE.1.0D0) GOTO 15
      GOTO 200
 160  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 170
      IF(DABS(CF0).GT.1.000000001D0) GOTO 200
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 180
 170  F0=PI05
 180  EB1(3)=(PI-F0)/PI180
      IF(EA(3).GT.0.18D3) EB1(3)=(PI2-F0)/PI180
      FB=F0/PI180
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
C
C   2-ND ARC
C
 200  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 217
      IF(DABS(CF0).GT.1.000000001D0) GOTO 400
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 220
 217  F0=PI05
 220  F0=PI2-F0
      EB2(3)=PI2-F0
      IF(EA(3).GT.0.18D3) EB2(3)=PI2+PI-F0
      CSOB=DCOS(EB2(3))
      SSOB=DSIN(EB2(3))
      EB2(3)=EB2(3)/PI180
      QKOB=SSOB*SSOB+CSOB*CSOB
      IF(DABS(QKOB-1.0D0).GT.1.0D-12) WRITE(*,230)
 230  FORMAT(/,'  ERROR IN SUBROUTINE `AMETH` (QKOB.NE.1).')
      IF(DABS(QKOB-1.0D0).GT.1.0D-12) STOP
      IF(DABS(SSOB).LT.1.0D-12.AND.DABS(SSO).GT.1.0D-9) GOTO 400
      IF(DABS(SSOB).LT.1.0D-12) GOTO 245
      SIB=SSO*SSK/SSOB
      IF(DABS(SIB).GT.1.0000000001D0) GOTO 400
      GOTO 260
 245  SIB=SSK
 260  CONTINUE
      IF(SIB.GT.1.0D0) SIB=1.0D0
      IF(SIB.LT.-1.0D0) SIB=-1.0D0
      CIB=DSQRT(1.0D0-SIB*SIB)
      JIB=1
 265  CONTINUE
      IF(DABS(CIB).LT.1.0D-12) EB2(5)=PI05
      IF(DABS(CIB).LT.1.0D-12) GOTO 270
      EB2(5)=DATAN(SIB/CIB)
      IF(CIB.LT.0.0D0) EB2(5)=PI+EB2(5)
 270  CONTINUE
      EB2(5)=EB2(5)/PI180
      WEVO=CSOB*CSOB+SSOB*SSOB*CIB*CIB
C      IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) DD2=0.999D99
      IF(DABS(CSOB).LT.1.0D-12.AND.DABS(CIB).LT.1.0D-12) GOTO 400
      QQS=(YQ*CSOB-XQ*SSOB*CIB)/WEVO
      QQC=(XQ*CSOB+YQ*SSOB*CIB)/WEVO
      QKONTR=QQS*QQS+QQC*QQC
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,280)
 280  FORMAT(/,'  ERROR IN SUBROUTINE `AMETH` ((SIN**2(EB(4)) + COS**(EB
     *(4))).GT.1).')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP     
      IF(DABS(QQC).LT.1.0D-12) EB2(4)=PI05
      IF(DABS(QQC).LT.1.0D-12) GOTO 290
      EB2(4)=DATAN(QQS/QQC)
      EB2(4)=DABS(EB2(4))
      IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) EB2(4)=PI2-EB2(4)
      IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) EB2(4)=PI-EB2(4)
      IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) EB2(4)=PI+EB2(4)
 290  CONTINUE
      EB2(4)=EB2(4)/PI180
      CALL DCRIT(EA,EB2,DD2)
      IF(JIB.EQ.2) GOTO 310
      SOBO=EB2(3)
      HOBO=EB2(4)
      SKBO=EB2(5)
      DDO=DD2
 300  JIB=2
      CIB=-CIB
      GOTO 265
 310  CONTINUE
      TA=TMIN
      TB=TMAX
      DT=DT0
      IF(DD2.LT.DDO) GOTO 320
      EB2(3)=SOBO
      EB2(4)=HOBO
      EB2(5)=SKBO
      DD2=DDO
C
 320  DWM=2.0D0*PI2
      VO0=EB2(4)*PI180
      IF(EB2(3).GT.0.18D3) VO0=EB2(4)*PI180+PI
      IF(VO0.GE.PI2) VO0=VO0-PI2
      DO 330 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-VO0).GT.DWM) GOTO 330
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=T20
      DWM=DABS(WL-VO0)
 330  CONTINUE
      IF(DT.LT.DT05) GOTO 350
      TA=R2(6)-DT
      TB=R2(6)+DT
      DT=0.1D-1*DT
      GOTO 320
C
 350  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))      
      IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 360
      RZ0=RZ
      CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2)
      IF(DABS(CF0).LE.1.0D0) GOTO 200
      GOTO 400
 360  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 370
      IF(DABS(CF0).GT.1.000000001D0) GOTO 400
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 380
 370  F0=PI05
 380  F0=PI2-F0
      EB2(3)=(PI2-F0)/PI180
      IF(EA(3).GT.0.18D3) EB2(3)=(PI2+PI-F0)/PI180
      FB=F0/PI180
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
C
 400  CONTINUE
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE HMETH(TMIN,TMAX,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (H) OMEGA-ADJUSTMENT
C     (HASEGAWA I.: 1990, PUBL. ASTRON. SOC. JAPAN 42, 175.)
C
C     "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C     POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0. NEXT
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "EB1(I)", AND "EB2(I)", I = 1, 2,
C     3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE MAIN PROGRAM. AS
C     THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1" AND "DD2" OF D-
C     DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-PERIHELION ARC
C     OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI180=PI/0.18D3
      PI05=PI/2.0D0
      PI2=2.0D0*PI
C
      GOTO 10
   5  DD1=0.999D99
      DD2=0.999D99
      RETURN
  10  EB1(1)=EA(1)
      EB1(2)=EA(2)
      EB2(1)=EA(1)
      EB2(2)=EA(2)
      CSK=DCOS(PI180*EA(5))
      SSK=DSIN(PI180*EA(5))
C
      DT0=(TMAX-TMIN)/3.6525D2
      DT05=DT0/2.0D0
C
      DO 13 T20=TMIN,TMAX,DT05
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      RZ0=WEL(1)*(1.0D0-WEL(2)*WEL(2))/(1.0D0+WEL(2)*DCOS(WEL(7)))
      CF0=(EA(1)*(1.0D0+EA(2))/RZ0-1.0D0)/EA(2)
      IF(DABS(CF0).LE.1.0D0) GOTO 15
 13   CONTINUE
      GOTO 5
C
 15   CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 17
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 20
  17  F0=PI05
  20  U=EA(3)*PI180+F0
      IF(U.GE.PI2) U=U-PI2
      CU=DCOS(U)
      SU=DSIN(U)
      IF(DABS(CU).LT.1.0D-12) HLVO=PI05
      IF(DABS(CU).LT.1.0D-12) GOTO 30
      HLVO=DATAN(SU*CSK/CU)
      HLVO=DABS(HLVO)
  30  CONTINUE
      IF(SU*CSK.GE.0.0D0.AND.CU.LT.0.0D0) HLVO=PI-HLVO
      IF(SU*CSK.LT.0.0D0.AND.CU.GE.0.0D0) HLVO=PI2-HLVO
      IF(SU*CSK.LT.0.0D0.AND.CU.LT.0.0D0) HLVO=PI+HLVO
      HL=HLVO+EA(4)*PI180
      IF(HL.GE.PI2) HL=HL-PI2
C
      TA=TMIN
      TB=TMAX
      DT=DT0
      DWM=2.0D0*PI2
C
  40  CONTINUE
      DO 50 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-HL).GT.DWM) GOTO 50
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=T20
      DWM=DABS(WL-HL)
  50  CONTINUE
      IF(DT.LT.DT05) GOTO 60
      TA=R1(6)-DT
      TB=R1(6)+DT
      DT=0.1D-1*DT
      GOTO 40
C
  60  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      CF0=(EA(1)*(1.0D0+EA(2))/RZ-1.0D0)/EA(2)
      IF(DABS(CF0).GT.1.000000001D0) GOTO 5
      IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 65
      RZ0=RZ
      GOTO 15
  65  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 67
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      GOTO 70
  67  F0=PI05
C
  70  CONTINUE
      EB1(3)=PI2-F0
      EB1(4)=HL
      U=U/PI180
      IF(U.LE.9.0D1.OR.U.GT.2.7D2) GOTO 80
      EB1(3)=EB1(3)-PI
      IF(EB1(3).LT.0.0D0) EB1(3)=EB1(3)+PI2
      EB1(4)=EB1(4)+PI
      IF(EB1(4).GE.PI2) EB1(4)=EB1(4)-PI2
  80  EB1(3)=EB1(3)/PI180
      EB1(4)=EB1(4)/PI180
C
      SIB=SSK*DCOS(HLVO)
      IF(DABS(SIB).GT.1.0000000001D0) GOTO 90
      GOTO 100
  90  CONTINUE
      WRITE(*,95)
  95  FORMAT(/,'  ERROR IN SUBROUTINE `HMETH` (ABS(SIN(EB(5))).GT.1).')
      STOP
 100  CONTINUE
      IF(SIB.GT.1.0D0) SIB=1.0D0
      IF(SIB.LT.-1.0D0) SIB=-1.0D0
      CIB=DSQRT(1.0D0-SIB*SIB)
      IF(DABS(CIB).LT.1.0D-12) EB1(5)=PI05
      IF(DABS(CIB).LT.1.0D-12) GOTO 110
      EB1(5)=DATAN(SIB/CIB)
      IF(EB1(5).LT.0.0D0) EB1(5)=PI+EB1(5)
 110  EB1(5)=EB1(5)/PI180
      IF(EA(5).GT.9.0D1.AND.EB1(5).LT.9.0D1) EB1(5)=1.8D2-EB1(5)
      IF(EA(5).LT.9.0D1.AND.EB1(5).GT.9.0D1) EB1(5)=1.8D2-EB1(5)
C
      CALL DCRIT(EA,EB1,DD1)
      FB=F0/PI180
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
C
C   2-ND ARC
C
 200  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 217
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      F0=PI2-F0
      GOTO 220
 217  F0=PI2-PI05
 220  U=EA(3)*PI180+F0
      IF(U.GE.PI2) U=U-PI2
      CU=DCOS(U)
      SU=DSIN(U)
      IF(DABS(CU).LT.1.0D-12) HLVO=PI05
      IF(DABS(CU).LT.1.0D-12) GOTO 230
      HLVO=DATAN(SU*CSK/CU)
      HLVO=DABS(HLVO)
 230  CONTINUE
      IF(SU*CSK.GE.0.0D0.AND.CU.LT.0.0D0) HLVO=PI-HLVO
      IF(SU*CSK.LT.0.0D0.AND.CU.GE.0.0D0) HLVO=PI2-HLVO
      IF(SU*CSK.LT.0.0D0.AND.CU.LT.0.0D0) HLVO=PI+HLVO
      HL=HLVO+EA(4)*PI180
      IF(HL.GE.PI2) HL=HL-PI2
C
      TA=TMIN
      TB=TMAX
      DT=DT0
      DWM=2.0D0*PI2
C
 240  CONTINUE
      DO 250 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      IF(DABS(WL-HL).GT.DWM) GOTO 250
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=T20
      DWM=DABS(WL-HL)
 250  CONTINUE
      IF(DT.LT.DT05) GOTO 260
      TA=R2(6)-DT
      TB=R2(6)+DT
      DT=0.1D-1*DT
      GOTO 240
C
 260  CONTINUE
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZM))
      CF0=(EA(1)*(1.0D0+EA(2))/RZ-1.0D0)/EA(2)
      IF(DABS(CF0).GT.1.000000001D0) GOTO 5
      IF(DABS(RZ-RZ0).LT.1.0D-4) GOTO 265
      RZ0=RZ
      GOTO 200
 265  CONTINUE
      IF(DABS(CF0).LT.1.0D-9) GOTO 267
      IF(CF0.GT.1.0D0) CF0=1.0D0
      IF(CF0.LT.-1.0D0) CF0=-1.0D0
      SF0=DSQRT(1.0D0-CF0*CF0)
      F0=DATAN(SF0/CF0)
      IF(F0.LT.0.0D0) F0=F0+PI
      F0=PI2-F0
      GOTO 270
 267  F0=PI2-PI05
C
 270  CONTINUE
      EB2(3)=PI2-F0
      EB2(4)=HL
      U=U/PI180
      IF(U.LE.9.0D1.OR.U.GT.2.7D2) GOTO 280
      EB2(3)=EB2(3)-PI
      IF(EB2(3).LT.0.0D0) EB2(3)=EB2(3)+PI2
      EB2(4)=EB2(4)+PI
      IF(EB2(4).GE.PI2) EB2(4)=EB2(4)-PI2
 280  EB2(3)=EB2(3)/PI180
      EB2(4)=EB2(4)/PI180
C
      SIB=SSK*DCOS(HLVO)
      IF(DABS(SIB).GT.1.0000000001D0) GOTO 290
      GOTO 300
 290  CONTINUE
      WRITE(*,95)
      STOP
 300  CONTINUE
      IF(SIB.GT.1.0D0) SIB=1.0D0
      IF(SIB.LT.-1.0D0) SIB=-1.0D0
      CIB=DSQRT(1.0D0-SIB*SIB)
      IF(DABS(CIB).LT.1.0D-12) EB2(5)=PI05
      IF(DABS(CIB).LT.1.0D-12) GOTO 310
      EB2(5)=DATAN(SIB/CIB)
      IF(EB2(5).LT.0.0D0) EB2(5)=PI+EB2(5)
 310  EB2(5)=EB2(5)/PI180
      IF(EA(5).GT.9.0D1.AND.EB2(5).LT.9.0D1) EB2(5)=1.8D2-EB2(5)
      IF(EA(5).LT.9.0D1.AND.EB2(5).GT.9.0D1) EB2(5)=1.8D2-EB2(5)
C
      CALL DCRIT(EA,EB2,DD2)
      FB=F0/PI180
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
C
 400  CONTINUE
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE PMETH(WM1,WM2,EA,R1,EB1,DD1,R2,EB2,DD2)
C ---------------------------------------------------------------------
C      - (P) PORTER'S METHOD
C     (PORTER J.G.: 1952, COMETS AND METEOR STREAMS, CHAPMAN AND HALL
C     LTD., LONDON.)
C
C     INPUT (FIELD OF VALUES "EA(I)", I = 1, 2,..., 6) AND OUTPUT (FIEL-
C     DS OF VALUES "R1(I)", "R2(I)", "WM1(I)", WM2(I)", "EB1(I)", AND
C     "EB2(I)", I = 1, 2, 3,..., 6 OR 5, RESPECTIVELY), PLEASE, SEE THE
C     MAIN PROGRAM. AS THE OUTPUT, THERE ARE MOREOVER THE VALUES "DD1"
C     AND "DD2" OF D-DISCRIMINANT CORRESPONDING TO THE POST- AND PRE-
C     PERIHELION ARC OF THE PARENT BODY ORBIT, RESPECTIVELY.
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB1(5),EB2(5),R1(6),R2(6),WM1(6),WM2(6)
      DIMENSION WEL(7),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI05=PI/2.0D0
      PI2=2.0D0*PI
      PI180=PI/0.18D3
      GAUSS=.17202098950D-1
      GS2=GAUSS*GAUSS
C
      CHO=DCOS(EA(4)*PI180)
      SHO=DSIN(EA(4)*PI180)
      CSK=DCOS(EA(5)*PI180)
      SSK=DSIN(EA(5)*PI180)
C
      DD1=0.999D99
      IF(WM1(2).GT.1.0D2.OR.WM1(5).LT.0.0D0) GOTO 200
      CVK=DCOS(WM1(5)*PI180)
      SVK=DSIN(WM1(5)*PI180)
      WVFI=DSQRT(1.0D0+EA(2)*EA(2)+2.0D0*EA(2)*CVK)
      CHFI=-EA(2)*SVK/WVFI
      SHFI=(1.0D0+EA(2)*CVK)/WVFI
      QKONTR=SHFI*SHFI+CHFI*CHFI
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,5)
   5  FORMAT(/,'  ERROR IN SUBROUTINE `PMETH` ((SIN**2(HFI) + COS**2(VFI
     *).NE.1).')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP     
      IF(SHFI.LT.-1.0D-12) WRITE(*,7)
   7  FORMAT(/,'  ERROR IN SUBROUTINE `PMETH` (ANGLE `HFI` > PI).')
      IF(SHFI.LT.-1.0D-12) STOP     
      IF(DABS(CHFI).LT.1.0D-12) HFI=PI05
      IF(DABS(CHFI).LT.1.0D-12) GOTO 10
      HFI=DATAN(SHFI/CHFI)
      IF(HFI.LT.0.0D0) HFI=PI+HFI
  10  TMN1=WM1(1)
      CALL ECLXYZ(TMN1,X,Y,Z,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R1(1)=WL/PI180
      R1(6)=TMN1
      RZ=DSQRT(X*X+Y*Y+Z*Z)
      IF(EA(1).GT.RZ) GOTO 20
      IF(EA(2).GE.0.9999999D0) GOTO 30
      QQK=EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2))
      IF(QQK.GE.RZ) GOTO 30
      AV=GAUSS*DSQRT(2.0D0/QQK-(1.0D0-EA(2))/EA(1))
      GOTO 40
  20  CONTINUE
      AV=GAUSS*DSQRT((1.0D0+EA(2))/EA(1))
      GOTO 40
  30  CONTINUE
      AV=GAUSS*DSQRT(2.0D0/RZ-(1.0D0-EA(2))/EA(1))
  40  AR3S=(WM1(5)+EA(3))*PI180-HFI
      C3S=DCOS(AR3S)
      S3S=DSIN(AR3S)
      VX=-AV*(C3S*CHO-S3S*CSK*SHO)
      VY=-AV*(C3S*SHO+S3S*CSK*CHO)
      VZ=-AV*S3S*SSK
C
      WMEN=(X*VY-Y*VX)*(X*VY-Y*VX)+(Y*VZ-Z*VY)*(Y*VZ-Z*VY)
      WMEN=DSQRT(WMEN+(Z*VX-X*VZ)*(Z*VX-X*VZ))
      CIB=(X*VY-Y*VX)/WMEN
      SIB=DSQRT(1.0D0-CIB*CIB)
      IF(DABS(CIB).LT.1.0D-12) GOTO 50
      EB1(5)=DATAN(SIB/CIB)
      IF(EB1(5).LT.0.0D0) EB1(5)=PI+EB1(5)
      GOTO 60
  50  EB1(5)=PI05
  60  WMEN=(Y*VZ-Z*VY)*(Y*VZ-Z*VY)+(Z*VX-X*VZ)*(Z*VX-X*VZ)
      WMEN=DSQRT(WMEN)
      CVO=(X*VZ-Z*VX)/WMEN
      SVO=(Y*VZ-Z*VY)/WMEN
      IF(DABS(CVO).LT.1.0D-12) GOTO 70
      EB1(4)=DATAN(SVO/CVO)
      EB1(4)=DABS(EB1(4))
      GOTO 80
  70  EB1(4)=PI05
  80  CONTINUE
      IF(CVO.GE.0.0D0.AND.SVO.LT.0.0D0) EB1(4)=PI2-EB1(4)
      IF(CVO.LT.0.0D0.AND.SVO.GE.0.0D0) EB1(4)=PI-EB1(4)
      IF(CVO.LT.0.0D0.AND.SVO.LT.0.0D0) EB1(4)=PI+EB1(4)
C
      AKI=2.0D0/RZ-AV*AV/GS2
      EB1(2)=DSQRT(1.0D0-RZ*AKI*(2.0D0-RZ*AKI)*SHFI*SHFI)
      IF(DABS(EB1(2)-1.0D0).LT.1.0D-9) GOTO 90
      EB1(1)=(1.0D0-EB1(2))/AKI
      GOTO 100
  90  EB1(1)=RZ*SHFI*SHFI
 100  CF=(EB1(1)*(1.0D0+EB1(2))/RZ-1.0D0)/EB1(2)
      IF(CF.GE.0.999999999D0) GOTO 115
      IF(CF.LE.-0.999999999D0) GOTO 117
      SF=DSQRT(1.0D0-CF*CF)
      IF(DABS(CF).LT.1.0D-12) GOTO 110
      F=DATAN(SF/CF)
      IF(F.LT.0.0D0) F=PI+F
      GOTO 120
 110  F=PI05
      GOTO 120
 115  F=0.0D0
      GOTO 120
 117  F=PI
 120  CONTINUE
      IF(HFI.LT.PI05) F=PI2-F
C
      C2S=(X*CVO+Y*SVO)/RZ
      IF(DABS(CIB).LT.1.0D-9) GOTO 130
      S2S=(Y*CVO-X*SVO)/RZ/CIB
      GOTO 140
 130  S2S=Z/RZ/SIB
 140  CONTINUE
      IF(DABS(C2S).LT.1.0D-12) GOTO 150
      U2S=DATAN(S2S/C2S)
      U2S=DABS(U2S)
      GOTO 160
 150  U2S=PI05
 160  CONTINUE
      IF(C2S.GE.0.0D0.AND.S2S.LT.0.0D0) U2S=PI2-U2S
      IF(C2S.LT.0.0D0.AND.S2S.GE.0.0D0) U2S=PI-U2S
      IF(C2S.LT.0.0D0.AND.S2S.LT.0.0D0) U2S=PI+U2S
      EB1(3)=U2S-F
      IF(EB1(3).LT.0.0D0) EB1(3)=EB1(3)+PI2
C
      EB1(3)=EB1(3)/PI180
      EB1(4)=EB1(4)/PI180
      EB1(5)=EB1(5)/PI180
      CALL DCRIT(EA,EB1,DD1)
      FB=F/PI180
      R1(1)=R1(1)-0.18D3
      IF(R1(1).LT.0.0D0) R1(1)=R1(1)+0.36D3
      CALL RADIANT(EB1,FB,EZ,FZM,R1)
C
C   2-ND ARC
C
 200  CONTINUE
      DD2=0.999D99
      IF(WM2(2).GT.1.0D2.OR.WM2(5).GT.0.36D3) RETURN
      CVK=DCOS(WM2(5)*PI180)
      SVK=DSIN(WM2(5)*PI180)
      WVFI=DSQRT(1.0D0+EA(2)*EA(2)+2.0D0*EA(2)*CVK)
      CHFI=-EA(2)*SVK/WVFI
      SHFI=(1.0D0+EA(2)*CVK)/WVFI
      QKONTR=SHFI*SHFI+CHFI*CHFI
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,5)
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP     
      IF(SHFI.LT.-1.0D-12) WRITE(*,7)
      IF(SHFI.LT.-1.0D-12) STOP     
      IF(DABS(CHFI).LT.1.0D-12) HFI=PI05
      IF(DABS(CHFI).LT.1.0D-12) GOTO 210
      HFI=DATAN(SHFI/CHFI)
      IF(HFI.LT.0.0D0) HFI=PI+HFI
 210  TMN2=WM2(1)
      CALL ECLXYZ(TMN2,X,Y,Z,WEL)
      WL=WEL(7)+WEL(4)
      IF(WL.GE.PI2) WL=WL-PI2
      EZ(1)=WEL(1)
      EZ(2)=WEL(2)
      EZ(3)=WEL(4)-WEL(5)
      IF(EZ(3).LT.0.0D0) EZ(3)=EZ(3)+PI2
      EZ(4)=WEL(5)
      EZ(5)=WEL(6)
      FZM=WEL(7)
      R2(1)=WL/PI180
      R2(6)=TMN2
      RZ=DSQRT(X*X+Y*Y+Z*Z)
      IF(EA(1).GT.RZ) GOTO 220
      IF(EA(2).GE.0.9999999D0) GOTO 230
      QQK=EA(1)*(1.0D0+EA(2))/(1.0D0-EA(2))
      IF(QQK.GE.RZ) GOTO 230
      AV=GAUSS*DSQRT(2.0D0/QQK-(1.0D0-EA(2))/EA(1))
      GOTO 240
 220  CONTINUE
      AV=GAUSS*DSQRT((1.0D0+EA(2))/EA(1))
      GOTO 240
 230  CONTINUE
      AV=GAUSS*DSQRT(2.0D0/RZ-(1.0D0-EA(2))/EA(1))
 240  AR3S=(WM2(5)+EA(3))*PI180-HFI
      C3S=DCOS(AR3S)
      S3S=DSIN(AR3S)
      VX=-AV*(C3S*CHO-S3S*CSK*SHO)
      VY=-AV*(C3S*SHO+S3S*CSK*CHO)
      VZ=-AV*S3S*SSK
C
      WMEN=(X*VY-Y*VX)*(X*VY-Y*VX)+(Y*VZ-Z*VY)*(Y*VZ-Z*VY)
      WMEN=DSQRT(WMEN+(Z*VX-X*VZ)*(Z*VX-X*VZ))
      CIB=(X*VY-Y*VX)/WMEN
      SIB=DSQRT(1.0D0-CIB*CIB)
      IF(DABS(CIB).LT.1.0D-12) GOTO 250
      EB2(5)=DATAN(SIB/CIB)
      IF(EB2(5).LT.0.0D0) EB2(5)=PI+EB2(5)
      GOTO 260
 250  EB2(5)=PI05
 260  WMEN=(Y*VZ-Z*VY)*(Y*VZ-Z*VY)+(Z*VX-X*VZ)*(Z*VX-X*VZ)
      WMEN=DSQRT(WMEN)
      CVO=(X*VZ-Z*VX)/WMEN
      SVO=(Y*VZ-Z*VY)/WMEN
      IF(DABS(CVO).LT.1.0D-12) GOTO 270
      EB2(4)=DATAN(SVO/CVO)
      EB2(4)=DABS(EB2(4))
      GOTO 280
 270  EB2(4)=PI05
 280  CONTINUE
      IF(CVO.GE.0.0D0.AND.SVO.LT.0.0D0) EB2(4)=PI2-EB2(4)
      IF(CVO.LT.0.0D0.AND.SVO.GE.0.0D0) EB2(4)=PI-EB2(4)
      IF(CVO.LT.0.0D0.AND.SVO.LT.0.0D0) EB2(4)=PI+EB2(4)
C
      AKI=2.0D0/RZ-AV*AV/GS2
      EB2(2)=DSQRT(1.0D0-RZ*AKI*(2.0D0-RZ*AKI)*SHFI*SHFI)
      IF(DABS(EB2(2)-1.0D0).LT.1.0D-9) GOTO 290
      EB2(1)=(1.0D0-EB2(2))/AKI
      GOTO 300
 290  EB2(1)=RZ*SHFI*SHFI
 300  CF=(EB2(1)*(1.0D0+EB2(2))/RZ-1.0D0)/EB2(2)
      IF(CF.GE.0.999999999D0) GOTO 315
      IF(CF.LE.-0.999999999D0) GOTO 317
      SF=DSQRT(1.0D0-CF*CF)
      IF(DABS(CF).LT.1.0D-12) GOTO 310
      F=DATAN(SF/CF)
      IF(F.LT.0.0D0) F=PI+F
      GOTO 320
 310  F=PI05
      GOTO 320
 315  F=0.0D0
      GOTO 320
 317  F=PI
 320  CONTINUE
      IF(HFI.LT.PI05) F=PI2-F
C
      C2S=(X*CVO+Y*SVO)/RZ
      IF(DABS(CIB).LT.1.0D-9) GOTO 330
      S2S=(Y*CVO-X*SVO)/RZ/CIB
      GOTO 340
 330  S2S=Z/RZ/SIB
 340  CONTINUE
      IF(DABS(C2S).LT.1.0D-12) GOTO 350
      U2S=DATAN(S2S/C2S)
      U2S=DABS(U2S)
      GOTO 360
 350  U2S=PI05
 360  CONTINUE
      IF(C2S.GE.0.0D0.AND.S2S.LT.0.0D0) U2S=PI2-U2S
      IF(C2S.LT.0.0D0.AND.S2S.GE.0.0D0) U2S=PI-U2S
      IF(C2S.LT.0.0D0.AND.S2S.LT.0.0D0) U2S=PI+U2S
      EB2(3)=U2S-F
      IF(EB2(3).LT.0.0D0) EB2(3)=EB2(3)+PI2
C
      EB2(3)=EB2(3)/PI180
      EB2(4)=EB2(4)/PI180
      EB2(5)=EB2(5)/PI180
      CALL DCRIT(EA,EB2,DD2)
      FB=F/PI180
      R2(1)=R2(1)-0.18D3
      IF(R2(1).LT.0.0D0) R2(1)=R2(1)+0.36D3
      CALL RADIANT(EB2,FB,EZ,FZM,R2)
C
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE RADIANT(EB,FB,EZ,FZ,RR)
C ---------------------------------------------------------------------
C      - CALCULATION OF METEOR SHOWER RADIANT, IF THE FOLLOWING INPUT
C     QUANTITIES OF THE ORBIT CROSSING THE EARTH'S ORBIT ARE WELL-KNOWN:
C
C        EB(1) - PERIHELION DISTANCE [AU]
C        EB(2) - ECCENTRICITY
C        EB(3) - ARGUMENT OF PERIHELION [DEGREES]
C        EB(4) - LONGITUDE OF NODE [DEGREES]
C        EB(5) - INCLINATION TO ECLIPTIC [DEGREES]
C        FB - TRUE ANOMALY OF THE POINT OF INTERSECTION OF CALCULATED
C            SHOWER AND EARTH ORBITS [DEGREES]
C        EZ(1) - SEMI-MAJOR AXIS OF THE EARTH'S ORBIT [AU]
C        EZ(2) - ECCENTRICITY ...
C        EZ(3) - ARGUMENT OF PERIHELION ... [RADIANS]
C        EZ(4) - LONGITUDE OF NODE ... [RADIANS]
C        EZ(5) - INCLINATION TO ECLIPTIC ... [RADIANS]
C        FZ - TRUE ANOMALY OF THE EARTH IN THE POINT ... [RADIANS]
C
C     OUTPUT:
C     THE QUANTITIES DESCRIBING THE RADIANT OF POTENTIAL METEOR SHOWER
C        RR(2) - RIGHT ASCENSION OF THE PREDICTED RADIANT [DEGREES]
C        RR(3) - DECLINATION OF THE RADIANT [DEGREES]
C        RR(4) - PREDICTED GEOCENTRIC VELOCITY OF METEORS [KM/S]
C        RR(5) - PREDICTED HELIOCENTRIC VELOCITY OF METEORS [KM/S]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION RR(6),EB(5),EZ(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI05=PI/2.0D0
      PI2=2.0D0*PI
      PI180=PI/0.18D3
      GAUSS=0.17202098950D-1
      AUKM=0.149597870D9
      D24=8.6400D4
      EPS=2.3439291D1
C
      CEPS=DCOS(PI180*EPS)
      SEPS=DSIN(PI180*EPS)
      CHO=DCOS(PI180*EB(4))
      SHO=DSIN(PI180*EB(4))
      CSK=DCOS(PI180*EB(5))
      SSK=DSIN(PI180*EB(5))
      RK=EB(1)*(1.0D0+EB(2))/(1.0D0+EB(2)*DCOS(FB*PI180))
      CHZ=DCOS(EZ(4))
      SHZ=DSIN(EZ(4))
      CSZ=DCOS(EZ(5))
      SSZ=DSIN(EZ(5))
      RZ=EZ(1)*(1.0D0-EZ(2)*EZ(2))/(1.0D0+EZ(2)*DCOS(FZ))
C    
      VHR=PI180*FB
      CVH=DCOS(VHR)
      SVH=DSIN(VHR)
      WHFI=DSQRT(1.0D0+EB(2)*EB(2)+2.0D0*EB(2)*CVH)
      CHFI=-EB(2)*SVH/WHFI
      SHFI=(1.0D0+EB(2)*CVH)/WHFI
      QKONTR=SHFI*SHFI+CHFI*CHFI
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,100)
 100  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (ABS(SIN**2(HFI) + COS**
     *2(HFI)).GT.1)')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP
      IF(DABS(CHFI).LT.1.0D-12) HFI=PI05
      IF(DABS(CHFI).LT.1.0D-12) GOTO 120
      HFI=DATAN(SHFI/CHFI)
      HFI=DABS(HFI)
      IF(CHFI.GE.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI2-HFI
      IF(CHFI.LT.0.0D0.AND.SHFI.GE.0.0D0) HFI=PI-HFI
      IF(CHFI.LT.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI+HFI
 120  CONTINUE
      CFSO=DCOS(PI180*(FB+EB(3)))
      SFSO=DSIN(PI180*(FB+EB(3)))
      XK=RK*(CFSO*CHO-SFSO*CSK*SHO)
      YK=RK*(CFSO*SHO+SFSO*CSK*CHO)
      ZK=RK*SFSO*SSK
      C3=DCOS(PI180*(FB+EB(3))-HFI)
      S3=DSIN(PI180*(FB+EB(3))-HFI)
      AV=GAUSS*DSQRT(2.0D0/RK-(1.0D0-EB(2))/EB(1))
      VXE=-AV*(C3*CHO-S3*CSK*SHO)
      VYE=-AV*(C3*SHO+S3*CSK*CHO)
      VZE=-AV*S3*SSK
      VE=DSQRT(VXE*VXE+VYE*VYE+VZE*VZE)
      IF(DABS(AV-VE).GT.1.0D-12) WRITE(*,140)
 140  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (AV.NE.VE).')
      IF(DABS(AV-VE).GT.1.0D-12) STOP
      VXK=VXE
      VYK=VYE*CEPS-VZE*SEPS
      VZK=VYE*SEPS+VZE*CEPS
      VH=DSQRT(VXK*VXK+VYK*VYK+VZK*VZK)
      IF(DABS(AV-VH).GT.1.0D-12) WRITE(*,160)
 160  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (AV.NE.VH).')
      IF(DABS(AV-VH).GT.1.0D-12) STOP
C
      CVH=DCOS(FZ)
      SVH=DSIN(FZ)
      WHFI=DSQRT(1.0D0+EZ(2)*EZ(2)+2.0D0*EZ(2)*CVH)
      CHFI=-EZ(2)*SVH/WHFI
      SHFI=(1.0D0+EZ(2)*CVH)/WHFI
      QKONTR=SHFI*SHFI+CHFI*CHFI
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,200)
 200  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (ABS(SIN**2(HFI) + COS**
     *2(HFI)).GT.1)')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP
      IF(DABS(CHFI).LT.1.0D-12) HFI=PI05
      IF(DABS(CHFI).LT.1.0D-12) GOTO 220
      HFI=DATAN(SHFI/CHFI)
      HFI=DABS(HFI)
      IF(CHFI.GE.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI2-HFI
      IF(CHFI.LT.0.0D0.AND.SHFI.GE.0.0D0) HFI=PI-HFI
      IF(CHFI.LT.0.0D0.AND.SHFI.LT.0.0D0) HFI=PI+HFI
 220  CONTINUE
      CFSO=DCOS(FZ+EZ(3))
      SFSO=DSIN(FZ+EZ(3))
      XZ=RZ*(CFSO*CHZ-SFSO*CSZ*SHZ)
      YZ=RZ*(CFSO*SHZ+SFSO*CSZ*CHZ)
      ZZ=RZ*SFSO*SSZ
      C3=DCOS(FZ+EZ(3)-HFI)
      S3=DSIN(FZ+EZ(3)-HFI)
      AVZ=GAUSS*DSQRT(2.0D0/RZ-1.0D0/EZ(1))
      VXE=-AVZ*(C3*CHZ-S3*CSZ*SHZ)
      VYE=-AVZ*(C3*SHZ+S3*CSZ*CHZ)
      VZE=-AVZ*S3*SSZ
      VEZ=DSQRT(VXE*VXE+VYE*VYE+VZE*VZE)
      IF(DABS(AVZ-VEZ).GT.1.0D-12) WRITE(*,240)
 240  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (AVZ.NE.VEZ).')
      IF(DABS(AVZ-VEZ).GT.1.0D-12) STOP
      VXZ=VXE
      VYZ=VYE*CEPS-VZE*SEPS
      VZZ=VYE*SEPS+VZE*CEPS
      VHZ=DSQRT(VXZ*VXZ+VYZ*VYZ+VZZ*VZZ)
      IF(DABS(AVZ-VHZ).GT.1.0D-12) WRITE(*,260)
 260  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (AVZ.NE.VHZ).')
      IF(DABS(AVZ-VHZ).GT.1.0D-12) STOP
C
      VX=-(VXK-VXZ)
      VY=-(VYK-VYZ)
      VZ=-(VZK-VZZ)
      VG=DSQRT(VX*VX+VY*VY+VZ*VZ)
      IF(DABS(VZ).LT.1.0D-12) D=0.0D0
      IF(DABS(VZ).LT.1.0D-12) GOTO 310
      HELVEC=DSQRT(1.0D0-VZ*VZ/VG/VG)
      D=DATAN(VZ/HELVEC/VG)
 310  CONTINUE
      IF(DABS(D-PI05).LT.1.0D-12.OR.DABS(D+PI05).LT.1.0D-12) A=0.0D0
      IF(DABS(D-PI05).LT.1.0D-12.OR.DABS(D+PI05).LT.1.0D-12) GOTO 340
      CDVEC=DCOS(D)
      QQS=VY/VG/CDVEC
      QQC=VX/VG/CDVEC
      QKONTR=QQS*QQS+QQC*QQC
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) WRITE(*,320)
 320  FORMAT(/,'  ERROR IN SUBROUTINE `RADIANT` (SIN**2(A) + COS**2(A).G
     *T.1).')
      IF(DABS(QKONTR-1.0D0).GT.1.0D-12) STOP
      IF(DABS(QQC).LT.1.0D-12) A=PI05
      IF(DABS(QQC).LT.1.0D-12) GOTO 340
      A=DATAN(QQS/QQC)
      A=DABS(A)
      IF(QQC.GE.0.0D0.AND.QQS.LT.0.0D0) A=PI2-A
      IF(QQC.LT.0.0D0.AND.QQS.GE.0.0D0) A=PI-A
      IF(QQC.LT.0.0D0.AND.QQS.LT.0.0D0) A=PI+A
 340  CONTINUE
      RR(2)=A/PI180
      RR(3)=D/PI180
      RR(4)=VG*AUKM/D24
      RR(5)=AV*AUKM/D24
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE DCRIT(EA,EB,DD)
C ---------------------------------------------------------------------
C      - COMPUTATION OF D-DISCRIMINANT OF THE SOUTHWORTH-HAWKINS' D-CRI-
C     TERION EVALUATING SIMILARITY OF TWO ORBITS (SOUTHWORTH R.B.,
C     HAWKINS G.S.: 1963, SMITHSON. CONTR. ASTROPHYS. 7, 261.)
C
C INPUT:
C     EA(1), EB(1) - PERIHELION DISTANCES OF THE ORBITS [AU]
C     EA(2), EB(2) - ECCENTRICITIES
C     EA(3), EB(3) - ARGUMENTS OF PERIHELIA [DEGREES]
C     EA(4), EB(4) - LONGITUDES OF NODES [DEGREES]
C     EA(5), EB(5) - INCLINATIONS TO ECLIPTIC [DEGREES]
C
C OUTPUT:
C     DD - RESULTANT VALUE OF D-DISCRIMINANT
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION EA(6),EB(5)
C
      PI=4.0D0*DATAN(1.0D0)
      PI180=PI/0.18D3
C
      DQ=EA(1)-EB(1)
      DE=EA(2)-EB(2)
      DSO=(EA(3)-EB(3))*PI180
      DHO=(EA(4)-EB(4))*PI180
      DSK=(EA(5)-EB(5))*PI180
      SSKA=DSIN(EA(5)*PI180)
      SSKB=DSIN(EB(5)*PI180)
      POMSK=DSIN(DSK/2.0D0)
      POMHO=DSIN(DHO/2.0D0)
      DAI=4.0D0*POMSK*POMSK+SSKA*SSKB*4.0D0*POMHO*POMHO
      IF(DAI.GE.3.9999999D0.AND.DAI.LT.4.0000001D0) DAI=3.9999999D0
      IF(DAI.LT.4.0D0) GOTO 700
      WRITE(*,600)
 600  FORMAT(/,'  ERROR IN SUBROUTINE `DCRIT` (QUANTITY DAI.GE.4).')
      WRITE(*,*) DAI
      STOP
 700  SECSK=1.0D0/SQRT(1.0D0-DAI/4.0D0)
      ARGS=DCOS((EA(5)+EB(5))*PI180/2.0D0)*POMHO*SECSK
      JSGAR=-1
      IF(ARGS.GE.0.0D0) JSGAR=1
      POM=DABS(DHO)
      JSIGN=2
      IF(POM.GT.PI) JSIGN=-2
      POM=DABS(ARGS)
      IF(POM.LT.1.0000001D0) GOTO 720
      WRITE(*,620)
 620  FORMAT(/,'  ERROR IN SUBROUTINE `DCRIT` (ABS(ARGS).GT.1).')
      STOP
 720  CONTINUE
      IF(POM.LT.1) GOTO 740
      DBP=DSO+JSIGN*JSGAR*PI/2.0D0
      GOTO 760
 740  ODM=ARGS/DSQRT(1.0D0-ARGS*ARGS)
      DBP=DSO+JSIGN*DATAN(ODM)
 760  POM=DSIN(DBP/2.0D0)
      DD=DSQRT(DQ*DQ+DE*DE+DAI+(EA(2)+EB(2))*(EA(2)+EB(2))*POM*POM)
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE APPROACH(TMIN,TMAX,EA,WM1,WM2)
C ---------------------------------------------------------------------
C      - CALCULATION OF SOME QUANTITIES CHARACTERIZING THE MEAN ORBIT OF
C        THE METEOR SHOWER
C        
C INPUT:
C    "TMIN", "TMAX" - BEGINNING AND END OF YEAR OF INVESTIGATION OF THE
C       POTENTIAL SHOWER GIVEN IN JULIAN CENTURIES FROM 2000.0.
C    ORBITAL ELEMENTS OF THE ORBIT OF THE PARENT BODY
C       EA(1) - PERIHELION DISTANCE [AU]
C       EA(2) - ECCENTRICITY
C       EA(3) - ARGUMENT OF PERIHELION [DEGREES]
C       EA(4) - LONGITUDE OF NODE [DEGREES]
C       EA(5) - INCLINATION TO ECLIPTIC [DEGREES]
C       YEAR, MONTH, DAY - TIME OF PERIHELION PASSAGE OF THE PARENT BODY
C       (CONVERTED TO: EA(6) - THE TIME IN JULIAN CENTURIES FROM
C                              JD = 2451545.0)
C
C OUTPUT:
C   SOME QUANTITIES DESCRIBING THEORETICAL RADIANTS AT BOTH ARCS (PRE-
C   PERIHELION AND POST-PERIHELION) OF THE ORBIT OF THE PARENT BODY;
C   THE OUTPUT QUANTITIES AT THE POST-PERIHELION ARC ARE:
C      WM1(1) - MOMENT OF MINIMUM DISTANCE BETWEEN THE EARTH AND THE
C               ORBIT OF THE PARENT BODY [JULIAN CENTURIES FROM 2000.0]
C      WM1(2) - MINIMUM DISTANCE BETWEEN POST-PERIHELION ARC OF THE
C               PARENT BODY ORBIT AND EARTH'S ORBIT [AU]
C      WM1(3) - DISTANCE OF THE PARENT BODY FROM THE EARTH AT THE MOMENT
C               OF MAXIMUM [AU]
C      WM1(4) - TIME INTERVAL BETWEEN THE PASSAGES OF THE PARENT BODY
C               AND EARTH ACROSS THE NEAREST POINTS OF POST-PERIHELION
C               ARC OF THE PARENT BODY ORBIT AND EARTH'S ORBIT [DAYS]
C      WM1(5) - TRUE ANOMALY OF THE PARENT BODY IN THE MOMENT OF ITS
C               NEAREST POST-PERIHELION APPROACH TO THE EARTH'S ORBIT
C               [DEGREES]
C      WM1(6) - TRUE ANOMALY OF THE EARTH IN THE MOMENT OF ITS NEAREST
C               APPROACH TO THE POST-PERIHELION ARC OF THE PARENT BODY
C               ORBIT [DEGREES]
C   THE OUTPUT QUANTITIES AT THE PRE-PERIHELION ARC ARE ANALOGOUSLY
C   STORED IN FIELD "WM2(I)", WHERE "I = 1, 2, 3,..., 6".
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION EA(6),WM1(6),WM2(6),WEL(7)
C
      PI=4.0D0*DATAN(1.0D0)
      PI180=PI/0.18D3
      GAUSS=.17202098950D-1
C
      TMN1=(TMIN+TMAX)/2.0D0
      TMN2=TMN1
      DT0=(TMAX-TMIN)/3.6525D2
C
      CHO=DCOS(PI180*EA(4))
      SHO=DSIN(PI180*EA(4))
      CSK=DCOS(PI180*EA(5))
      SSK=DSIN(PI180*EA(5))
      PRM=EA(1)*(1.0D0+EA(2))
      IF(EA(2).GT.0.99999999D0) GOTO 30
      A=EA(1)/(1.0D0-EA(2))
      PP05=PI*A*DSQRT(A)/GAUSS
      PP=2.0D0*PP05
C
C  === 1-ST ARC ===
C
  30  VKD=0.0D0
      VKU=0.179D3
      DVK=1.0D0
      TA=TMIN
      TB=TMAX
      DT=DT0
C
  50  WM1(2)=1.111D3
      DO 200 VK=VKD,VKU,DVK
      WEN=1.0D0+EA(2)*DCOS(PI180*VK)
      RK=1.0D6
      IF(DABS(WEN).LT.1.0D-12) GOTO 100
      RK=PRM/WEN
 100  CVS=DCOS(PI180*(VK+EA(3)))
      SVS=DSIN(PI180*(VK+EA(3)))
      XK=RK*(CVS*CHO-SVS*CSK*SHO)
      YK=RK*(CVS*SHO+SVS*CSK*CHO)
      ZK=RK*SVS*SSK
C
      DO 140 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ))
      IF(DR.GT.1.0D2) GOTO 240
      IF(DR.GE.WM1(2)) GOTO 140
      WM1(2)=DR
      WM1(5)=VK
      WM1(6)=WEL(7)
      TMN1=T20
 140  CONTINUE
 200  CONTINUE
C      BOTH LABELS 200 AND 240 ARE NECESSARY!
 240  CONTINUE
      IF(DVK.GT.0.5D0) GOTO 245
      IF(DABS(WM1(2)-RMOLD).LT.1.0D-5) GOTO 260
 245  VKD=WM1(5)-1.0D0*DVK
      VKU=WM1(5)+1.0D0*DVK
      TA=TMN1-1.0D0*DT
      TB=TMN1+1.0D0*DT
      DVK=0.1D-1*DVK
      DT=0.1D-1*DT
      RMOLD=WM1(2)
      GOTO 50
C
 260  CONTINUE
C
      IF(WM1(2).GT.1.0D2) GOTO 530
      DTM=3.6525D4*(TMN1-EA(6))
      EQX=EA(1)
      EEX=EA(2)
      CALL KEPLER(EQX,EEX,DTM,QWV)
      WENK=1.0D0+EA(2)*DCOS(QWV)
      RK=PRM/WENK
      CVS=DCOS(QWV+PI180*EA(3))
      SVS=DSIN(QWV+PI180*EA(3))
      XK=RK*(CVS*CHO-SVS*CSK*SHO)
      YK=RK*(CVS*SHO+SVS*CSK*CHO)
      ZK=RK*SVS*SSK
      CALL ECLXYZ(TMN1,XZ,YZ,ZZ,WEL)
      DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ))
      WM1(3)=DR
C
      WM1(4)=DTM
      QWV=WM1(5)*PI180
      CALL TIMEANOM(EQX,EEX,QWV,DTT)
      WM1(4)=WM1(4)-DTT
      IF(EA(2).GT.0.99999999D0) GOTO 270
 263  CONTINUE
      IF(WM1(4).LE.PP05) GOTO 265
      WM1(4)=WM1(4)-PP
      GOTO 263
 265  CONTINUE
      IF(WM1(4).GE.-PP05) GOTO 270
      WM1(4)=WM1(4)+PP
      GOTO 265
C
 270  WM1(1)=TMN1
C
C  === 2-ND ARC ===
C
 530  VKD=0.359D3
      VKU=0.181D3
      DVK=-1.0D0
      TA=TMIN
      TB=TMAX
      DT=DT0
C
 550  WM2(2)=1.111D3
      DO 700 VK=VKD,VKU,DVK
      WEN=1.0D0+EA(2)*DCOS(PI180*VK)
      RK=1.0D6
      IF(DABS(WEN).LT.1.0D-12) GOTO 600
      RK=PRM/WEN
 600  CVS=DCOS(PI180*(VK+EA(3)))
      SVS=DSIN(PI180*(VK+EA(3)))
      XK=RK*(CVS*CHO-SVS*CSK*SHO)
      YK=RK*(CVS*SHO+SVS*CSK*CHO)
      ZK=RK*SVS*SSK
C
      DO 640 T20=TA,TB,DT
      CALL ECLXYZ(T20,XZ,YZ,ZZ,WEL)
      DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ))
      IF(DR.GT.1.0D2) GOTO 740
      IF(DR.GE.WM2(2)) GOTO 640
      WM2(2)=DR
      WM2(5)=VK
      WM2(6)=WEL(7)
      TMN2=T20
 640  CONTINUE
 700  CONTINUE
C      BOTH LABELS 700 AND 740 ARE NECESSARY!
 740  CONTINUE
      IF(DVK.GT.0.5D0) GOTO 745
      IF(DABS(WM2(2)-RMOLD).LT.1.0D-5) GOTO 760
 745  VKD=WM2(5)-1.0D0*DVK
      VKU=WM2(5)+1.0D0*DVK
      TA=TMN2-1.0D0*DT
      TB=TMN2+1.0D0*DT
      DVK=0.1D-1*DVK
      DT=0.1D-1*DT
      RMOLD=WM2(2)
      GOTO 550
C
 760  CONTINUE
C
      IF(WM2(2).GT.1.0D2) GOTO 800
      DTM=3.6525D4*(TMN2-EA(6))
      EQX=EA(1)
      EEX=EA(2)
      CALL KEPLER(EQX,EEX,DTM,QWV)
      WENK=1.0D0+EA(2)*DCOS(QWV)
      RK=PRM/WENK
      CVS=DCOS(QWV+PI180*EA(3))
      SVS=DSIN(QWV+PI180*EA(3))
      XK=RK*(CVS*CHO-SVS*CSK*SHO)
      YK=RK*(CVS*SHO+SVS*CSK*CHO)
      ZK=RK*SVS*SSK
      CALL ECLXYZ(TMN2,XZ,YZ,ZZ,WEL)
      DR=DSQRT((XK-XZ)*(XK-XZ)+(YK-YZ)*(YK-YZ)+(ZK-ZZ)*(ZK-ZZ))
      WM2(3)=DR
C
      WM2(4)=DTM
      QWV=WM2(5)*PI180
      CALL TIMEANOM(EQX,EEX,QWV,DTT)
      WM2(4)=WM2(4)-DTT
      IF(EA(2).GT.0.99999999D0) GOTO 770
 763  CONTINUE
      IF(WM2(4).LE.PP05) GOTO 765
      WM2(4)=WM2(4)-PP
      GOTO 763
 765  CONTINUE
      IF(WM2(4).GE.-PP05) GOTO 770
      WM2(4)=WM2(4)+PP
      GOTO 765
C
 770  WM2(1)=TMN2
 800  CONTINUE
      RETURN
C
      END
C ----------------------------------------------------------------------
      SUBROUTINE ECLXYZ(T20,X,Y,Z,WEL)
C ----------------------------------------------------------------------
C     - CALCULATION OF QUANTITIES CHARACTERIZING THE POSITION OF THE
C       EARTH AND ITS ORBIT AT GIVEN TIME (RECTANGULAR COORDINATES OF
C       THE EARTH WERE CALCULATED BY: BRETAGNON P.: 1982, ASTRON.
C       ASTROPHYS. 114, 278.)
C
C INPUT:
C    T20 - THE TIME [JULIAN CENTURIES FROM JD = 2451545.0]
C
C OUTPUT = THE VALUES VALID IN TIME "T20" (EQUINOX 2000.0):
C    X, Y, Z - HELIOCENTRIC RECTANGULAR COORDINATES REFERRED TO THE
C              ECLIPTIC [AU]
C    WEL(1) - SEMI-MAJOR AXIS [AU]
C    WEL(2) - ECCENTRICITY
C    WEL(3) - ECLIPTIC MEAN LONGITUDE [RADIANS]
C    WEL(4) - LONGITUDE OF PERIHELION [RADIANS]
C    WEL(5) - LONGITUDE OF ASCENDING NODE [RADIANS]
C    WEL(6) - INCLINATION TO THE ECLIPTIC [RADIANS]
C    WEL(7) - TRUE ANOMALY [RADIANS]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (A,C,E,P-Z)
      DIMENSION WEL(7)
C
C
      PI=4.0D0*DATAN(1.0D0)
      PI2=2.0D0*PI
C
      WEL(1)=1.00000101778D0
C
      WEL(3)=1.75347031435D0+6.28307584918D2*T20
      WEL(3)=WEL(3)-9.9189D-8*T20*T20+7.3D-13*T20*T20*T20
   3  CONTINUE
      IF(WEL(3).GE.0.0D0) GOTO 5
      WEL(3)=WEL(3)+PI2
      GOTO 3
   5  CONTINUE
      IF(WEL(3).LT.PI2) GOTO 7
      WEL(3)=WEL(3)-PI2
      GOTO 5
C
   7  WK=-3.7408165D-3-8.2266699D-5*T20
      WK=WK+2.748939D-7*T20*T20+1.04217D-9*T20*T20*T20
      WH=1.628447663D-2-6.2030259D-5*T20
      WH=WH-3.353888D-7*T20*T20+7.1185D-10*T20*T20*T20
      WQ=-1.13469002D-4*T20+1.237314D-7*T20*T20+1.2705D-9*T20*T20*T20
      WP=1.0180391D-5*T20+4.701998D-7*T20*T20-5.3829D-10*T20*T20*T20
C 
      WEL(2)=DSQRT(WK*WK+WH*WH)
C
      IF(DABS(WK).LT.1.0D-12) GOTO 210
      WEL(4)=DATAN(WH/WK)
      WEL(4)=DABS(WEL(4))
      IF(WH.GE.0.0D0.AND.WK.LT.0.0D0) WEL(4)=PI-WEL(4)
      IF(WH.LT.0.0D0.AND.WK.LT.0.0D0) WEL(4)=PI+WEL(4)
      IF(WH.LT.0.0D0.AND.WK.GE.0.0D0) WEL(4)=PI2-WEL(4)
C
  10  CONTINUE
      IF(DABS(WQ).LT.1.0D-12) GOTO 220
      WEL(5)=DATAN(WP/WQ)
      WEL(5)=DABS(WEL(5))
      IF(WP.GE.0.0D0.AND.WQ.LT.0.0D0) WEL(5)=PI-WEL(5)
      IF(WP.LT.0.0D0.AND.WQ.LT.0.0D0) WEL(5)=PI+WEL(5)
      IF(WP.LT.0.0D0.AND.WQ.GE.0.0D0) WEL(5)=PI2-WEL(5)
  20  WG=DSQRT(WP*WP+WQ*WQ)
      POM=1.0D0-WG*WG
      IF(DABS(POM).LT.1.0D-12) GOTO 230
      POM=WG/DSQRT(POM)
      WEL(6)=2.0D0*DATAN(POM)
      IF(WEL(6).LT.0.0D0) WEL(6)=PI+WEL(6)
C
  30  CONTINUE
      WM=WEL(3)-WEL(4)
      WEL(7)=WM+(2.0D0-WEL(2)*WEL(2)/4.0D0)*WEL(2)*DSIN(WM)
      WEL(7)=WEL(7)+5.0D0*WEL(2)*WEL(2)/4.0D0*DSIN(2.0D0*WM)
      WEL(7)=WEL(7)+1.3D1*WEL(2)*WEL(2)*WEL(2)/1.2D1*DSIN(3.0D0*WM)
  40  CONTINUE
      IF(WEL(7).GE.0.0D0) GOTO 50
      WEL(7)=WEL(7)+PI2
      GOTO 40
  50  CONTINUE
      IF(WEL(7).LT.PI2) GOTO 60
      WEL(7)=WEL(7)-PI2
      GOTO 50
  60  RZ=WEL(1)*(1.0D0-WEL(2)*WEL(2))/(1.0D0+WEL(2)*DCOS(WEL(7)))
      W2=WEL(7)+WEL(4)-WEL(5)
      C2O=DCOS(W2)
      S2O=DSIN(W2)
      CVO=DCOS(WEL(5))
      SVO=DSIN(WEL(5))
      CSK=DCOS(WEL(6))
      SSK=DSIN(WEL(6))
      X=RZ*(C2O*CVO-S2O*CSK*SVO)
      Y=RZ*(C2O*SVO+S2O*CSK*CVO)
      Z=RZ*S2O*SSK
      RETURN
C
 210  WEL(4)=PI/2.0D0
      IF(WH.LT.0.0D0) WEL(4)=1.5D0*PI
      GOTO 10
 220  WEL(5)=PI/2.0D0
      IF(WP.LT.0.0D0) WEL(5)=1.5D0*PI
      GOTO 20
 230  WEL(6)=PI
      GOTO 30
C
      END
C --------------------------------------------------------------------
      SUBROUTINE JULIAN(IROK,IMES,DEN,T19,T20)
C --------------------------------------------------------------------
C     - CALCULATION OF JULIAN DATE FROM CIVIL DATE
C
C INPUT:
C    IROK, IMES, DEN - YEAR, MONTH, AND DAY OF THE DATE
C
C OUTPUT:
C    T19 - GIVEN TIME ACCOUNTED FROM THE BEGINNING OF EPOCH 1900.0
C          (JD = 2415020.0) AND EXPRESSED IN JULIAN CENTURIES
C    T20 - GIVEN TIME ACCOUNTED FROM THE BEGINNING OF EPOCH 2000.0
C          (JD = 2451545.0) AND EXPRESSED IN JULIAN CENTURIES
C
C    JULIAN DATE CAN BE CALCULATED FROM "T20" AS:
C    JD=36525.0*T20+2451545.0
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (A-H,K-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION JMESI(13)
      DATA JMESI/0,31,28,31,30,31,30,31,31,30,31,30,31/
C
      IF(IROK.GT.1582) GOTO 101
      IF(IROK.EQ.1582) GOTO 102
      IF(IROK.GT.0.AND.IROK.LT.1582) GOTO 103
      IF(IROK.GE.-4713.AND.IROK.LT.0) GOTO 104
C
      T20=(IROK+4713)*365
      T20=T20-INT((-IROK-4712)/4)
  100 T0JD=-0.5D0
      GOTO 210
  101 T20=(IROK-1583)*365
      T20=T20+INT((IROK-1581)/4)
      T20=T20-INT((IROK-1501)/100)+INT((IROK-1201)/400)
      T0JD=2.2992385D6
      GOTO 210
  102 T20=0.0D0
      IF(IMES.EQ.10.AND.DEN.GT.1.49999999D1) T20=-1.0D1
      IF(IMES.GT.10) T20=-1.0D1
      T0JD=2.2988835D6
      GOTO 210
  103 T20=(IROK-1)*365
      T20=T20+INT((IROK-1)/4)
      T0JD=1.7214235D6
      GOTO 210
  104 T20=(IROK+4713)*365
      T20=T20+INT((IROK+4715)/4)
      T0JD=-0.5D0
  210 YY=IROK/1.00D2-INT(IROK/100)
      ZZ=IROK/4.00D2-INT(IROK/400)
      IF(IROK.GT.1582.AND.YY.EQ.0.0D0.AND.ZZ.GT.0.0D0) GOTO 220
      YY=IROK/4.0D0-INT(IROK/4)
      IF(YY.EQ.0.0D0) JMESI(3)=29
      IF(IROK.EQ.-1) JMESI(3)=29
  220 CONTINUE
      DO 225 JJM=1,IMES-1
      T20=T20+JMESI(JJM+1)
  225 CONTINUE
      T20=T20+DEN-1.0D0+T0JD
      JMESI(3)=28
C
      T19=(T20-2.415020D6)/3.6525D4
      T20=(T20-2.451545D6)/3.6525D4
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE INVJUL(TJD,IRK,IMS,DNI)
C ---------------------------------------------------------------------
C     - CALCULATION OF CIVIL DATE FROM JULIAN DATE (INVERSION OPERATION
C       TO THAT EXECUTED BY SUBROUTINE "JULIAN")
C
C INPUT:
C    TJD - TIME FROM THE BEGINNING OF EPOCH 2000.0 (JD = 2451545.0)
C          IN JULIAN CENTURIES (IT CAN BE CALCULATED FROM JULIAN DATE
C          "JD" IN DAYS AS: TJD=(JD-2451545.0)/36525.0)
C
C OUTPUT:
C    IRK, IMS, DNI - YEAR, MONTH, AND DAY OF THE CORRESPONDING CIVIL
C                    DATE
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (D,T)
      IMPLICIT INTEGER (I)
      DIMENSION IMESI(13)
      DATA IMESI/0,31,28,31,30,31,30,31,31,30,31,30,31/
C
      DNI=TJD*3.6525D4+2.451545D6+1.0D0
      IF(DNI.LT.0.5D0) GOTO 16
      IF(DNI.GE.0.5D0.AND.DNI.LT.1.7214245D6) GOTO 15
      IF(DNI.GE.1.7214245D6.AND.DNI.LT.2.2991615D6) GOTO 14
      IF(DNI.GE.2.2991615D6.AND.DNI.LT.2.2992395D6) GOTO 13
      IF(DNI.GE.2.2992395D6.AND.DNI.LT.2.4515455D6) GOTO 12
C
C  FROM 2000-JAN-1.0
  11  CONTINUE
      DNI=DNI-2.4515445D6
      IRK=INT(DNI/3.65242198D2)
      DN0=365*IRK+INT((IRK+3)/4.0D0)-INT((IRK-1)/1.0D2)+INT((IRK-1)/4.0D
     *2)
      DNI=DNI-DN0
      IF(DNI.LT.1.0D0) GOTO 30
      IRK=IRK+2000
      GOTO 20
  30  IRK=IRK-1
      DNI=DNI+3.65D2
      DH4=IRK/4.0D0-INT(IRK/4.0D0)
      DH100=IRK/1.0D2-INT(IRK/1.0D2)
      IF(DABS(DH4).LT.1.0D-6.AND.DH100.GT.0.0D0) DNI=DNI+1.0D0
      DH400=IRK/4.0D2-INT(IRK/4.0D2)
      IF(DABS(DH400).LT.1.0D-6) DNI=DNI+1.0D0
      IRK=IRK+2000
      GOTO 20
C
C   FROM 1583-JAN-1.0 TO 2000-JAN-1.0
  12  CONTINUE 
      DNI=DNI-2.2992385D6
      IRK=INT(DNI/3.65242198D2)
      DN0=365*IRK+INT((IRK+2)/4.0D0)-INT((IRK+82)/1.0D2)+INT((IRK+382)/4
     *.0D2)
      DNI=DNI-DN0
      IF(DNI.LT.1.0D0) GOTO 35
      IRK=IRK+1583
      GOTO 20
  35  IRK=IRK-1
      DH4=(IRK+3)/4.0D0-INT((IRK+3)/4.0D0)
      DH100=(IRK+83)/1.0D2-INT((IRK+83)/1.0D2)
      DNI=DNI+3.65D2
      IF(DABS(DH4).LT.1.0D-6.AND.DH100.GT.0.0D0) DNI=DNI+1.0D0
      DH400=(IRK+383)/4.0D2-INT((IRK+383)/4.0D2)
      IF(DABS(DH400).LT.1.0D-6) DNI=DNI+1.0D0
      IRK=IRK+1583
      GOTO 20
C
C   FROM 1582-OCT-15.0 TO 1583-JAN-1.0
  13  CONTINUE 
      IRK=1582
      DNI=DNI-2.2991605D6+1.4D1
      IF(DNI.GE.3.2D1) GOTO 40
      IMS=10
      GOTO 26
  40  CONTINUE
      IF(DNI.GE.6.2D1) GOTO 45
      DNI=DNI-3.1D1
      IMS=11
      GOTO 26
  45  DNI=DNI-6.1D1
      IMS=12
      GOTO 26
C
C   FROM 1-JAN-1.0 TO 1582-OCT-15.0
  14  CONTINUE
      DNI=DNI-1.7214235D6
      IRK=INT(DNI/3.6525D2)+1
      DN0=3.65D2*(IRK-1)+INT((IRK-1)/4.0D0)
      DNI=DNI-DN0
      IF(DNI.GE.1.0D0) GOTO 20
      IRK=IRK-1
      DH4=IRK/4.0D0-INT(IRK/4.0D0)
      DNI=DNI+3.65D2
      IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0
      GOTO 20
C
C   FROM -4713-JAN-1.0 TO 1-JAN-1.0
  15  CONTINUE 
      DNI=DNI+0.5D0
      IRK=INT(DNI/3.6525D2)
      DN0=365*IRK+INT((IRK+2)/4.0D0)
      IF(DNI.GT.1.721118D6) DNI=DNI-1.0D0
      DNI=DNI-DN0
      IF(DNI.LT.1.0D0) GOTO 50
      IRK=IRK-4713
      GOTO 20
  50  IRK=IRK-1
      DH4=(IRK+3)/4.0D0-INT((IRK+3)/4.0D0)
      DNI=DNI+3.65D2
      IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0
      IRK=IRK-4713
      GOTO 20
C
C   TO -4713-JAN-1.0
  16  CONTINUE
      DNI=DNI+1.9307115D6
      IRK=INT(DNI/3.6525D2)
      DN0=365*IRK+INT(IRK/4.0D0)
      DNI=DNI-DN0
      IF(DNI.LT.1.0D0) GOTO 55
      IRK=IRK-9999
      GOTO 20
  55  IRK=IRK-1
      DH4=(IRK+1)/4.0D0-INT((IRK+1)/4.0D0)
      DNI=DNI+3.65D2
      IF(DABS(DH4).LT.1.0D-6) DNI=DNI+1.0D0
      IRK=IRK-9999
C
  20  CONTINUE
      DH4=IRK/4.0D0-INT(IRK/4.0D0)
      DH100=IRK/1.0D2-INT(IRK/1.0D2)
      DH400=IRK/4.0D2-INT(IRK/4.0D2)
      IF(DABS(DH4).LT.1.0D-6) IMESI(3)=29
      IF(IRK.GT.1582.AND.DABS(DH100).LT.1.0D-6) IMESI(3)=28
      IF(IRK.GT.1582.AND.DABS(DH400).LT.1.0D-6) IMESI(3)=29
      IF(IRK.EQ.0) IRK=-1
      IMS=1
  23  CONTINUE
      IF(IMS.NE.13) GOTO 60
      IMS=12
      DNI=DNI-1.0D0
  60  CONTINUE
      IF(IMESI(IMS+1)+1.GT.DNI) GOTO 26
      DNI=DNI-IMESI(IMS+1)
      IMS=IMS+1
      GOTO 23
  26  CONTINUE
      IMESI(3)=28
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE TIMEANOM(EQX,EEX,QWV,DTT)
C ---------------------------------------------------------------------
C     - CALCULATION OF THE TIME ELAPSED FROM THE MOMENT OF PERIHELION
C       PASSAGE, IF TRUE ANOMALY IS WELL-KNOWN
C
C INPUT:
C    EQX - PERIHELION DISTANCE [AU]
C    EEX - ECCENTRICITY
C    QWV - TRUE ANOMALY [RADIANS; IN INTERVAL FROM 0 TO 2*PI]
C
C OUTPUT:
C    DTT - THE TIME ELAPSED FROM THE MOMENT OF PERIHELION PASSAGE
C          (NEGATIVE AT PRE-PERIHELION ARC) [DAYS]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (D,E,G,P,Q,T)
C
      PI=4.0D0*DATAN(1.0D0)
      PI2=2.0D0*PI
      GAUSS=1.7202098950D-2
C
      IF(QWV.GE.PI) QWV=QWV-PI2
 300  TGV=DTAN(QWV/2.0D0)
      IF(DABS(EEX-1.0D0).LT.1.0D-9) GOTO 400
      QA=EQX/DABS(1.0D0-EEX)
      QNI=GAUSS/QA/DSQRT(QA)
      IF(EEX.GT.1.0D0) GOTO 500
      TGE=DSQRT((1.0D0-EEX)/(1.0D0+EEX))*TGV
      QEE=2.0D0*DATAN(TGE)
      DTT=(QEE-EEX*DSIN(QEE))/QNI
      GOTO 550
 400  DTT=EQX*DSQRT(2.0D0*EQX)*(1.0D0+TGV*TGV/3.0D0)*TGV/GAUSS
      GOTO 550
 500  THF=DSQRT((EEX-1.0D0)/(EEX+1.0D0))*TGV
      QFF=DLOG((1.0D0+THF)/(1.0D0-THF))
      DTT=(EEX*(DEXP(QFF)-1.0D0/DEXP(QFF))/2.0D0-QFF)/QNI
 550  CONTINUE
      IF(QWV.LT.0.0D0) QWV=QWV+PI2
      RETURN
C
      END
C ---------------------------------------------------------------------
      SUBROUTINE KEPLER(EQX,EEX,DTM,QWV)
C ---------------------------------------------------------------------
C     - SOLUTION OF THE KEPLER'S EQUATION (CALCULATION OF TRUE ANOMALY,
C       IF THE TIME ELAPSED FROM THE MOMENT OF PERIHELION PASSAGE IS
C       WELL-KNOWN)
C
C INPUT:
C    EQX - PERIHELION DISTANCE [AU]
C    EEX - ECCENTRICITY
C    DTM - THE TIME ELAPSED FROM THE MOMENT OF PERIHELION PASSAGE
C          (NEGATIVE AT PRE-PERIHELION ARC) [DAYS]
C          
C
C OUTPUT:
C    QWV - TRUE ANOMALY [RADIANS; IN INTERVAL FROM 0 TO 2*PI]
C
C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C
      IMPLICIT REAL*8 (D,E,G,P,Q,U)
C
      PI=4.0D0*DATAN(1.0D0)
      GAUSS=0.17202098950D-1
C
      IF(DABS(EEX-1.0D0).LT.1.0D-9) GOTO 730
      URN=EQX/DABS(1.0D0-EEX)
      URN=GAUSS/URN/DSQRT(URN)
      EMX=URN*DTM
      URD0=EMX
      IF(DABS(URD0).GT.5.0D2) URD0=URD0*5.0D2/DABS(URD0)
      IF(EEX.GT.1.0D0) GOTO 720
 710  URPOM=1.0D0-EEX*DCOS(URD0)
      IF(DABS(URPOM).LT.1.0D-15) URPOM=URPOM+1.0D-15
      URD=URD0+(EMX-URD0+EEX*DSIN(URD0))/URPOM
      IF(DABS(URD-URD0).LT.1.0D-9) GOTO 715
      URD0=URD
      GOTO 710
 715  CONTINUE
      IF(DABS(URD-PI).LT.1.0D-9.OR.DABS(URD+PI).LT.1.0D-9) GOTO 717
      URPOM=DSQRT((1.0D0+EEX)/(1.0D0-EEX))*DTAN(URD/2.0D0)
      QWV=2.0D0*DATAN(URPOM)
      IF(QWV.LT.0.0D0) QWV=2.0D0*PI+QWV
      GOTO 790
 717  QWV=PI
      GOTO 790
 720  URPOM=EEX*(DEXP(URD0)+1.0D0/DEXP(URD0))/2.0D0-1.0D0
      IF(DABS(URPOM).LT.1.0D-15) URPOM=URPOM+1.0D-15
      URD=(EMX-EEX*(DEXP(URD0)-1.0D0/DEXP(URD0))/2.0D0+URD0)/URPOM
      URD=URD0+URD
      IF(DABS(URD-URD0).LT.1.0D-9) GOTO 725
      URD0=URD
      GOTO 720
 725  URPOM=DSQRT((EEX+1.0D0)/(EEX-1.0D0))
      URPOM=URPOM*(DEXP(URD/2.0D0)-1.0D0/DEXP(URD/2.0D0))
      URPOM=URPOM/(DEXP(URD/2.0D0)+1.0D0/DEXP(URD/2.0D0))
      QWV=2.0D0*DATAN(URPOM)
      IF(QWV.LT.0.0D0) QWV=2.0D0*PI+QWV
      GOTO 790
 730  URN=GAUSS/EQX/DSQRT(8.0D0*EQX)
      URD0=0.0D0
 740  URD=6.0D0*URN*DTM-3.0D0*URD0-URD0*URD0*URD0
      URD=URD0+URD/(3.0D0+3.0D0*URD0*URD0)
      IF(DABS(URD-URD0).LT.1.0D-9) GOTO 745
      URD0=URD
      GOTO 740
 745  QWV=2.0D0*DATAN(URD)
      IF(QWV.LT.0.0D0) QWV=2.0D0*PI+QWV
 790  CONTINUE
      RETURN
C
      END