



      SUBROUTINE CLEAR(A,N)                                             CLEAR  2

C     SET N ELEMENTS OF ARRAY A TO ZERO                                 CLEAR  3

      REAL*8 A(N)                                                       CLEAR  4

      IF(N.LE.0) RETURN                                                 CLEAR  5

      DO    10 I=1,N                                                    CLEAR  6

10    A(I)=0.0                                                          CLEAR  7

      RETURN                                                            CLEAR  8

      END                                                               CLEAR  9











      FUNCTION ND2B(INDEC)                                              ND2B   2

C                                                                       ND2B   3

C     CONVERT A NUMERIC DECIMAL DIGIT STRING TO A BIT STRING            ND2B   4

C     WITH EACH BIT SET BY A CORRESPONDING DECIMAL DIGIT.               ND2B   5

C                                                                       ND2B   6

      ND2B=0                                                            ND2B   7

      IF(INDEC.LE.0) RETURN                                             ND2B   8

      M=INDEC                                                           ND2B   9

      MB=1                                                              ND2B  10

C                                                                       ND2B  11

10    IF(MOD(M,10).NE.0) ND2B=ND2B+MB                                   ND2B  12

      MB=MB*2                                                           ND2B  13

      M=M/10                                                            ND2B  14

      IF(M.GT.0) GO TO 10                                               ND2B  15

      END                                                               ND2B  16











      LOGICAL FUNCTION UCON(JC1,U,CONV)                                 UCON   2

c

      implicit real*8 (a-h,o-z)

c

      LOGICAL SETUCON                                                   UCON   3

C     PROVIDES KEYWORD UNITS CONVERSION FOR W-ARRAY INPUT.              UCON   4

C     UNITS SPECIFICATION MUST BE IN THE FORM 'DV UN' WHERE DV IS THE   UCON   5

C     TYPE OF VARIABLE(SUCH AS AN FOR ANGLE) AND UN IS THE UNITS CHOICE UCON   6

C     SUCH AS RD FOR RADIANS). RETURN VALUE IS THE CONVERSION FACTOR.   UCON   7

C     COMMON DECK "CUCON" INSERTED HERE                                 CUCON  2

      COMMON/UCONV/CNVV(4,4)                                            CUCON  4

      CHARACTER PCV*3,CNVC*2                                            CUCON  5

      COMMON/UCONC/PCV(4),CNVC(4,4)                                     CUCON  6

C     COMMON DECK "CONST" INSERTED HERE                                 CCONST 2

      COMMON/PCONST/CREF,RGAS,GAMMA                                     CCONST 4

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10                         CCONST 5

C     COMMON DECK "WWR" INSERTED HERE                                   CWWR   2

      PARAMETER (NWARSZ=1000)                                           CWW1   3

      real*8 id

      COMMON/WW/ID(10),MAXW,W(NWARSZ)                                   CWW1   4

      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON                              CWW2   2

      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),    CWW2   3

     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),      CWW2   4

     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),          CWW2   5

     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),         CWW2   6

     8 (RCVRH,W(20)),                                                   CWW2   7

     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))CWW2   8

     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),                      CWW2   9

     6 (HMIN,W(27)),(RGMAX,W(28)),                                      CWW2  10

     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),                     CWW2  11

     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),       CWW2  12

     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))        CWW2  13

     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),        CWW2  14

     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))              CWW2  15

     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))                   CWW2  16

      CHARACTER U*(*),PU*3,PV*2                                         UCON  11

C                                                                       UCON  12

      DATA NPV,NPU/3,4/                                                 UCON  13

      UCON=.FALSE.                                                      UCON  14

      CONV=1.0                                                          UCON  15

C                                                                       UCON  16

      READ(U,100) PU,PV                                                 UCON  17

100   FORMAT(A,A)                                                       UCON  18

C                                                                       UCON  19

      DO 2010 J=1,NPV                                                   UCON  20

2010     IF(PU.EQ.PCV(J)) GO TO 2015                                    UCON  21

      JC1=4                                                             UCON  22

      RETURN                                                            UCON  23

2015  JC1=J                                                             UCON  24

      DO 2020 K=1,NPU                                                   UCON  25

2020     IF(PV.EQ.CNVC(J,K)) GO TO 2025                                 UCON  26

      JC1=4                                                             UCON  27

      RETURN                                                            UCON  28

2025  UCON=.TRUE.                                                       UCON  29

      CONV=CNVV(JC1,K)                                                  UCON  30

      if(conv.eq.-1.0) conv = 1.0/earthr                                UCON305

      RETURN                                                            UCON  31

C                                                                       UCON  32

C     INITIAL CONVERSION CONSTANTS                                      UCON  33

      ENTRY SETUCON(JC1,U,CONV)                                         UCON  34

      SETUCON=.TRUE.                                                    UCON  35

C                                                                       UCON  36

      CNVV(1,2) = RAD                                                   UCON  37

      CNVV(1,3) = -1.0                                                  UCON380

      CNVV(2,2) = 1.0E-3                                                UCON  39

      CNVV(2,3) = 3.048006096E-4                                        UCON  40

      CNVV(2,4) = 1.852                                                 UCON  41

      CNVV(3,2) = PIT2                                                  UCON  42

      CNVV(3,3) = PIT2*1.0E3                                            UCON  43

      CNVV(3,4) = PIT2*RAD                                              UCON  44

      END                                                               UCON  45











      SUBROUTINE TRACE                                                  TRACE  3

c

      implicit real*8 (a-h,o-z)

c

C                                                 VERSION 3/4/86        TRACE  4

C     CALCULATES ONE RAYPATH FOR THE REQUESTED NUMBER OF HOPS           TRACE  5

C     (CROSSINGS OR CLOSEST APPROACH TO RECEIVER SURFACE; A CLOSEST     TRACE  6

C     APPROACH COUNTS AS TWO HOPS).                                     TRACE  7

C     THIS VERSION IS LIMITED TO USE WITH ATMOSPHERIC MODELS            TRACE  8

      REAL*8 NORYST                                                     TRACE  9

      PARAMETER (NORYST=0.0,MXPRGS=200)                                 TRACE 10

C                                                                       TRACE 11

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       TRACE 12

C     COMMON DECK "FLAG" INSERTED HERE                                  TRACE 13

      LOGICAL NEWWR,NEWWP,NEWTRC,PENET                                  TRACE 14

      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH      TRACE 15

C     COMMON DECK "TRAC" INSERTED HERE                                  CTRAC  2

      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY                    CTRAC  4

      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT     CTRAC  5

      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1             CTRAC  6

C     COMMON DECK "RR" INSERTED HERE                                    TRACE 17

      real*8 modrec

c     REAL MODREC                                                       TRACE 18

      COMMON/RR/ MODREC(4)                                              TRACE 19

      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH                                  TRACE 20

      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME            TRACE 21

C     COMMON DECK "GG" INSERTED HERE                                    TRACE 22

      real*8 modg

c     REAL MODG                                                         TRACE 23

      COMMON/GG/MODG(4)                                                 TRACE 24

      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH                                  TRACE 25

      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME            TRACE 26

C     COMMON DECK "RKTIME" INSERTED HERE                                TRACE 27

      COMMON/CRKTIME/RKTIME                                             TRACE 28

C     COMMON DECK "RINPLEX" INSERTED HERE                               TRACE 29

      REAL*8 KAY2,KAY2I                                                 TRACE 30

      COMPLEX*16 PNP,POLAR,LPOLAR                                       TRACE 31

      LOGICAL SPACE                                                     TRACE 32

