      SUBROUTINE REACH                                                  REAC001 
	character*8 id,ndate,modx,mody,modz,modrin,koll,modsav
C          CALCULATES THE STRAIGHT LINE RAY PATH BETWEEN THE EARTH      REAC002 
C          AND THE IONOSPHERE OR BETWEEN IONOSPHERIC LAYERS             REAC003 
      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       REAC004 
      COMMON /TRAC/ GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT               REAC005 
      COMMON /COORD/ S                                                  REAC006 
      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR,     REAC007 
     1             LPOLAR                                               REAC008 
      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  REAC009 
      COMMON R(20),T,STP,DRDT(20) /WW/ ID(10),W0,W(400)                 REAC010 
      EQUIVALENCE (EARTHR,W(2)),(XMTRH,W(3)),(RCVRH,W(20))              REAC011 
      LOGICAL CROSS,CROSSG,CROSSR,SPACE,GROUND,PERIGE,THERE,MINDIS,     REAC012 
     1        NEWRAY,RSPACE                                             REAC013 
      REAL N2,N2I                                                       REAC 14 
      COMPLEX PNP,POLAR,LPOLAR                                          REAC 15 
	data nstep / 500 /
      CALL HAMLTN                                                       REAC 17 
      H=R(1)-EARTHR                                                     REAC 18 
      IF (.NOT.NEWRAY.AND..NOT.RSPACE) CALL PRINTR(8HEXIT ION,0.)       REAC 19 
      NEWRAY=.FALSE.                                                    REAC 20 
      V=SQRT(R(4)**2+R(5)**2+R(6)**2)                                   REAC 21 
C********* NORMALIZE THE WAVE NORMAL DIRECTION TO ONE                   REAC 22 
      R(4)=R(4)/V                                                       REAC 23 
      R(5)=R(5)/V                                                       REAC 24 
      R(6)=R(6)/V                                                       REAC 25 
C********* NEGATIVE OF DISTANCE ALONG RAY TO CLOSEST APPROACH TO CENTER REAC 26 
C********* OF EARTH                                                     REAC 27 
      UP=R(1)*R(4)                                                      REAC 28 
      RADG=EARTHR**2-R(1)**2*(R(5)**2+R(6)**2)                          REAC 29 
      DISTG=SQRT (AMAX1(0.,RADG))                                       REAC 30 
C********* DISTANCE ALONG RAY TO FIRST INTERSECTION WITH OR CLOSEST     REAC 31 
C********* APPROACH TO THE EARTH                                        REAC 32 
      SG=-UP-DISTG                                                      REAC 33 
C********* CROSSG IS TRUE IF THE RAY WILL INTERSECT OR TOUCH THE EARTH  REAC 34 
      CROSSG=UP.LT.0..AND.RADG.GE.0.                                    REAC 35 
      RADR=(EARTHR+RCVRH)**2-R(1)**2*(R(5)**2+R(6)**2)                  REAC 36 
      DISTR=SQRT (AMAX1(0.,RADR))                                       REAC 37 
C********* DISTANCE ALONG RAY TO THE FIRST INTERSECTION WITH OR CLOSEST REAC 38 
C********* APPROACH TO THE RECEIVER HEIGHT                              REAC 39 
      SR=DISTR-UP                                                       REAC 40 
      IF (UP.LT.0..AND.DISTR.LT.-UP.AND.R(1).NE.EARTHR+RCVRH) SR=-DISTR REAC 41 
     1 -UP                                                              REAC 42 
C********* CROSSR IS TRUE IF THE RAY WILL INTERSECT WITH OR MAKE A      REAC 43 
C********* CLOSEST APPROACH TO THE RECEIVER HEIGHT                      REAC 44 
      CROSSR=R(4).LT.0..OR.R(1).LT.(EARTHR+RCVRH)                       REAC 45 
      CROSS=CROSSG.OR.CROSSR                                            REAC 46 
C********* MAXIMUM DISTANCE IN WHICH TO LOOK FOR THE IONOSPHERE         REAC 47 
      S1=AMIN1(SR,SG)                                                   REAC 48 
      IF(.NOT.CROSSG) S1=SR                                             REAC 49 
      IF (UP.GE.0.) GO TO 15                                            REAC 50 
      CROSS=.TRUE.                                                      REAC 51 