c

      real*8 modrin   

      COMMON /RIN/ MODRIN(8),RAYNAME(2,3),TYPE(3),SPACE                 TRACE 33

      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I                             TRACE 34

      COMMON/RIN/PNP(10),POLAR,LPOLAR,SGN                               TRACE 35

      COMMON R(20),T,STP,DRDT(20)                                       TRACE 36

C     COMMON DECK "WWR" INSERTED HERE                                   TRACE 37

      PARAMETER (NWARSZ=1000)                                           TRACE 38

      real*8 id

      COMMON/WW/ID(10),MAXW,W(NWARSZ)                                   TRACE 39

      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON                              TRACE 40

      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),    TRACE 41

     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),      TRACE 42

     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),          TRACE 43

     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),         TRACE 44

     8 (RCVRH,W(20)),                                                   TRACE 45

     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))TRACE 46

     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),                      TRACE 47

     6 (HMIN,W(27)),(RGMAX,W(28)),                                      TRACE 48

     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),                     TRACE 49

     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),       TRACE 50

     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))        TRACE 51

     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),        TRACE 52

     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))              TRACE 53

     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))                   TRACE 54

C    THE FOLLOWING IS A LOCAL COMMON ONLY                               TRACE 55

C     COMMON DECK "TRLOCAL" INSERTED HERE                               TRACE 56

      COMMON/TRLOCAL/RSIGN,HOME,FDOT,GDOT,GOLD,GDOLD                    TRACE 57

      LOGICAL PCROSS,HOME,WASNT                                         TRACE 58

      LOGICAL LAUNCH,APOGEE,PLTENB,EXTON                                TRACE 59

      EQUIVALENCE (R(2),TH),(R(3),PHI)                                  TRACE 60

C                                                                       TRACE 61

      CHARACTER*10 XCOND,DEFCND                                         TRACE 62

      PARAMETER (DEFCND=' MAX HOPS')                                    TRACE 63

C                                                                       TRACE 64

      REAL*8 PRGHST(3,MXPRGS)                                           TRACE 65

      EQUIVALENCE(RKTIME,IRKTIME), (APOG,PRIGEE,W(80))                  TRACE 66

c

      real*8 gtoo(12),ftoo(12)

      equivalence (gtoo(1),g),(ftoo(1),f) 

C                                                                       TRACE 68

C     INCREMENT EVENT COUNT                                             TRACE 69

      rktime = rktime + 1.

c     IRKTIME=IRKTIME+1                                                 TRACE 70

      PLTENB=APOG.EQ.0.0                                                TRACE 71

C     ENABLE EXTINCTION TEST IF LIMIT IS SET AND VARIABLE IS BEING INTE TRACE 72

      EXTON=(EXTINC.GT.0.0) .AND. (W(58).NE.0.0)                        TRACE 73

C     POINT INDEX TO ABSORTION VARIABLE                                 TRACE 74

      NR=7                                                              TRACE 75

      IF(W(57).NE.0.0) NR=8                                             TRACE 76

      RMAX=HMAX+EARTHR                                                  TRACE 77

      RMIN=HMIN+EARTHR                                                  TRACE 78

      NHOP=HOP                                                          TRACE 79

      MAX=MAXSTP                                                        TRACE 80

      NSKIP=SKIP                                                        TRACE 81

C********* INITIALIZE PARAMETERS FOR INTEGRATION SUBROUTINE RKAM        TRACE 82

      MODE=INTYP                                                        TRACE 83

      STEP=STEP1                                                        TRACE 84

      E1MAX=MAXERR                                                      TRACE 85

      E1MIN=MAXERR/ERATIO                                               TRACE 86

      E2MAX=STPMAX                                                      TRACE 87

      E2MIN=STPMIN                                                      TRACE 88

      FACT=FACTR                                                        TRACE 89

      RSTART=1.                                                         TRACE 90

      CALL HAMLTN                                                       TRACE 91

      fdot = dot(ftoo)

c     FDOT=DOT(F)                                                       TRACE 92

C                                                                       TRACE 93

C     CHECK FOR EQUALITY WITH RECEIVER HEIGHT WITHIN MACHINE PRECISION  TRACE 94

      THERE=F.EQ.0.0                                                    TRACE 95

C                                                                       TRACE 96

      HOME=FDOT*F.GE.0.                                                 TRACE 97

C********* IHOP=0 TELLS PRINTR TO PRINT HEADING AND PUNCH A TRANSMITTER TRACE 98

C********* RAYSET AND NEWTRC=TRUE TELLS RAYPLT TO START A NEW RAY       TRACE 99

      NEWTRC=.TRUE.                                                     TRACE100

      IHOP=0                                                            TRACE101

C********* RESET PERIGEE PLOT COUNTER                                   TRACE102

      NPRGS=0                                                           TRACE103

C     USE CURRENT RELATIVE POSITION OF RAY TO PREDICT                   TRACE104

C     SIGN OF NEXT INTERSECTION OF RAY WITH RECEIVER HEIGHT             TRACE105

C     (WILL SOON BE REVERSED).                                          TRACE106

c

      xxx = 1.0

      RSIGN=DSIGN(xxx,F)                                                TRACE107

C     IF RAY IS LAUNCHED FROM RECEIVER HEIGHT, USE DIRECTION INSTEAD    TRACE108

      IF(THERE) RSIGN=DSIGN(xxx,FDOT)                                   TRACE109

      gdot = dot(gtoo)

c     GDOT=DOT(G)                                                       TRACE110

C     TEST CASE FOR RECEIVER ON THE GROUND                              TRACE111

      IF(F.EQ.G) RSIGN=1.0                                              TRACE112

C     IF RECEIVER IS ON TERRAIN THEN DIRECTION OF FIRST                 TRACE113

C     RECEIVER CROSSING IS DOWNWARD                                     TRACE114

      CALL PRINTR('TXMTR',RAYSET)                                       TRACE115

      IF (PLTENB) CALL RAYPLT                                           TRACE116

      NEWRAY=.TRUE.                                                     TRACE117

      STH0=DSIN(TH)                                                     TRACE118

      CTH0=DCOS(TH)                                                     TRACE119

      PHI0=PHI                                                          TRACE120

C                                                                       TRACE121

C     SET DEFAULT RAY EVENTS FOR 'PRINTR'                               TRACE122

      XCOND=DEFCND                                                      TRACE123

      XSET=NORYST                                                       TRACE124

C                                                                       TRACE125

C********* LOOP ON NUMBER OF HOPS                                       TRACE126

10    IHOP=IHOP+1                                                       TRACE127

C      REVERSE SIGN AT EACH CROSSING OF RECEIVER HEIGHT                 TRACE128

      RSIGN=-RSIGN                                                      TRACE129

      IF (IHOP.GT.NHOP) GO TO 100                                       TRACE130

      PENET=.FALSE.                                                     TRACE131

C********* LOOP ON MAXIMUM NUMBER OF STEPS PER HOP                      TRACE132

      DO 79 J=1,MAX                                                     TRACE133

C     LAUNCH=TRUE ONLY WHEN TRANSMITTER ON GROUND WITH DOWNWARD         TRACE134

C     TRANSMISSION(POSSIBLE ONLY ON 1ST STEP).                          TRACE135

      LAUNCH=G.EQ.0.0 .AND. GDOT.LT.0.0                                 TRACE136

C                                                                       TRACE137

C     SAVE CURRENT STATE VALUES FOR LATER COMPARISONS                   TRACE138

   12 DO 13 L=1,6                                                       TRACE139

      ROLD(L)=R(L)                                                      TRACE140

   13 DROLD(L)=DRDT(L)                                                  TRACE141

      FDOLD=FDOT                                                        TRACE142

      GDOLD=GDOT                                                        TRACE143

      FOLD=F                                                            TRACE144

      GOLD=G                                                            TRACE145

      TOLD=T                                                            TRACE146

C                                                                       TRACE147

C     PROCESS NEXT RAY POINT                                            TRACE148

      CALL RKAM                                                         TRACE149

      fdot = dot(ftoo)

c     FDOT=DOT(F)                                                       TRACE150

      WASNT=.NOT.HOME                                                   TRACE151

      HOME=FDOT*F.GE.0.                                                 TRACE152

C                                                                       TRACE153

C                                                                       TRACE154

C     LOOK FOR DOWNGOING CROSSING OF RECEIVER SURFACE                   TRACE155

      IF(.NOT.LAUNCH.AND.F.LT.0.0.AND.RSIGN.LT.0.0) GO TO 50            TRACE156

C     CHECK FOR CASE OF CLOSEST APPROACH                                TRACE157

      IF (HOME.AND.WASNT) GO TO 30                                      TRACE158

C     DETERMINE IF A GROUND CROSSING HAS OCCURED                        TRACE159

      IF(PCROSS(gtoo,GDOT)) GO TO 20                                    TRACE160

C     LOOK FOR UPGOING CROSSING OF RECEIVER SURFACE                     TRACE161

      IF(.NOT.LAUNCH.AND.F.GT.0.0.AND.RSIGN.GT.0.0) GO TO 50            TRACE162

      GO TO 18                                                          TRACE163

C                                                                       TRACE164

C     CHECK FOR PERIGEE CONDITION                                       TRACE165

18    IF (DROLD(1).GE.0..OR.DRDT(1).LE.0.) GO TO 25                     TRACE166

      CALL PRINTR(' PERIGEE ',NORYST)                                   TRACE167

      IF(PRIGEE.EQ.0.0) GO TO 25                                        TRACE168

C                                                                       TRACE169

      IF(NPRGS.GE.MXPRGS) CALL STOPIT('PRG LIMIT')                      TRACE170

      NPRGS=NPRGS+1                                                     TRACE171

      CALL RMOVE(PRGHST(1,NPRGS),R,3)                                   TRACE172

C                                                                       TRACE173

25    APOGEE=DROLD(1).GT.0..AND.DRDT(1).LT.0.                           TRACE174

      PLTENB=APOG.EQ.0.0.OR.APOGEE                                      TRACE175

C                                                                       TRACE176

      IF(APOGEE) CALL PRINTR(' APOGEE  ',NORYST)                        TRACE177

      IF (DROLD(2)*DRDT(2).LT.0.) CALL PRINTR(' MAX LAT ',NORYST)       TRACE178

      IF (DROLD(3)*DRDT(3).LT.0.) CALL PRINTR(' MAX LONG',NORYST)       TRACE179

      DO 14 I=4,6                                                       TRACE180

      IF (ROLD(I)*R(I).LT.0.) CALL PRINTR(' WAVE REV',NORYST)           TRACE181

   14 CONTINUE                                                          TRACE182

      GO TO 75                                                          TRACE183

C********* RAY WENT UNDERGROUND                                         TRACE184

C     USE 'FULL' REGULAR PARABOLIC FIT(ENTRY 'FIT')                     TRACE185

20    continue

      avl = -1.0

      CALL RCROSS(avl,gtoo,GDOT,'GGRND REF',GROUND)                     TRACE186

      GO TO 75                                                          TRACE187

C********* RAY MAY HAVE MADE A CLOSEST APPROACH                         TRACE188

C     USE 'FULL' REGULAR PARABOLIC FIT                                  TRACE189

30    CALL FIT(ftoo,FOLD,FDOLD)                                         TRACE190

      IF(RAD.GE.0.0) GO TO 501                                          TRACE191

C     ESTIMATE TIME OF CLOSEST APPROACH                                 TRACE192

      TP=T-ZDOT/D2Z                                                     TRACE193

      CALL GRAZE(ftoo,RSIGN,TP)                                         TRACE194

      FDOT=ZDOT                                                         TRACE195

      IF (THERE) GO TO 51                                               TRACE196

C                                                                       TRACE197

C     SET DRDT(1)=0 TO AVOID INCORRECT APOGEE OR PERIGEE PRINTOUT, A NEWTRACE198

C     VALUE WILL BE COMPUTED BEFORE FURTHER ANALYSIS IS DONE            TRACE199

40    DRDT(1)=0.                                                        TRACE200

      HPUNCH=R(1)-EARTHR                                                TRACE201

      CALL PRINTR('MMIN DIST',RAYSET)                                   TRACE202

      IF(PLTENB) CALL RAYPLT                                            TRACE203

      IF (IHOP.GE.NHOP) GO TO 100                                       TRACE204

      IHOP=IHOP+1                                                       TRACE205

      RSIGN=-RSIGN                                                      TRACE206

      CALL PRINTR ('MMIN DIST',RAYSET)                                  TRACE207

      GO TO 89                                                          TRACE208

C********* RAY CROSSED RECEIVER HEIGHT                                  TRACE209

50    CALL FIT(F,FOLD,FDOLD)                                            TRACE210

C     ESTIMATE GROUP TIME WHEN RAY CROSSES RECEIVER HEIGHT              TRACE211

C     (IN THE CORRECT DIRECTION).                                       TRACE212

501   TC=T-2.0*F/(ZDOT+DSIGN(RAD1,RSIGN))                               TRACE213

      CALL BACKUP(ftoo,RSIGN,TC)                                        TRACE214

      FDOT=ZDOT                                                         TRACE215

      IF(.NOT.THERE) GO TO 40                                           TRACE216

51    HPUNCH=R(1)-EARTHR                                                TRACE217

      CALL PRINTR('RRCVR    ',RAYSET)                                   TRACE218

      IF (PLTENB) CALL RAYPLT                                           TRACE219

      IF (GET(ftoo).NE.GET(gtoo)) GO TO 89                              TRACE220

C     RECEIVER IS ON TERRAINE                                           TRACE221

      IF (IHOP.GE.NHOP) GO TO 100                                       TRACE222

      IHOP=IHOP+1                                                       TRACE223

      RSIGN=-RSIGN                                                      TRACE224

C********* GROUND REFLECT                                               TRACE225

      reefer = reflect(gtoo)

c     CALL REFLECT(G)                                                   TRACE226

      CALL HAMLTN                                                       TRACE227

      FDOT=DOT(ftoo)                                                    TRACE228

      RSTART=1.                                                         TRACE229

      HPUNCH=R(1)-EARTHR                                                TRACE230

      CALL PRINTR('GGRND REF',RAYSET)                                   TRACE231

      THERE=.TRUE.                                                      TRACE232

      HPUNCH=R(1)-EARTHR                                                TRACE233

      CALL PRINTR ('RRCVR    ',RAYSET)                                  TRACE234

      GO TO 89                                                          TRACE235

C*********                                                              TRACE236