C********* IF RAY IS GOING DOWN, S1 IS AT MOST THE DISTANCE TO A PERIGEEREAC 52 
      S1=AMIN1(S1,-UP)                                                  REAC 53 
C********* CONVERT THE POSITION AND DIRECTION OF THE RAY TO CARTESIAN   REAC 54 
C********* COORDINATES                                                  REAC 55 
   15 CALL POL CAR                                                      REAC 56 
      SSTEP=100.                                                        REAC 57 
      S=SSTEP                                                           REAC 58 
      DO 20 I=1,NSTEP                                                   REAC 59 
      IF ((S-SSTEP).GT.S1.AND.CROSS) GO TO 25                           REAC 60 
C********* CONVERT POSITION AND DIRECTION TO SPHERICAL POLAR COORDINATESREAC 61 
C********* AT A DISTANCE S ALONG THE RAY                                REAC 62 
      CALL CAR POL                                                      REAC 63 
      CALL ELECTX                                                       REAC 64 
C********* FREE SPACE                                                   REAC 65 
      IF (X.EQ.0.) GO TO 20                                             REAC 66 
      CALL RINDEX                                                       REAC 67 
C********* EFFECTIVELY FREE SPACE                                       REAC 68 
      IF (SPACE) GO TO 20                                               REAC 69 
      IF (SSTEP.LT.0.5E-4) GO TO 25                                     REAC 70 
C********* RAY IN THE IONOSPHERE.  STEP BACK OUT                        REAC 71 
      S=S-SSTEP                                                         REAC 72 
C********* DECREASE STEP SIZE                                           REAC 73 
      SSTEP=SSTEP/10.                                                   REAC 74 
   20 S=S+SSTEP                                                         REAC 75 
      PRINT 2000, NSTEP                                                 REAC 76 
 2000 FORMAT (9H EXCEEDED,I5,26H STEPS IN SUBROUTINE REACH)             REAC 77 
      CALL EXIT                                                         REAC 78 
   25 IF(CROSS) S=AMIN1(S,S1)                                           REAC 79 
C********* CONVERT POSITION AND DIRECTION TO SPHERICAL POLAR COORDINATESREAC 80 
C********* AT A DISTANCE S ALONG THE RAY                                REAC 81 
      CALL CAR POL                                                      REAC 82 
C********* AVOID THE RAY BEING SLIGHTLY UNDERGROUND                     REAC 83 
      R(1)=AMAX1(R(1),EARTHR)                                           REAC 84 
C********* ONE STEP INTEGRATION                                         REAC 85 
      IF (N.LT.7) GO TO 31                                              REAC 86 
      DO 30 NN=7,N                                                      REAC 87 
   30 R(NN)=R(NN)+S*DRDT(NN)                                            REAC 88 
   31 T=T+S                                                             REAC 89 
      CALL RINDEX                                                       REAC 90 
C********* AT A PERIGEE                                                 REAC 91 
      PERIGE=S.EQ.(-UP)                                                 REAC 92 
C********* CORRECT MINOR ERRORS                                         REAC 93 
      IF (PERIGE) R(4)=0.                                               REAC 94 
C********* KEEP CONSISTENCY AFTER CORRECTING MINOR ERRORS               REAC 95 
      DRDT(1)=R(4)                                                      REAC 96 
C********* ON THE GROUND                                                REAC 97 
      GROUND=S.EQ.SG.AND.CROSSG                                         REAC 98 
C********* AT THE RECEIVER HEIGHT                                       REAC 99 
      THERE=S.EQ.SR.AND.CROSSR.AND..NOT.PERIGE                          REAC100 
C********* AT A CLOSEST APPROACH TO THE RECEIVER HEIGHT                 REAC101 
      MINDIS=PERIGE.AND.S.EQ.SR.AND.CROSSR                              REAC102 
      RSPACE=SPACE                                                      REAC103 
      V=SQRT(N2/(R(4)**2+R(5)**2+R(6)**2))                              REAC104 
C********* RENORMALIZE THE WAVE NORMAL DIRECTION TO = SQRT(REAL(N**2))  REAC105 
      R(4)=R(4)*V                                                       REAC106 
      R(5)=R(5)*V                                                       REAC107 
      R(6)=R(6)*V                                                       REAC108 
      RSTART=1.                                                         REAC109 
      IF (.NOT.SPACE) CALL PRINTR (8HENTR ION,0.)                       REAC110 
      RETURN                                                            REAC111 
         END                                                            REAC112-