75    IF (PLTENB) CALL RAYPLT                                           TRACE237

      PLTENB=APOG.EQ.0.0                                                TRACE238

      IF(EXTON) THEN                                                    TRACE239

         IF (R(NR).GE.EXTINC) XCOND='EEXTINC'                           TRACE240

      ENDIF                                                             TRACE241

      IF (R(1).GT.RMAX.AND.F.GT.0.0.AND.FDOT.GT.0.) XCOND='UMAX HT'     TRACE242

      IF (R(1).LT.RMIN.AND.F.LT.0.0.AND.FDOT.LT.0.) XCOND='DMIN HT'     TRACE243

      RANGE=EARTHR*DACOS(DCOS(TH)*CTH0+DSIN(TH)*STH0*DCOS(PHI-PHI0))    TRACE244

      IF(RGMAX.GT.0.0.AND.RGMAX.LT.RANGE) XCOND='FMAX RANGE'            TRACE245

      IF(XCOND.NE.DEFCND) GO TO 90                                      TRACE246

C                                                                       TRACE247

      IF (MOD(J,NSKIP).EQ.0) CALL PRINTR('        ',NORYST)             TRACE248

   79 CONTINUE                                                          TRACE249

C********* EXCEEDED MAXIMUM NUMBER OF STEPS                             TRACE250

      HPUNCH=R(1)-EARTHR                                                TRACE251

      CALL PRINTR('SSTEP MAX',RAYSET)                                   TRACE252

      GO TO 100                                                         TRACE253

C*********                                                              TRACE254

   89 HOME=.TRUE.                                                       TRACE255

      GDOT=DOT(gtoo)                                                    TRACE256

      GO TO 10                                                          TRACE257

C********* RAY PENETRATED COMPUTATIONAL AREA BOUNDARY AND WAS HEADING OUTRACE258

90    PENET=.TRUE.                                                      TRACE259

      HPUNCH=R(1)-EARTHR                                                TRACE260

      XSET=RAYSET                                                       TRACE261

C                                                                       TRACE262

C     NORMAL EXIT                                                       TRACE263

100   CALL PRINTR(XCOND,XSET)                                           TRACE264

      IF(NPRGS.LE.0) RETURN                                             TRACE265

C                                                                       TRACE266

      CALL DASH                                                         TRACE267

      NEWTRC=.TRUE.                                                     TRACE268

      DO 110 I=1,NPRGS                                                  TRACE269

      CALL RMOVE(R,PRGHST(1,I),3)                                       TRACE270

110   CALL RAYPLT                                                       TRACE271

C                                                                       TRACE272

      CALL RESET('DASH')                                                TRACE273

      END                                                               TRACE274











      LOGICAL FUNCTION PCROSS(Z,ZDT)                                    TRACE275

c

      implicit real*8 (a-h,o-z)

c

C     RETURNS 'TRUE' IF PARABOLIC 'FIT2'(SEE 'FIT' ROUTINE)             TRACE276

C     INDICATES CROSSING OF SURFACE 'Z'.                                TRACE277

      real*8 z(12)

      REAL*8 ZDT(3)                                                     TRACE278

C     COMMON DECK "TRAC" INSERTED HERE                                  CTRAC  2

      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY                    CTRAC  4

      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT     CTRAC  5

      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1             CTRAC  6

C                                                                       TRACE280

C     USE SPECIAL PARABOLIC FIT GIVING 'SMT' AND 'OSMT'                 TRACE281

      ZDT(1)=DOT(Z)                                                     TRACE282

      PCROSS=.TRUE.                                                     TRACE283

      IF(Z(1).LT.0.0) GO TO 20                                          TRACE284

      IF(ZDT(1).LE.0.0.OR.ZDT(3).GE.0.0) GO TO 15                       TRACE285

      CALL FIT2(Z,ZDT(2),ZDT(3))                                        TRACE286

      IF(SMT.GT.Z(1).OR.OSMT.GT.ZDT(2)) GO TO 20                        TRACE287

C                                                                       TRACE288

15    PCROSS=.FALSE.                                                    TRACE289

20    RETURN                                                            TRACE290

      END                                                               TRACE291











      SUBROUTINE RCROSS(S,Z,ZDT,EVENT,QMODE)                            TRACE292

c

      implicit real*8 (a-h,o-z)

      real*8 z(12)

c

C     FIND ESTIMATED CROSSING POINT OF SURFACE 'Z' THEN USE             TRACE293

C     ROUTINE 'BACKUP' TO GO THERE AND PERFORM A RAY REFLECTION         TRACE294

C     WITH ROUTINE 'REFLECT'.                                           TRACE295

C     EXTENDS RAY PROPAGATION TO INTERSECTION WITH REFLECTING SURFACE   TRACE296

C     Z AND CALLS 'REFLECT' TO OBTAIN VECTOR COMPONENTS.                TRACE297

C                                                                       TRACE298

      CHARACTER EVENT*8                                                 TRACE299

      REAL*8 ZDT(3)                                                     TRACE300

      LOGICAL QMODE                                                     TRACE301

C     COMMON DECK "TRAC" INSERTED HERE                                  CTRAC  2

      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY                    CTRAC  4

      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT     CTRAC  5

      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1             CTRAC  6

C     COMMON DECK "RR" INSERTED HERE                                    TRACE303

      real*8 modrec

c     REAL MODREC                                                       TRACE304

      COMMON/RR/ MODREC(4)                                              TRACE305

      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH                                  TRACE306

      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME            TRACE307

C     COMMON DECK "WWR" INSERTED HERE                                   TRACE308

      PARAMETER (NWARSZ=1000)                                           TRACE309

      real*8 id

      COMMON/WW/ID(10),MAXW,W(NWARSZ)                                   TRACE310

      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON                              TRACE311

      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),    TRACE312

     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),      TRACE313

     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),          TRACE314

     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),         TRACE315

     8 (RCVRH,W(20)),                                                   TRACE316

     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))TRACE317

     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),                      TRACE318

     6 (HMIN,W(27)),(RGMAX,W(28)),                                      TRACE319

     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),                     TRACE320

     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),       TRACE321

     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))        TRACE322

     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),        TRACE323

     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))              TRACE324

     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))                   TRACE325

C     COMMON DECK "FLAG" INSERTED HERE                                  TRACE326

      LOGICAL NEWWR,NEWWP,NEWTRC,PENET                                  TRACE327

      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH      TRACE328

      COMMON R(20),T,STP,DRDT(20)                                       TRACE329

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       TRACE330

      LOGICAL HOME                                                      TRACE331

C     COMMON DECK "TRLOCAL" INSERTED HERE                               TRACE332

      COMMON/TRLOCAL/RSIGN,HOME,FDOT,GDOT,GOLD,GDOLD                    TRACE333

c

      real*8 ftoo(12)

      equivalence (ftoo(1),f)

C                                                                       TRACE334

      IF(QMODE) GO TO 60                                                TRACE335

      CALL FIT(Z,ZDT(2),ZDT(3))                                         TRACE336

      TC=T-2.0*Z(1)/(ZDOT-RAD1)                                         TRACE337

      avl = -1.0

      CALL BACKUP(Z,avl,TC)                                             TRACE338

C********* GROUND REFLECT                                               TRACE339

60    reefer = reflect(z)

c60   CALL REFLECT(Z)                                                   TRACE340

      CALL HAMLTN                                                       TRACE341

      ZDT(1)=DOT(Z)                                                     TRACE342

      FDOT=DOT(ftoo)                                                    TRACE343

      HOME=FDOT*F.GE.0.0                                                TRACE344

      RSTART=1.                                                         TRACE345

      HPUNCH=R(1)-EARTHR                                                TRACE346

      CALL PRINTR(EVENT,RAYSET)                                         TRACE347

      IF(F.NE.Z(1)) RETURN                                              TRACE348

C     AVOID RECEIVER CROSSING AT THE TRANSMITTER IF RECEIVER ON TERRAIN TRACE349

      RSIGN=S                                                           TRACE350

      HOME=.TRUE.                                                       TRACE351

      END                                                               TRACE352











      SUBROUTINE HAMLTN                                                 HAMLTN 2

c

      implicit real*8 (a-h,o-z)

c

C********* CALCULATES HAMILTONS EQUATIONS FOR RAY TRACING               HAMLTN 3

C     AND OTHER QUANTITIES TO BE INTEGRATED                             HAMLTN 4

C     COMMON DECK "GG" INSERTED HERE                                    CGG    2

      real*8 modg

c     REAL MODG                                                         CGG    4

      COMMON/GG/MODG(4)                                                 CGG    5

      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH                                  CGG    6

      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME            CGG    7

C     COMMON DECK "CONST" INSERTED HERE                                 CCONST 2

      COMMON/PCONST/CREF,RGAS,GAMMA                                     CCONST 4

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10                         CCONST 5

C     COMMON DECK "RINREAL" INSERTED HERE                               CRINREA2

      LOGICAL SPACE                                                     CRINREA4

      REAL*8 LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I                        CRINREA5

      CHARACTER DISPM*6                                                 CRINREA6

      COMMON/RINPL/DISPM                                                CRINREA7

      real*8 modrin

      COMMON /RIN/ MODRIN(8),RAYNAME(2,3),TYPE(3),SPACE                 CRINREA8

      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I,                            CRINREA9

     1              H,HI,PHT,PHTI,PHR,PHRI,PHTH,PHTHI,PHPH,PHPHI        CRINRE10

     2, PHOW,PHOWI,PHKR,PHKRI,PHKTH,PHKTI, PHKPH,PHKPI                  CRINRE11

     3 ,KPHK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN                       CRINRE12

      COMMON R(20),T,STP,DRDT(20)                                       HAMLTN 8

      PARAMETER (NWARSZ=1000)                                           CWW1   3

      real*8 id

      COMMON/WW/ID(10),MAXW,W(NWARSZ)                                   CWW1   4

      EQUIVALENCE (TH,R(2)),(PH,R(3)),(KR,R(4)),(KTH,R(5)),(KPH,R(6)),  HAMLTN10

     1 (DTHDT,DRDT(2)),(DPHDT,DRDT(3)),(DKRDT,DRDT(4)),(DKTHDT,DRDT(5)),HAMLTN11

     2 (DKPHDT,DRDT(6)),(HMAX,W(26)),(OW,W(6))                          HAMLTN12

      REAL*8 KR,KTH,KPH                                                 HAMLTN13

      STH=DSIN(TH)                                                      HAMLTN14

      CTH=DSIN(PID2-TH)                                                 HAMLTN15

      RSTH=R(1)*STH                                                     HAMLTN16

      RCTH=R(1)*CTH                                                     HAMLTN17

      CALL DISPER                                                       HAMLTN18

      DENPHC=1.0/(PHOW*CREF)                                            HAMLTN19

      DRDT(1)=-PHKR*DENPHC                                              HAMLTN20

      DTHDT=-PHKTH*DENPHC/R(1)                                          HAMLTN21

      DPHDT=-PHKPH*DENPHC/RSTH                                          HAMLTN22

      DKRDT=PHR*DENPHC+KTH*DTHDT+KPH*STH*DPHDT                          HAMLTN23

      DKTHDT=(PHTH*DENPHC-KTH*DRDT(1)+KPH*RCTH*DPHDT)/R(1)              HAMLTN24

      DKPHDT=(PHPH*DENPHC-KPH*STH*DRDT(1)-KPH*RCTH*DTHDT)/RSTH          HAMLTN25

      NR=6                                                              HAMLTN26

C********* PHASE PATH                                                   HAMLTN27

      IF (W(57).EQ.0.) GO TO 10                                         HAMLTN28

      NR=NR+1                                                           HAMLTN29

      DRDT(NR)=-     KPHK/PHOW/OW                                       HAMLTN30

C********* ABSORPTION                                                   HAMLTN31

   10 IF (W(58).EQ.0.) GO TO 15                                         HAMLTN32

      NR=NR+1                                                           HAMLTN33

      DRDT(NR)= 10./ALN10*KPHK*KAY2I/(KR*KR+KTH*KTH+KPH*KPH)*DENPHC     HAMLTN34

C********* DOPPLER SHIFT                                                HAMLTN35

   15 IF (W(59).EQ.0.) GO TO 20                                         HAMLTN36

      NR=NR+1                                                           HAMLTN37

      DRDT(NR)=-PHT*DENPHC/PIT2                                         HAMLTN38

C********* GEOMETRICAL PATH LENGTH                                      HAMLTN39

   20 IF (W(60).EQ.0.) GO TO 25                                         HAMLTN40

      NR=NR+1                                                           HAMLTN41

      DRDT(NR)=DSQRT(PHKR**2+PHKTH**2+PHKPH**2)*DABS(DENPHC)            HAMLTN42

C********* TERRAIN FUNCTION AS COMPLETE INTEGRAL                        HAMLTN43

25    IF(DABS(W(61))+DABS(W(62))+DABS(W(63))+DABS(W(64)).EQ.0.0)

     1 GO TO 45                                                         HAMLTN44

      CALL TOPOG                                                        HAMLTN45

C                                                                       HAMLTN46

      IF(W(61).EQ.0.) GO TO 30                                          HAMLTN47

      NR=NR+1                                                           HAMLTN48

      DRDT(NR)=PGR*DRDT(1)+PGTH*DTHDT+PGPH*DPHDT                        HAMLTN49

C********* TERRAIN TIME DERIVATIVES AS COMPLETE INTEGRALS               HAMLTN50

   30 IF(W(62).EQ.0.) GO TO 35                                          HAMLTN51

      NR=NR+1                                                           HAMLTN52

C     DERIVATIVE WITH RESPECT TO R                                      HAMLTN53

      DRDT(NR)=PGRR*DRDT(1)+PGRTH*DTHDT+PGRPH*DPHDT                     HAMLTN54

   35 IF(W(63).EQ.0.) GO TO 40                                          HAMLTN55

      NR=NR+1                                                           HAMLTN56

C     DERIVATIVE WITH RESPECT TO THETA                                  HAMLTN57

      DRDT(NR)=PGRTH*DRDT(1)+PGTHTH*DTHDT+PGTHPH*DPHDT                  HAMLTN58

   40 IF(W(64).EQ.0.) GO TO 45                                          HAMLTN59

      NR=NR+1                                                           HAMLTN60

C     DERIVATIVE WITH RESPECT TO PHI                                    HAMLTN61

      DRDT(NR)=PGRPH*DRDT(1)+PGTHPH*DTHDT+PGPHPH*DPHDT                  HAMLTN62

C********* OTHER CALCULATIONS                                           HAMLTN63

45    CONTINUE                                                          HAMLTN64

      RETURN                                                            HAMLTN65

         END                                                            HAMLTN66











      SUBROUTINE RKAM                                                   RKAM   2

c

      implicit real*8 (a-h,o-z)

c

C     KEEPS TRACK OF INTERNAL INTEGRATION STEPS PERFORMED BY THE RKAM1  RKAM   3

C     ROUTINE AND MAKES THEM AVAILABLE TO CALLING ROUTINES.             RKAM   4

      COMMON/RKAMS/XV(5),FV(4,20),YU(5,20),EPM,ALPHA,MM                 RKAM   5

C     COMMON DECK "RK" INSERTED HERE                                    CRK    2

C     DEFINE SIZE REQUIRED FOR RAY STATE SAVE ARRAY                     CRK    4

      PARAMETER (LRKAMS=87+2*100,NXRKMS=12+LRKAMS,MXEQPT=21)            CRK    5

      PARAMETER (NRKSAV=NXRKMS+MXEQPT-1)                                CRK    6

      COMMON /RK/ NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART    CRK    7

C     COMMON DECK "RKTIME" INSERTED HERE                                CRKTIME2

      COMMON/CRKTIME/RKTIME                                             CRKTIME4

      COMMON Y(20),T,SPACE,DYDT(20)                                     RKAM   8

      REAL*8 RV(1)                                                      RKAM   9

      DOUBLE PRECISION YU                                               RKAM  10

      EQUIVALENCE(RKTIME,IRKTIME),(RV,R)                                RKAM  11

C                                                                       RKAM  12

      REAL*8 SVBUF(NRKSAV)                                              RKAM  13

      LOGICAL SAVED                                                     RKAM  14

C                                                                       RKAM  15

      DATA SAVED/.FALSE./                                               RKAM  16

C     PERFORM CLOCK ADVANCE(SEE 'GET' ROUTINE)                          RKAM  17

      rktime = rktime + 1.0

c     IRKTIME=IRKTIME+1                                                 RKAM  18

C                                                                       RKAM  19

      IF(RSTART.NE.0.0 .OR. MODE.LE.2) GO TO 250                        RKAM  20

C                                                                       RKAM  21

      NV=NV+1                                                           RKAM  22

      IF(NV.LE.4) GO TO 300                                             RKAM  23

C                                                                       RKAM  24

250   TOLD=T                                                            RKAM  25

C     PERFORM NUMERICAL INTEGRATION TO TIME 'T'                         RKAM  26

      CALL RKAM1                                                        RKAM  27

      IF(MODE.LE.2) GO TO 400                                           RKAM  28

      IF(RSTART.NE.0.0) GO TO 260                                       RKAM  29

C     SEARCH FOR REQUESTED VALUE IN PIPELINE                            RKAM  30

C     TAKING SIGN OF 'STEP' INTO ACCOUNT.                               RKAM  31

      TUP=TOLD+.25*STEP                                                 RKAM  32

      DO 150 NV=1,3                                                     RKAM  33

      IF(SPACE.GT.0.0.AND.TUP.LT.XV(NV)) GO TO 300                      RKAM  34

      IF(SPACE.LT.0.0.AND.TUP.GT.XV(NV)) GO TO 300                      RKAM  35

150   CONTINUE                                                          RKAM  36

C     GOT IT ALREADY, RETURN                                            RKAM  37

      GO TO 400                                                         RKAM  38

C                                                                       RKAM  39

260   NV=1                                                              RKAM  40

C     RETRIEVE REQUESTED VALUES FROM PIPELINE                           RKAM  41

300   T=XV(NV)                                                          RKAM  42

      DO 350 I=1,NEQS                                                   RKAM  43

      Y(I)=YU(NV,I)                                                     RKAM  44

350   DYDT(I)=FV(NV,I)                                                  RKAM  45

C                                                                       RKAM  46

C     STANDARD EXIT SEQUENCE                                            RKAM  47

400   RETURN                                                            RKAM  48

C                                                                       RKAM  49

      ENTRY RKSAVE(SVBUF)                                               RKAM  50

C                                                                       RKAM  51

      SVBUF(1)=NV                                                       RKAM  52

      SVBUF(2)=NEQS                                                     RKAM  53

      SVBUF(3)=MM                                                       RKAM  54

      SVBUF(4)=SPACE                                                    RKAM  55

      SVBUF(5)=MODE                                                     RKAM  56

      CALL RMOVE(SVBUF(6),E1MAX,6)                                      RKAM  57

      CALL RMOVE(SVBUF(12),XV,LRKAMS)                                   RKAM  58

      CALL RMOVE(SVBUF(NXRKMS),Y,MXEQPT)                                RKAM  59

C                                                                       RKAM  60

      SAVED=.TRUE.                                                      RKAM  61

      RETURN                                                            RKAM  62

      ENTRY RKRSTR(SVBUF)                                               RKAM  63

C                                                                       RKAM  64

      IF(.NOT.SAVED) RETURN                                             RKAM  65

C                                                                       RKAM  66

      NV=SVBUF(1)                                                       RKAM  67

      NEQS=SVBUF(2)                                                     RKAM  68

      MM=SVBUF(3)                                                       RKAM  69

      SPACE=SVBUF(4)                                                    RKAM  70

      MODE=SVBUF(5)                                                     RKAM  71

      CALL RMOVE(E1MAX,SVBUF(6),6)                                      RKAM  72

      CALL RMOVE(XV,SVBUF(12),LRKAMS)                                   RKAM  73

      CALL RMOVE(Y,SVBUF(NXRKMS),MXEQPT)                                RKAM  74

C                                                                       RKAM  75

      WRITE(3,*) 'RESTORING NV,NEQS,MM,ALPHA,T='                        RKAM  76

     1   ,NV,NEQS,MM,ALPHA,T                                            RKAM  77

      RETURN                                                            RKAM  78

      END                                                               RKAM  79











      SUBROUTINE RKAM1                                                  RKAM1  2

c

      implicit real*8 (a-h,o-z)

c 

C        NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS                RKAM1  3

C     THIS ROUTINE IS A MODIFICATION OF RKAMSUB, WHICH WAS WRITTEN      RKAM1  4

C     BY G.J. LASTMAN AND IS AVAILABLE THROUGH THE CDC CO-OP LIBRARY    RKAM1  5

C     AS 'D2 UTEX RKAMSUB'.                                             RKAM1  6

      COMMON /RK/ NN,SPACE,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART     RKAM1  7

      COMMON Y(20),T,STEP,DYDT(20)                                      RKAM1  8

      COMMON/RKAMS/XV(5),FV(4,20),YU(5,20),EPM,ALPHA,MM                 RKAM1  9

      DIMENSION DELY(4,20),BET(4)                                       RKAM1 10

      DOUBLE PRECISION YU                                               RKAM1 11

C                                                                       RKAM1 12

      IF (RSTART.EQ.0.) GO TO 1000                                      RKAM1 13

      LL=1                                                              RKAM1 14

      MM=1                                                              RKAM1 15

      IF (MODE.EQ.1) MM=4                                               RKAM1 16

      ALPHA=T                                                           RKAM1 17

      EPM=0.0                                                           RKAM1 18

      BET(1)=0.5                                                        RKAM1 19

      BET(2)=0.5                                                        RKAM1 20

      BET(3)=1.0                                                        RKAM1 21

      BET(4)=0.0                                                        RKAM1 22

      STEP=SPACE                                                        RKAM1 23

      R=19.0/270.0                                                      RKAM1 24

      XV(MM)=T                                                          RKAM1 25

      IF (E1MIN.LE.0.) E1MIN=E1MAX/50.                                  RKAM1 26

      IF (FACT.LE.0.) FACT=0.5                                          RKAM1 27

      CALL HAMLTN                                                       RKAM1 28

      DO 320 I=1,NN                                                     RKAM1 29

      FV(MM,I)=DYDT(I)                                                  RKAM1 30

  320 YU(MM,I)=Y(I)                                                     RKAM1 31

      RSTART=0.                                                         RKAM1 32

      GO TO 1001                                                        RKAM1 33

 1000 IF (MODE.NE.1) GO TO 2000                                         RKAM1 34

C                                                                       RKAM1 35

C          RUNGE-KUTTA                                                  RKAM1 36

 1001 DO 1034 K=1,4                                                     RKAM1 37

      DO 1350 I=1,NN                                                    RKAM1 38

      DELY(K,I)=STEP*FV(MM,I)                                           RKAM1 39

      Z=YU(MM,I)                                                        RKAM1 40

 1350 Y(I)=Z+BET(K)*DELY(K,I)                                           RKAM1 41

      T=BET(K)*STEP+XV(MM)                                              RKAM1 42

      CALL HAMLTN                                                       RKAM1 43

      DO 1034 I=1,NN                                                    RKAM1 44

 1034 FV(MM,I)=DYDT(I)                                                  RKAM1 45

      DO 1039 I=1,NN                                                    RKAM1 46

      DEL=(DELY(1,I)+2.0*DELY(2,I)+2.0*DELY(3,I)+DELY(4,I))/6.0         RKAM1 47

 1039 YU(MM+1,I)=YU(MM,I)+DEL                                           RKAM1 48

      MM=MM+1                                                           RKAM1 49

      XV(MM)=XV(MM-1)+STEP                                              RKAM1 50

      DO 1400 I=1,NN                                                    RKAM1 51

 1400 Y(I)=YU(MM,I)                                                     RKAM1 52

      T=XV(MM)                                                          RKAM1 53

      CALL HAMLTN                                                       RKAM1 54

      IF (MODE.EQ.1) GO TO 42                                           RKAM1 55

      DO 150 I=1,NN                                                     RKAM1 56

  150 FV(MM,I)=DYDT(I)                                                  RKAM1 57

      IF (MM.LE.3) GO TO 1001                                           RKAM1 58

C                                                                       RKAM1 59

C          ADAMS-MOULTON                                                RKAM1 60

 2000 DO 2048 I=1,NN                                                    RKAM1 61

      DEL=STEP*(55.*FV(4,I)-59.*FV(3,I)+37.*FV(2,I)-9.*FV(1,I))/24.     RKAM1 62

      Y(I)=YU(4,I)+DEL                                                  RKAM1 63

 2048 DELY(1,I)=Y(I)                                                    RKAM1 64

      T=XV(4)+STEP                                                      RKAM1 65

      CALL HAMLTN                                                       RKAM1 66

      XV(5)=T                                                           RKAM1 67

      DO 2051 I=1,NN                                                    RKAM1 68

      DEL=STEP*(9.*DYDT(I)+19.*FV(4,I)-5.*FV(3,I)+FV(2,I))/24.          RKAM1 69

      YU(5,I)=YU(4,I)+DEL                                               RKAM1 70

 2051 Y(I)=YU(5,I)                                                      RKAM1 71

      CALL HAMLTN                                                       RKAM1 72

      IF (MODE.LE.2) GO TO 42                                           RKAM1 73

C                                                                       RKAM1 74

C          ERROR ANALYSIS                                               RKAM1 75

      SSE=0.0                                                           RKAM1 76

      DO 3033 I=1,NN                                                    RKAM1 77

      EPSIL=R*DABS(Y(I)-DELY(1,I))                                      RKAM1 78

      IF (MODE.EQ.3.AND.Y(I).NE.0.) EPSIL=EPSIL/DABS(Y(I))              RKAM1 79

      IF (SSE.LT.EPSIL) SSE=EPSIL                                       RKAM1 80

 3033 CONTINUE                                                          RKAM1 81

      IF (E1MAX.GT.SSE) GO TO 3035                                      RKAM1 82

      IF (DABS(STEP).LE.E2MIN) GO TO 42                                 RKAM1 83

      LL=1                                                              RKAM1 84

      MM=1                                                              RKAM1 85

      STEP=STEP*FACT                                                    RKAM1 86

      GO TO 1001                                                        RKAM1 87

 3035 IF (LL.LE.1.OR.SSE.GE.E1MIN.OR.E2MAX.LE.DABS(STEP)) GO TO 42      RKAM1 88

      LL=2                                                              RKAM1 89

      MM=3                                                              RKAM1 90

      XV(2)=XV(3)                                                       RKAM1 91

      XV(3)=XV(5)                                                       RKAM1 92

      DO 5363 I=1,NN                                                    RKAM1 93

      FV(2,I)=FV(3,I)                                                   RKAM1 94

      FV(3,I)=DYDT(I)                                                   RKAM1 95

      YU(2,I)=YU(3,I)                                                   RKAM1 96

 5363 YU(3,I)=YU(5,I)                                                   RKAM1 97

      STEP=2.0*STEP                                                     RKAM1 98

      GO TO 1001                                                        RKAM1 99

C                                                                       RKAM1100

C          EXIT ROUTINE                                                 RKAM1101

   42 LL=2                                                              RKAM1102

      MM=4                                                              RKAM1103

      DO 12 K=1,3                                                       RKAM1104

      XV(K)=XV(K+1)                                                     RKAM1105

      DO 12 I=1,NN                                                      RKAM1106

      FV(K,I)=FV(K+1,I)                                                 RKAM1107

   12 YU(K,I)=YU(K+1,I)                                                 RKAM1108

      XV(4)=XV(5)                                                       RKAM1109

      DO 52 I=1,NN                                                      RKAM1110

      FV(4,I)=DYDT(I)                                                   RKAM1111

   52 YU(4,I)=YU(5,I)                                                   RKAM1112

      IF (MODE.LE.2) RETURN                                             RKAM1113

      E=DABS(XV(4)-ALPHA)                                               RKAM1114

      IF (E.LE.EPM+.25*STEP) GO TO 2000                                 RKAM1115

      EPM=E                                                             RKAM1116

      RETURN                                                            RKAM1117

         END                                                            RKAM1118











      SUBROUTINE BACKUP(Z,RSIGN,TC)                                     BACKUP 2

c

      implicit real*8 (a-h,o-z)

      real*8 z(12)

      real*8 kparlel

c

C     MOVES THE RAY TO THE CLOSEST INTERSECTION WITH THE RECEIVER       BACKUP 3

C     OR TERRAIN SURFACE(VARIABLES 'FSELECT' OR 'GSELECT' IN LABELED    BACKUP 4

C     COMMONS /RR/ OR /GG/, RESPECTIVELY, TELL WHICH KIND OF SURFACE).  BACKUP 5

C                                                                       BACKUP 6

      CHARACTER*9 NBAK,NGRAZ,NTRY                                       BACKUP 7

C                                                                       BACKUP 8

      PARAMETER (PRNZTL=0.5E-4,PRNDZTL=1.E-6)                           BACKUP 9

C                                                                       BACKUP10

      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART       BACKUP11

C     COMMON DECK "TRAC" INSERTED HERE                                  CTRAC  2

      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY                    CTRAC  4

      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT     CTRAC  5

      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1             CTRAC  6

      COMMON R(20),T,STP,DRDT(20)                                       BACKUP13

      PARAMETER (NWARSZ=1000)                                           CWW1   3

      real*8 id

      COMMON/WW/ID(10),MAXW,W(NWARSZ)                                   CWW1   4

      EQUIVALENCE (EARTHR,W(1)),(INTYP,W(41)),(STEP1,W(44))             BACKUP15

      REAL*8 KR                                                         BACKUP16

      EQUIVALENCE (KR,R(4)),(DKRDT,DRDT(4))                             BACKUP17

      REAL*8 INTYP                                                      BACKUP18

      LOGICAL PCNTRL                                                    BACKUP19

      DATA NBAK,NGRAZ/' BACK UP0', ' GRAZE 1'/                          BACKUP20

      data avl/0.0/

C                                                                       BACKUP21

      NTRY=NBAK                                                         BACKUP22

      GO TO 100                                                         BACKUP23

C                                                                       BACKUP24

      ENTRY GRAZE(Z,RSIGN,TC)                                           BACKUP25

      NTRY=NGRAZ                                                        BACKUP26

C                                                                       BACKUP27

C     DEFINE BASE STEP SIZE                                             BACKUP28

100   STPB=STP                                                          BACKUP29

C                                                                       BACKUP30

C     DEFINE BACKUP LOCATION TOLERANCES BASED ON INTEGRATION MODE       BACKUP31

C                                                                       BACKUP32

      IF(MODE.LT.3) THEN                                                BACKUP33

         TOL=STEP1                                                      BACKUP34

      ELSEIF(MODE.EQ.3) THEN                                            BACKUP35

         TOL=DABS(E1MAX*STPB)                                           BACKUP36

      ELSE                                                              BACKUP37

         TOL=E1MAX                                                      BACKUP38

      ENDIF                                                             BACKUP39

C     USE THIS OR MINIMUM STEP SIZE WHICHEVER IS LESS                   BACKUP40

      TOL=DMIN1(TOL,FACT*E2MIN)                                         BACKUP41

C                                                                       BACKUP42

      PCNTRL=W(73).NE.0.0                                               BACKUP43

      THERE=.TRUE.                                                      BACKUP44

C      STEP TO ESTIMATED CROSSING AT TC.                                BACKUP45

C********* DIAGNOSTIC PRINTOUT                                          BACKUP46

      IF(PCNTRL) CALL PRINTR (NTRY,avl)                                 BACKUP47

      STEP=TC-T                                                         BACKUP48

      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)                     BACKUP49

      MODE=1                                                            BACKUP50

      RSTART=1                                                          BACKUP51

      TOLD=T                                                            BACKUP52

      CALL RMOVE(DROLD,DRDT,3)                                          BACKUP53

      ZDOLD=DOT(Z)                                                      BACKUP54

      CALL RKAM                                                         BACKUP55

      RSTART=1.0                                                        BACKUP56

      ZDOT=DOT(Z)                                                       BACKUP57

      IF(NTRY.EQ.NGRAZ .OR. RSIGN*ZDOT.LE.0.0) GO TO 12                 BACKUP58

C                                                                       BACKUP59

C********* FIND NEAREST INTERSECTION OF RAY WITH THE HEIGHT HS          BACKUP60

      DO 10 I=1,10                                                      BACKUP61

      STEP=-Z(1)/ZDOT                                                   BACKUP62

      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)                     BACKUP63

      IF (DABS(Z(1)).LT.PRNZTL .AND. DABS(STEP).LT.TOL) GO TO 60        BACKUP64

C********* DIAGNOSTIC PRINTOUT                                          BACKUP65

      IF(PCNTRL) CALL PRINTR(' BACK UP1',avl)                           BACKUP66

      MODE=1                                                            BACKUP67

      RSTART=1.                                                         BACKUP68

      TOLD=T                                                            BACKUP69

      CALL RMOVE(DROLD,DRDT,3)                                          BACKUP70

      ZDOLD=ZDOT                                                        BACKUP71

      CALL RKAM                                                         BACKUP72

      ZDOT=DOT(Z)                                                       BACKUP73

   10 RSTART=1.                                                         BACKUP74

C                                                                       BACKUP75

C********* FIND NEAREST CLOSEST APPROACH OF RAY TO THE HEIGHT HS        BACKUP76

12    THERE=.FALSE.                                                     BACKUP77

      DO 20 I=1,10                                                      BACKUP78

C     DO 'LOCAL' PARABOLIC FIT                                          BACKUP79

      CALL FIT3(Z,ZOLD,ZDOLD)                                           BACKUP80

      STEP=-ZDOT/D2Z                                                    BACKUP81

      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)                     BACKUP82

      IF (DABS(ZDOT).LE.PRNDZTL .AND. DABS(STEP).LT.TOL) GO TO 60       BACKUP83

C********* DIAGNOSTIC PRINTOUT                                          BACKUP84

      IF(PCNTRL) CALL PRINTR (' GRAZE 1 ',avl)                          BACKUP85

      MODE=1                                                            BACKUP86

      RSTART=1.                                                         BACKUP87

      TOLD=T                                                            BACKUP88

      CALL RMOVE(DROLD,DRDT,3)                                          BACKUP89

      ZOLD=Z(1)                                                         BACKUP90

      ZDOLD=ZDOT                                                        BACKUP91

      CALL RKAM                                                         BACKUP92

      RSTART=1.                                                         BACKUP93

      IF (D2Z*Z(1).LT.0.) GO TO 30                                      BACKUP94

      IF(KPARLEL(Z).EQ.0.0) GO TO 60                                    BACKUP95

   20 CONTINUE                                                          BACKUP96

      WRITE(3,350)                                                      BACKUP97

350   FORMAT(' ****** COULDN''T FIND CLOSEST APPROACH IN 10 STEPS')     BACKUP98

      GO TO 60                                                          BACKUP99

C                                                                       BACKU100

   30 CONTINUE                                                          BACKU101

C********* DIAGNOSTIC PRINTOUT                                          BACKU102

      IF(PCNTRL) CALL PRINTR (' BACK UP2',avl)                          BACKU103

      MODE=1                                                            BACKU104

C********* ESTIMATE DISTANCE TO NEAREST INTERSECTION GOING THE RIGHT    BACKU105

C********* DIRECTION OF RAY WITH HEIGHT HS                              BACKU106

      CALL FIT3(Z,ZOLD,ZDOLD)                                           BACKU107

      STEP=(-ZDOT+RSIGN*RAD1)/D2Z                                       BACKU108

      RSTART=1.                                                         BACKU109

      CALL RKAM                                                         BACKU110

      RSTART=1.                                                         BACKU111

C********* FIND NEAREST INTERSECTION OF RAY WITH HEIGHT HS              BACKU112

      DO 40 I=1,10                                                      BACKU113

      ZDOT=DOT(Z)                                                       BACKU114

      STEP=-Z(1)/ZDOT                                                   BACKU115

      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)                     BACKU116

      IF (DABS(Z(1)).LT.PRNZTL .AND. DABS(STEP).LT.TOL) GO TO 50        BACKU117

C********* DIAGNOSTIC PRINTOUT                                          BACKU118

      IF(PCNTRL) CALL PRINTR (' BACK UP3',avl)                          BACKU119

      MODE=1                                                            BACKU120

      RSTART=1.                                                         BACKU121

      CALL RKAM                                                         BACKU122

   40 RSTART=1.                                                         BACKU123

   50 THERE=.TRUE.                                                      BACKU124

C********* RESET STANDARD MODE AND INTEGRATION TYPE                     BACKU125

   60 MODE=INTYP                                                        BACKU126

      STEP=STPB                                                         BACKU127

      RETURN                                                            BACKU128

         END                                                            BACKU129



