      PROGRAM NITIAL                                                    NITI001

C          SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE     NITI002

      DIMENSION MFLD(2)                                                 NITI003

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN                   NITI004

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             NITI005

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR,     NITI006

     1             LPOLAR,SGN                                           NITI007

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

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  NITI009

      COMMON /YY/ MODY,Y(16) /ZZ/ MODZ,Z(4)                             NITI010

      COMMON R(20),T,STP,DRDT(20) /WW/ ID(10),W0,W(400)                 NITI011

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

     1 (TLON,W(5)),(F,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),       NITI013

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

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

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

     5,(INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),(STEP1,W(44)),       NITI017

     6 (STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),(SKIP,W(71)),        NITI018

     7 (RAYSET,W(72)),(PLT,W(81)),(PERT,W(150))                         NITI019

      LOGICAL SPACE,NEWWR,NEWWP,PENET                                   NITI020

      REAL N2,N2I,LOGTEN,K,MAXSTP,INTYP,MAXERR,MU                       NITI021

      COMPLEX PNP,POLAR,LPOLAR                                          NITI022

      NDATE=IDATE(0)                                                    NITI025

      SECOND=KLOCK(0)*.001                                              NITI026

      KOLL=4H  NO                                                       NITI027

      IF (COLL.NE.0.) KOLL=4HWITH                                       NITI028

C********* CONSTANTS                                                    NITI029

      PI=3.1415926536                                                   NITI030

      PIT2=2.*PI                                                        NITI031

      PID2=PI/2.                                                        NITI032

      DEGS=180./PI                                                      NITI033

      RAD=PI/180.                                                       NITI034

      C=2.997925E5                                                      NITI035

      K=2.81785E-15*C**2/PI                                             NITI036

      LOGTEN=ALOG(10.)                                                  NITI037

C********* INITIALIZE SOME VARIABLES IN THE W ARRAY                     NITI038

      DO 5 NW=1,400                                                     NITI039

    5 W(NW)=0.                                                          NITI040

      PLON=0.                                                           NITI041

      PLAT=PID2                                                         NITI042

      EARTHR=6370.                                                      NITI043

      INTYP=3.                                                          NITI044

      MAXERR=1.E-4                                                      NITI045

      ERATIO=50.                                                        NITI046

      STEP1=1.                                                          NITI047

      STPMAX=100.                                                       NITI048

      STPMIN=1.E-8                                                      NITI049

      FACTR=0.5                                                         NITI050

      MAXSTP=1000.                                                      NITI051

      HOP=1.                                                            NITI052

C********* READ W ARRAY AND PRINT NON-ZERO VALUES                       NITI053

   10 CALL READ W                                                       NITI054

      F=BETA=AZ1=0.                                                     NITI055

      IF (SKIP.EQ.0.) SKIP=MAXSTP                                       NITI056

  12  RAY=SIGN(1.,RAY)                                                  NITI057

      NTYP=2.+FIELD*RAY                                                 NITI058

      GO TO (13,14,15), NTYP                                            NITI059

   13 MFLD(1)=8HEXTRAORD                                                NITI060

      MFLD(2)=5HINARY                                                   NITI061

      GO TO 16                                                          NITI062

   14 MFLD(1)=8HNO FIELD                                                NITI063

      MFLD(2)=1H                                                        NITI064

      GO TO 16                                                          NITI065

   15 MFLD(1)=8HORDINARY                                                NITI066

      MFLD(2)=8H                                                        NITI067

   16 MODSAV=MODX(2)                                                    NITI068

      IF (PERT.EQ.0.) MODX(2)=6H                                        NITI069

      IF (RAYSET.NE.0.) PUNCH 2000, ID,MODX(1),(W(NW),NW=101,107),      NITI070

     1 MODX(2),(W(NW),NW=151,157),MODY,(W(NW),NW=201,207),              NITI071

     2 MODZ,(W(NW),NW=251,257)                                          NITI072

 2000 FORMAT (10A8,4(/A8,2X,7E10.3))                                    NITI073

      PRINT 1000, ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL              NITI074

 1000 FORMAT (1H1,10A8,25X,A8/4(1X,A8),24X,3A8,1X,A8,A5,1X,A4,          NITI075

     1 11H COLLISIONS/)                                                 NITI076

      PRINT 1050                                                        NITI077

 1050 FORMAT (85H INITIAL VALUES FOR THE W ARRAY -- ALL ANGLES IN RADIANNITI078

     1S, ONLY NONZERO VALUES PRINTED/)                                  NITI079

      DO 17 NW=1,400                                                    NITI080

      IF (W(NW).NE.0.) PRINT 1700, NW,W(NW)                             NITI081

 1700 FORMAT (I4,E19.11)                                                NITI082

   17 CONTINUE                                                          NITI083

C********* LET SUBROUTINES PRINTR AND RAYPLT KNOW THERE IS A NEW W ARRAYNITI084

      NEWWP=.TRUE.                                                      NITI085

      NEWWR=.TRUE.                                                      NITI086

C********* INITIALIZE PARAMETERS FOR INTEGRATION SUBROUTINE RKAM        NITI087

      N=6                                                               NITI088

      DO 20 NR=7,20                                                     NITI089

      IF (W(50+NR).NE.0.) N=N+1                                         NITI090

   20 CONTINUE                                                          NITI091

      MODE=INTYP                                                        NITI092

      STEP=STEP1                                                        NITI093

      E1MAX=MAXERR                                                      NITI094

      E1MIN=MAXERR/ERATIO                                               NITI095

      E2MAX=STPMAX                                                      NITI096

      E2MIN=STPMIN                                                      NITI097

      FACT=FACTR                                                        NITI098

C********* DETERMINE TRANSMITTER LOCATION IN COMPUTATIONAL COORDINATE   NITI099

C********* SYSTEM (GEOMAGNETIC COORDINATES IF DIPOLE FIELD IS USED)     NITI100

      R0=EARTHR+XMTRH                                                   NITI101

      SP=SIN (PLAT)                                                     NITI102

      CP=SIN (PID2-PLAT)                                                NITI103

      SDPH=SIN (TLON-PLON)                                              NITI104

      CDPH=SIN (PID2-(TLON-PLON))                                       NITI105

      SL=SIN (TLAT)                                                     NITI106

      CL=SIN (PID2-TLAT)                                                NITI107

      ALPHA=ATAN2(-SDPH*CP,-CDPH*CP*SL+SP*CL)                           NITI108

      TH0=ACOS (CDPH*CP*CL+SP*SL)                                       NITI109

      PH0=ATAN2(SDPH*CL,CDPH*SP*CL-CP*SL)                               NITI110

C********* LOOP ON FREQUENCY, AZIMUTH ANGLE, AND ELEVATION ANGLE        NITI111

      NFREQ=1                                                           NITI112

      IF (FSTEP.NE.0.) NFREQ=(FEND-FBEG)/FSTEP+1.5                      NITI113

      NAZ=1                                                             NITI114

      IF (AZSTEP.NE.0.) NAZ=(AZEND-AZBEG)/AZSTEP+1.5                    NITI115

      NBETA=1                                                           NITI116

      IF (ELSTEP.NE.0.) NBETA=(ELEND-ELBEG)/ELSTEP+1.5                  NITI117

      DO 50 NF=1,NFREQ                                                  NITI118

      F=FBEG+(NF-1)*FSTEP                                               NITI119

      DO 45 J=1,NAZ                                                     NITI120

      AZ1=AZBEG+(J-1)*AZSTEP                                            NITI121

      AZA=AZ1*DEGS                                                      NITI122

      GAMMA=PI-AZ1+ALPHA                                                NITI123

      SGAMMA=SIN (GAMMA)                                                NITI124

      CGAMMA=SIN (PID2-GAMMA)                                           NITI125

      DO 40 I=1,NBETA                                                   NITI126

      BETA=ELBEG+(I-1)*ELSTEP                                           NITI127

      EL=BETA*DEGS                                                      NITI128

      CBETA=SIN (PID2-BETA)                                             NITI129

      R(1)=R0                                                           NITI130

      R(2)=TH0                                                          NITI131

      R(3)=PH0                                                          NITI132

      R(4)=SIN (BETA)                                                   NITI133

      R(5)=CBETA*CGAMMA                                                 NITI134

      R(6)=CBETA*SGAMMA                                                 NITI135

      T=0.                                                              NITI136

      RSTART=1.                                                         NITI137

C     SGN=1.  (NEED FOR RAY TRACING IN COMPLEX SPACE.)                  NITI138

C********* ALLOW IONOSPHERIC MODEL SUBROUTINES TO READ AND PRINT DATA   NITI139

      CALL RINDEX                                                       NITI140

      IF (I.NE.1.AND.NPAGE.LT.3.AND.LINES.LE.17) GO TO 25               NITI141

      NPAGE=LINES=0                                                     NITI142

      PRINT 1000, ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL              NITI143

      PRINT 2400, F,AZA                                                 NITI144

 2400 FORMAT (18X,11HFREQUENCY =,F12.6,37H MHZ, AZIMUTH ANGLE OF TRANSMINITI145

     1SSION =,F12.6,4H DEG)                                             NITI146

   25 NPAGE=NPAGE+1                                                     NITI147

      PRINT 2500, EL                                                    NITI148

 2500 FORMAT (/31X,33HELEVATION ANGLE OF TRANSMISSION =,F12.6,4H DEG/)  NITI149

      IF (N2.GT.0.) GO TO 30                                            NITI150

      CALL ELECTX                                                       NITI151

      FN=SIGN (SQRT (ABS (X))*F,X)                                      NITI152

      PRINT 2900, FN                                                    NITI153

 2900 FORMAT (58H0TRANSMITTER IN EVANESCENT REGION, TRANSMISSION IMPOSSINITI154

     1BLE/20H0PLASMA FREQUENCY = ,E17.10)                               NITI155

      GO TO 44                                                          NITI156

   30 MU=SQRT (N2/(R(4)**2+R(5)**2+R(6)**2))                            NITI157

      DO 34 NN=4,6                                                      NITI158

   34 R(NN)=R(NN)*MU                                                    NITI159

      DO 35 NN=7,N                                                      NITI160

   35 R(NN)=0.                                                          NITI161

      CALL TRACE                                                        NITI162

      OSEC=SECOND                                                       NITI163

      SECOND=KLOCK(0)*.001                                              NITI164

      DIFF=SECOND-OSEC                                                  NITI165

      PRINT 3500, DIFF                                                  NITI166

 3500 FORMAT (36X,26HTHIS RAY CALCULATION TOOK ,F8.3,4H SEC)            NITI167

      IF (PENET.AND.ONLY.NE.0..AND.IHOP.EQ.1) GO TO 44                  NITI168

   40 CONTINUE                                                          NITI169

   44 IF (PLT.NE.0.) CALL ENDPLT                                        NITI170

   45 CONTINUE                                                          NITI171

      IF(PENET.AND.ONLY.NE.0..AND.IHOP.EQ.1.AND.NAZ.EQ.1.AND.NBETA.EQ.1)NITI172

     1 GO TO 55                                                         NITI173

   50 CONTINUE                                                          NITI174

  55  IF (RAYSET.NE.0.) PUNCH 5000                                      NITI175

 5000 FORMAT (78X,1H-)                                                  NITI176

      MODX(2)=MODSAV                                                    NITI177

      GO TO 10                                                          NITI178

         END                                                            NITI179-

      SUBROUTINE READ W                                                 READ001

C                        READS W ARRAY                                  READ002

C       A 1 IN THE FOLLOWING COLUMNS WILL MAKE THE DESCRIBED CONVERSIONSREAD003

COL 18   DEGREES TO RADIANS                                             READ004

COL 19   GREAT CIRCLE DISTANCE IN KM TO RADIANS                         READ005

COL 20   NAUTICAL MILES TO KM                                           READ006

COL 21   FEET TO KM                                                     READ007

C                                                                       READ008

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,DUM(3)                       READ009

      COMMON /WW/ ID(10),W0,W(400)                                      READ010

      EQUIVALENCE (EARTHR,W(2))                                         READ011

      INTEGER DEG,FEET                                                  READ012

      READ 1000, ID                                                     READ013

 1000 FORMAT (10A8)                                                     READ014

      IF (EOF,60) 3,4                                                   READ015

    3 CALL EXIT                                                         READ016

    4 READ 1100, NW,W(NW),DEG,KM,NM,FEET                                READ017

 1100 FORMAT (I3,E14.7,5I1)                                             READ018

      IF (NW.EQ.0) GO TO 10                                             READ019

      IF (NW.GT.0.AND.NW.LE.400) GO TO 5                                READ020

      PRINT 4000, NW                                                    READ021

 4000 FORMAT(15H1THE SUBSCRIPT ,I3,77H ON THE W-ARRAY INPUT IS OUT OF BOREAD022

     1UNDS.  ALLOWABLE VALUES ARE 1 THROUGH 400. )                      READ023

      CALL EXIT                                                         READ024

    5 IF (DEG.NE.0.) W(NW)=W(NW)*RAD                                    READ025

      IF (KM.NE.0) W(NW)=W(NW)/EARTHR                                   READ026

      IF (NM.NE.0) W(NW)=W(NW)*1.852                                    READ027

      IF (FEET.NE.0) W(NW)=W(NW)*3.048006096E-4                         READ028

      GO TO 4                                                           READ029

   10 RETURN                                                            READ030

         END                                                            READ031-

      SUBROUTINE TRACE                                                  TRAC001

C          CALCULATES THE RAY PATH                                      TRAC002

      DIMENSION ROLD(20),DROLD(20)                                      TRAC003

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

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             TRAC005

      COMMON /TRAC/ GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT               TRAC006

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,N2,PNP(10),POLAR,LPOLAR   TRAC007

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  TRAC008

      COMMON R(20),T,STP,DRDT(20) /WW/ ID(10),W0,W(400)                 TRAC009

      LOGICAL SPACE,HOME,WASNT,UNDRGD,GROUND,PERIGE,THERE,MINDIS,NEWWR, TRAC010

     1        NEWWP,PENET,NEWRAY,WAS                                    TRAC011

      REAL MAXSTP                                                       TRAC012

      COMPLEX N2,PNP,POLAR,LPOLAR                                       TRAC013

      EQUIVALENCE (EARTHR,W(2)),(RCVRH,W(20)),(HOP,W(22)),(MAXSTP,W(23))TRAC014

     1,(SKIP,W(71)),(RAYSET,W(72)),(PLT,W(81))                          TRAC015

      NHOP=HOP                                                          TRAC016

      MAX=MAXSTP                                                        TRAC017

      NSKIP=SKIP                                                        TRAC018

      RSTART=1.                                                         TRAC019

      CALL HAMLTN                                                       TRAC020

      HOME=DRDT(1)*(R(1)-EARTHR-RCVRH).GE.0.                            TRAC021

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

C********* RAYSET AND TELLS RAYPLT TO START A NEW RAY                   TRAC023

      IHOP=0                                                            TRAC024

      CALL PRINTR (8HXMTR    ,0.)                                       TRAC025

      IF (PLT.NE.0.) CALL RAYPLT                                        TRAC026

      HTMAX=0.                                                          TRAC027

      NEWRAY=.TRUE.                                                     TRAC028

      THERE=R(1)-EARTHR.EQ.RCVRH                                        TRAC029

C********* LOOP ON NUMBER OF HOPS                                       TRAC030

   10 IHOP=IHOP+1                                                       TRAC031

      IF (IHOP.GT.NHOP) RETURN                                          TRAC032

      PENET=.FALSE.                                                     TRAC033

      APHT=RCVRH                                                        TRAC034

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

      DO 79 J=1,MAX                                                     TRAC036

      H=R(1)-EARTHR                                                     TRAC037

      IF (ABS(H-RCVRH).GT.ABS(APHT-RCVRH)) APHT=H                       TRAC038

      HTMAX=AMAX1(H,HTMAX)                                              TRAC039

      IF (.NOT.SPACE) GO TO 12                                          TRAC040

      CALL REACH                                                        TRAC041

      RSTART=1.                                                         TRAC042

      H=R(1)-EARTHR                                                     TRAC043

      IF (ABS(H-RCVRH).GT.ABS(APHT-RCVRH)) APHT=H                       TRAC044

      HTMAX=AMAX1(H,HTMAX)                                              TRAC045

      IF (.NOT.SPACE) GO TO 12                                          TRAC046

      IF (PERIGE) CALL PRINTR (8HPERIGEE ,0.)                           TRAC047

      IF (THERE) GO TO 51                                               TRAC048

      IF (MINDIS) GO TO 40                                              TRAC049

      IF (GROUND) GO TO 60                                              TRAC050

      IF (PLT.NE.0.) CALL RAYPLT                                        TRAC051

      IF (PERIGE) GO TO 79                                              TRAC052

   12 DO 13 L=1,N                                                       TRAC053

      ROLD(L)=R(L)                                                      TRAC054

   13 DROLD(L)=DRDT(L)                                                  TRAC055

      TOLD=T                                                            TRAC056

      WAS=THERE                                                         TRAC057

      CALL RKAM                                                         TRAC058

      H=R(1)-EARTHR                                                     TRAC059

      THERE=.FALSE.                                                     TRAC060

      WASNT=.NOT.HOME                                                   TRAC061

      HOME=DRDT(1)*(H-RCVRH).GE.0.                                      TRAC062

      TMP=(DRDT(1)-DROLD(1))*(T-TOLD)                                   TRAC063

      SMT=0.                                                            TRAC064

      IF (TMP.NE.0.) SMT=0.5*(R(1)-ROLD(1)+0.5*TMP)**2/ABS(TMP)         TRAC065

      IF (((H-RCVRH)*(ROLD(1)-EARTHR-RCVRH).LT.0..AND..NOT.WAS).OR.     TRAC066

     1 (WAS.AND.DRDT(1)*DROLD(1).LT.0..AND.HOME)) GO TO 50              TRAC067

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

      IF (H.LT.0..OR.DRDT(1).GT.0..AND.DROLD(1).LT.0..AND.SMT.GT.H)     TRAC069

     1 GO TO 20                                                         TRAC070

      IF (DROLD(1).LT.0..AND.DRDT(1).GT.0.) CALL PRINTR(8HPERIGEE ,0.)  TRAC071

      IF (DROLD(1).GT.0..AND.DRDT(1).LT.0.) CALL PRINTR(8HAPOGEE  ,0.)  TRAC072

      IF (DROLD(2)*DRDT(2).LT.0.) CALL PRINTR(8HMAX LAT ,0.)            TRAC073

      IF (DROLD(3)*DRDT(3).LT.0.) CALL PRINTR(8HMAX LONG,0.)            TRAC074

      DO 14 I=4,6                                                       TRAC075

      IF (ROLD(I)*R(I).LT.0.) CALL PRINTR(8HWAVE REV,0.)                TRAC076

   14 CONTINUE                                                          TRAC077

      GO TO 75                                                          TRAC078

C********* RAY WENT UNDERGROUND                                         TRAC079

   20 CALL BACK UP(0.)                                                  TRAC080

      GO TO 60                                                          TRAC081

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

   30 CALL GRAZE(RCVRH)                                                 TRAC083

      IF (THERE) GO TO 51                                               TRAC084

   40 DRDT(1)=0.                                                        TRAC085

      HPUNCH=R(1)-EARTHR                                                TRAC086

      CALL PRINTR(8HMIN DIST,RAYSET)                                    TRAC087

      IF (PLT.NE.0.) CALL RAYPLT                                        TRAC088

      IF (IHOP.GE.NHOP) RETURN                                          TRAC089

      IHOP=IHOP+1                                                       TRAC090

      CALL PRINTR (8HMIN DIST,RAYSET)                                   TRAC091

      GO TO 89                                                          TRAC092

C********* RAY CROSSED RECEIVER HEIGHT                                  TRAC093

   50 CALL BACK UP(RCVRH)                                               TRAC094

      THERE=.TRUE.                                                      TRAC095

   51 R(1)=EARTHR+RCVRH                                                 TRAC096

      HTMAX=AMAX1(RCVRH,HTMAX)                                          TRAC097

      HPUNCH=APHT                                                       TRAC098

      CALL PRINTR(8HRCVR    ,RAYSET)                                    TRAC099

      IF (PLT.NE.0.) CALL RAYPLT                                        TRAC100

      IF (RCVRH.NE.0.) GO TO 89                                         TRAC101

      IF (IHOP.GE.NHOP) RETURN                                          TRAC102

      IHOP=IHOP+1                                                       TRAC103

      APHT=RCVRH                                                        TRAC104

C********* GROUND REFLECT                                               TRAC105

   60 R(1)=EARTHR                                                       TRAC106

      IF (ABS(RCVRH).GT.ABS(APHT-RCVRH)) APHT=0.                        TRAC107

      R(4)=ABS (R(4))                                                   TRAC108

      DRDT(1)=ABS (DRDT(1))                                             TRAC109

      RSTART=1.                                                         TRAC110

      HPUNCH=HTMAX                                                      TRAC111

      CALL PRINTR(8HGRND REF,RAYSET)                                    TRAC112

      HTMAX=0.                                                          TRAC113

      IF (RCVRH.NE.0.) GO TO 65                                         TRAC114

      THERE=.TRUE.                                                      TRAC115

      HPUNCH=APHT                                                       TRAC116

      CALL PRINTR (8HRCVR    ,RAYSET)                                   TRAC117

      GO TO 89                                                          TRAC118

   65 H=0.                                                              TRAC119

      THERE=.FALSE.                                                     TRAC120

C*********                                                              TRAC121

   75 IF (PLT.NE.0.) CALL RAYPLT                                        TRAC122

      IF (H.GT.HMAX.AND.H.GT.RCVRH.AND.DRDT(1).GT.0.) GO TO 90          TRAC123

      IF (MOD(J,NSKIP).EQ.0) CALL PRINTR(8H        ,0.)                 TRAC124

   79 CONTINUE                                                          TRAC125

C********* EXCEEDED MAXIMUM NUMBER OF STEPS                             TRAC126

      HPUNCH=H                                                          TRAC127

      CALL PRINTR(8HSTEP MAX,RAYSET)                                    TRAC128

      RETURN                                                            TRAC129

C*********                                                              TRAC130

   89 HOME=.TRUE.                                                       TRAC131

      GO TO 10                                                          TRAC132

C********* RAY PENETRATED                                               TRAC133

   90 PENET=.TRUE.                                                      TRAC134

      HPUNCH=H                                                          TRAC135

      CALL PRINTR(8HPENETRAT,RAYSET)                                    TRAC136

      RETURN                                                            TRAC137

         END                                                            TRAC138

      SUBROUTINE REACH                                                  REAC001

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)                                                  REAC 16

      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-

      SUBROUTINE BACK UP(HS)                                            BACK001

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

      COMMON /TRAC/ GROUND,PERIGE,THERE,MINDIS,NEWRAY,SMT               BACK003

      COMMON R(20),T,STP,DRDT(20) /WW/ ID(10),W0,W(400)                 BACK004

      EQUIVALENCE (EARTHR,W(2)),(INTYP,W(41)),(STEP1,W(44))             BACK005

      REAL INTYP                                                        BACK006

      LOGICAL GROUND,PERIGE,THERE,MINDIS,NEWRAY,HOME                    BACK007

C********* DIAGNOSTIC PRINTOUT                                          BACK008

C        CALL PRINTR (8HBACK UP0,0.)                                    BACK009

C********* GOING AWAY FROM THE HEIGHT HS                                BACK010

      HOME=DRDT(1)*(R(1)-EARTHR-HS).GE.0.                               BACK011

      IF (HS.GT.0..AND..NOT.HOME.OR.HS.EQ.0..AND.DRDT(1).GT.0.) GO TO 30BACK012

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

      DO 10 I=1,10                                                      BACK014

      STEP=-(R(1)-EARTHR-HS)/DRDT(1)                                    BACK015

      STEP=SIGN(AMIN1(ABS(STP),ABS(STEP)),STEP)                         BACK016

      IF (ABS(R(1)-EARTHR-HS).LT..5E-4.AND.ABS(STEP).LT.1.) GO TO 60    BACK017

C********* DIAGNOSTIC PRINTOUT                                          BACK018

C        CALL PRINTR(8HBACK UP1,0.)                                     BACK019

      MODE=1                                                            BACK020

      RSTART=1.                                                         BACK021

      CALL RKAM                                                         BACK022

   10 RSTART=1.                                                         BACK023

C                                                                       BACK024

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

      ENTRY GRAZE                                                       BACK026

      THERE=.FALSE.                                                     BACK027

C********* DIAGNOSTIC PRINTOUT                                          BACK028

C        CALL PRINTR (8HGRAZE 0 ,0.)                                    BACK029

      IF (SMT.GT.ABS(R(1)-EARTHR-HS)) GO TO 30                          BACK030

      DO 20 I=1,10                                                      BACK031

      STEP=-R(4)/DRDT(4)                                                BACK032

      STEP=SIGN(AMIN1(ABS(STP),ABS(STEP)),STEP)                         BACK033

      IF (ABS(R(4)).LE.1.E-6.AND.ABS(STEP).LT.1.) GO TO 60              BACK034

C********* DIAGNOSTIC PRINTOUT                                          BACK035

C        CALL PRINTR (8HGRAZE 1 ,0.)                                    BACK036

      MODE=1                                                            BACK037

      RSTART=1.                                                         BACK038

      CALL RKAM                                                         BACK039

      RSTART=1.                                                         BACK040

      IF (DRDT(4)*(R(1)-EARTHR-HS).LT.0.) GO TO 30                      BACK041

      IF(R(5).EQ.0..AND.R(6).EQ.0.) GO TO 60                            BACK042

   20 CONTINUE                                                          BACK043

C********* IF A CLOSEST APPROACH COULD NOT BE FOUND IN 10 STEPS, IT     BACK044

C********* PROBABLY MEANS THAT THE RÁÙ INTERSECTS THE HEIGHT HS         BACK045

   30 CONTINUE                                                          BACK046

C********* DIAGNOSTIC PRINTOUT                                          BACK047

C        CALL PRINTR (8HBACK UP2,0.)                                    BACK048

      MODE=1                                                            BACK049

C********* ESTIMATE DISTANCE TO NEAREST INTERSECTION OF RAY WITH HEIGHT BACK050

C********* HS BEHIND THE PRESENT RAY POINT                              BACK051

      H=R(1)-EARTHR                                                     BACK0519

      STEP=(-R(4)-SQRT(R(4)**2-2.*(H-HS)*R(4)/DRDT(1)*DRDT(4)))/DRDT(4) BACK052

C     STEP=(-R(4)-SQRT(R(4)**2-2.*(R(1)-EARTHR-HS)*DRDT(4)))/DRDT(4)    BACK052

      RSTART=1.                                                         BACK053

      CALL RKAM                                                         BACK054

      RSTART=1.                                                         BACK055

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

      DO 40 I=1,10                                                      BACK057

      STEP=-(R(1)-EARTHR-HS)/DRDT(1)                                    BACK058

      STEP=SIGN(AMIN1(ABS(STP),ABS(STEP)),STEP)                         BACK059

      IF(ABS(R(1)-EARTHR-HS).LT..5E-4.AND.ABS(STEP).LT.1.) GO TO 50     BACK060

C     IF (ABS(R(1)-EARTHR-HS).LT..5E-4.AND.ABS(STEP).LT.1.) GO TO 60    BACK060

C********* DIAGNOSTIC PRINTOUT                                          BACK061

C        CALL PRINTR (8HBACK UP3,0.)                                    BACK062

      MODE=1                                                            BACK063

      RSTART=1.                                                         BACK064

      CALL RKAM                                                         BACK065

   40 RSTART=1.                                                         BACK066

   50 THERE=.TRUE.                                                      BACK067

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

   60 MODE=INTYP                                                        BACK069

      STEP=STEP1                                                        BACK070

      RETURN                                                            BACK071

         END                                                            BACK072-

      SUBROUTINE POL CAR                                                POLC001

      DIMENSION X0(6),X(6),R0(4)                                        POLC002

      COMMON R(6) /COORD/ S                                             POLC003

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                POLC004

C                                                                       POLC005

C          CONVERTS SPHERICAL COORDINATES TO CARTESIAN                  POLC006

      IF (R(5).EQ.0..AND.R(6).EQ.0.) GO TO 1                            POLC007

      VERT=0.                                                           POLC008

      SINA=SIN(R(2))                                                    POLC009

      COSA=SIN(PID2-R(2))                                               POLC010

      SINP=SIN(R(3))                                                    POLC011

      COSP=SIN(PID2-R(3))                                               POLC012

      X0(1)=R(1)*SINA*COSP                                              POLC013

      X0(2)=R(1)*SINA*SINP                                              POLC014

      X0(3)=R(1)*COSA                                                   POLC015

      X(4)=R(4)*SINA*COSP+R(5)*COSA*COSP-R(6)*SINP                      POLC016

      X(5)=R(4)*SINA*SINP+R(5)*COSA*SINP+R(6)*COSP                      POLC017

      X(6)=R(4)*COSA-R(5)*SINA                                          POLC018

      RETURN                                                            POLC019

C          VERTICAL INCIDENCE                                           POLC020

    1 VERT=1.                                                           POLC021

      R0(1)=R(1)                                                        POLC022

      R0(2)=R(2)                                                        POLC023

      R0(3)=R(3)                                                        POLC024

      R0(4)=SIGN (1.,R(4))                                              POLC025

      RETURN                                                            POLC026

C                                                                       POLC027

C          STEPS THE RAY A DISTANCE S, AND THEN                         POLC028

C          CONVERTS CARTESIAN COORDINATES TO SPHERICAL COORDINATES      POLC029

      ENTRY CAR POL                                                     POLC030

      IF (VERT.NE.0.) GO TO 2                                           POLC031

      X(1)=X0(1)+S*X(4)                                                 POLC032

      X(2)=X0(2)+S*X(5)                                                 POLC033

      X(3)=X0(3)+S*X(6)                                                 POLC034

      TEMP=SQRT(X(1)**2+X(2)**2)                                        POLC035

      R(1)=SQRT(X(1)**2+X(2)**2+X(3)**2)                                POLC036

      R(2)=ATAN2(TEMP,X(3))                                             POLC037

      R(3)=ATAN2(X(2),X(1))                                             POLC038

      R(4)=(X(1)*X(4)+X(2)*X(5)+X(3)*X(6))/R(1)                         POLC039

      R(5)=(X(3)*(X(1)*X(4)+X(2)*X(5))-(X(1)**2+X(2)**2)*X(6))/         POLC040

     1 (R(1)*TEMP)                                                      POLC041

      R(6)=(X(1)*X(5)-X(2)*X(4))/TEMP                                   POLC042

      RETURN                                                            POLC043

C          VERTICAL INCIDENCE                                           POLC044

    2 R(1)=R0(1)+R0(4)*S                                                POLC045

      R(2)=R0(2)                                                        POLC046

      R(3)=R0(3)                                                        POLC047

      R(4)=R0(4)                                                        POLC048

      R(5)=0.                                                           POLC049

      R(6)=0.                                                           POLC050

      RETURN                                                            POLC051

         END                                                            POLC 52-

      SUBROUTINE PRINTR(NWHY,CARD)                                      PRIN001

C          PRINTS OUTPUT AND PUNCHES RAYSETS WHEN REQUESTED             PRIN002

      DIMENSION G(3,3),G1(3,3),TYPE(3),HEADR1(20),HEADR2(20),UNITS(20), PRIN003

     1          HEAD1(20),HEAD2(20),UNIT(20),RPRINT(20),NPR(20)         PRIN004

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,DUM(3)                       PRIN005

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             PRIN006

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR(2),  PRIN007

     1      LPOLAR(2),SGN                                               PRIN008

C    1             LPOLAR(2)                                            PRIN008

      COMMON R(20),T /WW/ ID(10),W0,W(400)                              PRIN009

      EQUIVALENCE (THETA,R(2)),(PHI,R(3))                               PRIN010

      EQUIVALENCE (EARTHR,W(2)),(XMTRH,W(3)),(TLAT,W(4)),(TLON,W(5)),   PRIN011

     1 (F,W(6)),(AZ1,W(10)),(BETA,W(14)),(RCVRH,W(20)),(HOP,W(22)),     PRIN012

     2 (PLAT,W(24)),(PLON,W(25)),(RAYSET,W(72))                         PRIN013

      LOGICAL SPACE,NEWWR,NEWWP,PENET                                   PRIN014

      REAL N2,N2I,LPOLAR                                                PRIN015

      COMPLEX PNP                                                       PRIN016

      DATA (TYPE=1HX,1HN,1HO)                                           PRIN017

     2,(HEADR1(7)=6H  PHAS),(HEADR2(7)=6HE PATH),(UNITS(7)=6H   KM ),   PRIN018

     3 (HEADR1(8)=6H  ABSO),(HEADR2(8)=6HRPTION),(UNITS(8)=6H   DB ),   PRIN019

     4 (HEADR1(9)=6H   DOP),(HEADR2(9)=6HPLER  ),(UNITS(9)=6H  C/S ),   PRIN020

     5 (HEADR1(10)=6H PATH ),(HEADR2(10)=6HLENGTH),(UNITS(10)=6H   KM ) PRIN023

      CALL RINDEX                                                       PRIN024

      IF (.NOT.NEWWP) GO TO 10                                          PRIN025

C********* NEW W ARRAY -- REINITIALIZE                                  PRIN026

      NEWWP=.FALSE.                                                     PRIN027

      SPL=SIN (PLON-TLON)                                               PRIN028

      CPL=SIN (PID2-(PLON-TLON))                                        PRIN029

      SP=SIN (PLAT)                                                     PRIN030

      CP=SIN (PID2-PLAT)                                                PRIN031

      SL=SIN (TLAT)                                                     PRIN032

      CL=SIN (PID2-TLAT)                                                PRIN033

C********* MATRIX TO ROTATE COORDINATES                                 PRIN034

      G(1,1)=CPL*SP*CL-CP*SL                                            PRIN035

      G(1,2)=SPL*SP                                                     PRIN036

      G(1,3)=-SL*SP*CPL-CL*CP                                           PRIN037

      G(2,1)=-SPL*CL                                                    PRIN038

      G(2,2)=CPL                                                        PRIN039

      G(2,3)=SL*SPL                                                     PRIN040

      G(3,1)=CL*CP*CPL+SP*SL                                            PRIN041

      G(3,2)=CP*SPL                                                     PRIN042

      G(3,3)=-SL*CP*CPL+SP*CL                                           PRIN043

      DENM=G(1,1)*G(2,2)*G(3,3)+G(1,2)*G(3,1)*G(2,3)+G(2,1)*G(3,2)*G(1,3PRIN044

     1)-G(2,2)*G(3,1)*G(1,3)-G(1,2)*G(2,1)*G(3,3)-G(1,1)*G(3,2)*G(2,3)  PRIN045

C********* THE MATRIX G1 IS THE INVERSE OF THE MATRIX G                 PRIN046

      G1(1,1)=(G(2,2)*G(3,3)-G(3,2)*G(2,3))/DENM                        PRIN047

      G1(1,2)=(G(3,2)*G(1,3)-G(1,2)*G(3,3))/DENM                        PRIN048

      G1(1,3)=(G(1,2)*G(2,3)-G(2,2)*G(1,3))/DENM                        PRIN049

      G1(2,1)=(G(3,1)*G(2,3)-G(2,1)*G(3,3))/DENM                        PRIN050

      G1(2,2)=(G(1,1)*G(3,3)-G(3,1)*G(1,3))/DENM                        PRIN051

      G1(2,3)=(G(2,1)*G(1,3)-G(1,1)*G(2,3))/DENM                        PRIN052

      G1(3,1)=(G(2,1)*G(3,2)-G(3,1)*G(2,2))/DENM                        PRIN053

      G1(3,2)=(G(3,1)*G(1,2)-G(1,1)*G(3,2))/DENM                        PRIN054

      G1(3,3)=(G(1,1)*G(2,2)-G(2,1)*G(1,2))/DENM                        PRIN055

      R0=EARTHR+XMTRH                                                   PRIN056

C********* CARTESIAN COORDINATES OF TRANSMITTER                         PRIN057

      XR=R0*G(1,1)                                                      PRIN058

      YR=R0*G(2,1)                                                      PRIN059

      ZR=R0*G(3,1)                                                      PRIN060

      CTHR=G(3,1)                                                       PRIN061

      STHR=SIN (ACOS (CTHR))                                            PRIN062

      PHIR=ATAN2(YR,XR)                                                 PRIN063

      ALPH=ATAN2(G(3,2),G(3,3))                                         PRIN064

C*********                                                              PRIN065

      NR=6                                                              PRIN066

      NP=0                                                              PRIN067

      DO 7 NN=7,20                                                      PRIN068

      IF (W(NN+50).EQ.0.) GO TO 7                                       PRIN069

C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED             PRIN070

C********* NR IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED     PRIN071

      NR=NR+1                                                           PRIN072

      IF (W(NN+50).NE.2.) GO TO 7                                       PRIN073

C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED AND PRINTED.PRIN074

C********* NP IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED AND PRIN075

C********* PRINTED                                                      PRIN076

      NP=NP+1                                                           PRIN077

C********* SAVE THE INDEX OF THE DEPENDENT VARIABLE TO PRINT            PRIN078

      NPR(NP)=NR                                                        PRIN079

      HEAD1(NP)=HEADR1(NN)                                              PRIN080

      HEAD2(NP)=HEADR2(NN)                                              PRIN081

      UNIT(NP)=UNITS(NN)                                                PRIN082

    7 CONTINUE                                                          PRIN083

      NP1=MIN0(NP,3)                                                    PRIN084

      PDEV=ABSORB=DOPP=0.                                               PRIN085

C********* PRINT COLUMN HEADINGS AT THE BEGINNING OF EACH RAY           PRIN086

   10 IF (IHOP.NE.0) GO TO 12                                           PRIN087

      PRINT 1100, (HEAD1(NN),HEAD2(NN),NN=1,NP1)                        PRIN088

 1100 FORMAT (44X,7HAZIMUTH/43X,9HDEVIATION,8X,9HELEVATION/             PRIN089

     1 19X,16HHEIGHT     RANGE,1X,2(5X,12HXMTR   LOCAL),5X,26HPOLARIZATIPRIN090

     2ON    GROUP PATH,5A6,A5)                                          PRIN091

      PRINT 1150, (UNIT(NN),NN=1,NP1)                                   PRIN092

 1150 FORMAT (13X,2(8X,2HKM),2X,2(6X,3HDEG,5X,3HDEG),6X,12HREAL    IMAG,PRIN093

     1 7X,2HKM,4X,3(4X,A6,2X))                                          PRIN094

      IF (RAYSET.EQ.0.) GO TO 12                                        PRIN095

C********* PUNCH A TRANSMITTER RAYSET                                   PRIN096

      TLOND=TLON*DEGS                                                   PRIN097

      IF (TLOND.LT.0.) TLOND=TLOND+360.                                 PRIN098

      TLATD=TLAT*DEGS                                                   PRIN099

      IF (TLATD.LT.0.) TLATD=TLATD+360.                                 PRIN100

      AZ=AZ1*DEGS                                                       PRIN101

      EL=BETA*DEGS                                                      PRIN102

      NHOP=HOP                                                          PRIN103

      PUNCH 1200, ID(1),TYPE(NTYP),XMTRH,TLATD,TLOND,RCVRH,F,AZ,EL,POLARPRIN104

     1,NHOP,1HT                                                         PRIN105

 1200 FORMAT (A3,A1,4PF9.0,3P2F6.0,4P2F9.0,5P2F10.0,5X,2P2F5.0,I1,A1)   PRIN106

C*********                                                              PRIN107

   12 V=0.                                                              PRIN108

      IF (N2.NE.0.) V=(R(4)**2+R(5)**2+R(6)**2)/N2-1.                   PRIN109

      H=R(1)-EARTHR                                                     PRIN110

      STH=SIN (THETA)                                                   PRIN111

      CTH=SIN (PID2-THETA)                                              PRIN112

C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER    PRIN113

      XP=R(1)*STH*SIN (PID2-PHI)-XR                                     PRIN114

      YP=R(1)*STH*SIN (PHI)-YR                                          PRIN115

      ZP=R(1)*CTH-ZR                                                    PRIN116

C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER ANDPRIN117

C********* ROTATED                                                      PRIN118

      EPS=XP*G1(1,1)+YP*G1(1,2)+ZP*G1(1,3)                              PRIN119

      ETA=XP*G1(2,1)+YP*G1(2,2)+ZP*G1(2,3)                              PRIN120

      ZETA=XP*G1(3,1)+YP*G1(3,2)+ZP*G1(3,3)                             PRIN121

      RCE2=ETA**2+ZETA**2                                               PRIN122

      RCE=SQRT (RCE2)                                                   PRIN123

C********* GROUND RANGE                                                 PRIN124

      RANGE=EARTHR*ATAN2(RCE,EARTHR+EPS+XMTRH)                          PRIN125 

C********* ANGLE OF WAVE NORMAL WITH LOCAL HORIZONTAL                   PRIN126

      ELL=ATAN2(R(4),SQRT (R(5)**2+R(6)**2))*DEGS                       PRIN127

C********* STRAIGHT LINE DISTANCE FROM TRANSMITTER TO RAY POINT         PRIN128

      SR=SQRT (RCE2+EPS**2)                                             PRIN129

      IF (NP.LT.1) GO TO 16                                             PRIN130

      DO 15 I=1,NP                                                      PRIN131

      NN=NPR(I)                                                         PRIN132

   15 RPRINT(I)=R(NN)                                                   PRIN133

   16 IF (SR.GE.1.E-6) GO TO 20                                         PRIN134

C********* TOO CLOSE TO TRANSMITTER TO CALCULATE DIRECTION FROM         PRIN135

C********* TRANSMITTER                                                  PRIN136

      PRINT 1500, V,NWHY,H,RANGE,ELL,POLAR,T,(RPRINT(NN),NN=1,NP1)      PRIN137

 1500 FORMAT (1X,E6.0,1X,A8,F10.4,F11.4,26X,F8.3,F9.3,F8.3,4F12.4)      PRIN138

      GO TO 40                                                          PRIN139

C********* ELEVATION ANGLE OF RAY POINT FROM TRANSMITTER                PRIN140

   20 EL=ATAN2(EPS,RCE)*DEGS                                            PRIN141

      IF (RCE.GE.1.E-6) GO TO 30                                        PRIN142

C********* NEARLY DIRECTLY ABOVE OR BELOW TRANSMITTER. CAN NOT CALCULATEPRIN143

C********* AZIMUTH DIRECTION FROM TRANSMITTER ACCURATELY                PRIN144

      PRINT 2500, V,NWHY,H,RANGE,EL,ELL,POLAR,T,(RPRINT(NN),NN=1,NP1)   PRIN145

 2500 FORMAT (1X,E6.0,1X,A8,F10.4,F11.4,17X,F9.3,F8.3,F9.3,F8.3,        PRIN146

     1 4F12.4)                                                          PRIN147

      GO TO 40                                                          PRIN148

C********* AZIMUTH ANGLE OF RAY POINT FROM TRANSMITTER                  PRIN149

   30 ANGA=ATAN2(ETA,ZETA)                                              PRIN150

      AZDEV=180.-AMOD(540.-(AZ1-ANGA)*DEGS,360.)                        PRIN151

      IF (R(5).NE.0..OR.R(6).NE.0.) GO TO 34                            PRIN152

C********* WAVE NORMAL IS VERTICAL, SO AZIMUTH DIRECTION CANNOT BE      PRIN153

C********* CALCULATED                                                   PRIN154

      PRINT 3000, V,NWHY,H,RANGE,AZDEV,EL,ELL,POLAR,T,(RPRINT(NN),NN=1, PRIN155

     1 NP1)                                                             PRIN156

 3000 FORMAT (1X,E6.0,1X,A8,F10.4,F11.4,F9.3,8X,F9.3,F8.3,F9.3,F8.3,    PRIN157

     1 4F12.4)                                                          PRIN158

      GO TO 40                                                          PRIN159

   34 ANA=ANGA-ALPH                                                     PRIN160

      SANA=SIN (ANA)                                                    PRIN161

      SPHI=SANA*STHR/STH                                                PRIN162

      CPHI=-SIN (PID2-ANA)*SIN (PID2-(PHI-PHIR))+SANA*SIN (PHI-PHIR)    PRIN163

     1 *CTHR                                                            PRIN164

      AZA=180.-AMOD (540.-(ATAN2(SPHI,CPHI)-ATAN2(R(6),R(5)))*DEGS,360.)PRIN165

      PRINT 3500, V,NWHY,H,RANGE,AZDEV,AZA,EL,ELL,POLAR,T,(RPRINT(NN),NNPRIN166

     1 =1,NP1)                                                          PRIN167

 3500 FORMAT (1X,E6.0,1X,A8,F10.4,F11.4,2(F9.3,F8.3),F9.3,F8.3,         PRIN168

     1 4F12.4)                                                          PRIN169

C*********                                                              PRIN170

   40 LINES=LINES+1                                                     PRIN171

      IF (NP.LE.3) GO TO 45                                             PRIN172

C********* ADDITIONAL LINE TO PRINT REMAINING DEPENDENT INTEGRATION     PRIN173

C********* VARIABLES                                                    PRIN174

      PRINT 4000, (RPRINT(NN),NN=4,NP)                                  PRIN175

 4000 FORMAT (99X,3F12.4)                                               PRIN176

      LINES=LINES+1                                                     PRIN177

   45 IF (CARD.EQ.0.) RETURN                                            PRIN178

C                                                                       PRIN179

C********* PUNCH A RAYSET                                               PRIN180

      IF (AZDEV.LT.-90.) AZDEV=AZDEV+360.                               PRIN181

      IF (AZA.LT.-90.) AZA=AZA+360.                                     PRIN182

      TDEV=T-SR                                                         PRIN183

      NR=6                                                              PRIN184

      IF (W(57).EQ.0.) GO TO 47                                         PRIN185

C********* PHASE PATH                                                   PRIN186

      NR=NR+1                                                           PRIN187

      PDEV=R(NR)-SR                                                     PRIN188

   47 IF (W(58).EQ.0.) GO TO 48                                         PRIN189

C********* ABSORPTION                                                   PRIN190

      NR=NR+1                                                           PRIN191

      ABSORB=R(NR)                                                      PRIN192

C********* DOPPLER SHIFT                                                PRIN193

   48 IF (W(59).NE.0.) DOPP=R(NR+1)                                     PRIN194

      PUNCH 4500, HPUNCH,RANGE,AZDEV,AZA,ELL,SR,TDEV,PDEV,ABSORB,DOPP,  PRIN195

     1 POLAR,IHOP,NWHY                                                  PRIN196

 4500 FORMAT (4P2F9.0,3P3F6.0,3PF8.0,3P4F6.0,2P2F5.0,I1,A1)             PRIN197

      RETURN                                                            PRIN198

         END                                                            PRIN199-

      SUBROUTINE RAYPLT                                                 YPLT001

C        REPLACES SUBROUTINES RAYPLT,PLOT, AND LABPLT IF PLOTS ARE      YPLT002

C        NOT WANTED OR IF A PLOTTER IS NOT AVAILABLE                    YPLT003

      COMMON /WW/ ID(10),W0,W(400)                                      YPLT004

      EQUIVALENCE (PLT,W(81))                                           YPLT005

      PLT=0.                                                            YPLT006

      ENTRY ENDPLT                                                      YPLT007

      RETURN                                                            YPLT008

         END                                                            YPLT  9-

      SUBROUTINE RAYPLT                                                 RAYP001

C        W(81)=1. PLOTS PROJECTION OF RAYPATH ON VERTICAL PLANE         RAYP002

C             =2. PLOTS PROJECTION OF RAYPATH ON GROUND                 RAYP003

      COMMON /PLT/ XL,XR,YB,YT,RESET                                    RAYP004

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                RAYP005

      COMMON /FLG/ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCH             RAYP006

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 RAYP007

      EQUIVALENCE (TH,R(2)),(PH,R(3))                                   RAYP008

      EQUIVALENCE (EARTHR,W(2)),(PLAT,W(24)),(PLON,W(25)),(PLT,W(81)),  RAYP009

     1 (FACTR,W(82)),(LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))RAYP010

     2,(TIC,W(87)),(HB,W(88))                                           RAYP011

      REAL LLAT,LLON,LTIC                                               RAYP012

      LOGICAL NEWWR,NEWWP,PENET                                         RAYP013

      IF (.NOT.NEWWR) GO TO 5                                           RAYP014

C                                                                       RAYP015

C          NEW W ARRAY -- REINITIALIZE                                  RAYP016

      NEWWR=.FALSE.                                                     RAYP017

      RESET=1.                                                          RAYP018

CONVERT COORDINATES OF VERTICAL PLANE FROM GEOGRAPHIC TO GEOMAGNETIC    RAYP019

      SW=SIN (PLAT)                                                     RAYP020

      CW=SIN (PID2-PLAT)                                                RAYP021

      SLM=SIN (LLAT)                                                    RAYP022

      CLM=SIN (PID2-LLAT)                                               RAYP023

      SRM=SIN (RLAT)                                                    RAYP024

      CRM=SIN (PID2-RLAT)                                               RAYP025

      CDPHI=SIN (PID2-(LLON-PLON))                                      RAYP026

      PHL=ATAN2(SIN (LLON-PLON)*CLM,CDPHI*SW*CLM-CW*SLM)                RAYP027

      CTHL=CDPHI*CW*CLM+SW*SLM                                          RAYP028

      STHL=SIN (ACOS (CTHL))                                            RAYP029

      CDPHI=SIN (PID2-(RLON-PLON))                                      RAYP030

      PHR=ATAN2(SIN (RLON-PLON)*CRM,CDPHI*SW*CRM-CW*SRM)                RAYP031

      CTHR=CDPHI*CW*CRM+SW*SRM                                          RAYP032

      STHR=SIN (ACOS (CTHR))                                            RAYP033

      CLR=CTHL*CTHR+STHL*STHR*SIN (PID2-(PHL-PHR))                      RAYP034

      SLR=SQRT (1.-CLR**2)                                              RAYP035

      IF (PLT.EQ.2.) GO TO 3                                            RAYP036

      FACTR=1.                                                          RAYP037

      R0=EARTHR+HB                                                      RAYP038

      ALPHA=.5*ACOS (CLR)                                               RAYP039

      XR=R0*SIN (ALPHA)                                                 RAYP040

      XL=-XR                                                            RAYP041

      YB=R0*SIN (PID2-ALPHA)                                            RAYP042

      YT=YB+2.*XR                                                       RAYP043

      GO TO 5                                                           RAYP044

    3 IF (FACTR.EQ.0.) FACTR=1.                                         RAYP045

      ALPH1=ATAN2(STHR*SIN (PHR-PHL),(CTHR-CTHL*CLR)/STHL)              RAYP046

      XL=0.                                                             RAYP047

      XR=EARTHR*ACOS (CLR)                                              RAYP048

      YT=0.5*XR/FACTR                                                   RAYP049

      YB=-YT                                                            RAYP050

C                                                                       RAYP051

    5 STH=SIN (TH)                                                      RAYP052

      CTH=SIN (PID2-TH)                                                 RAYP053

      CR=CTHR*CTH+STHR*STH*SIN (PID2-(PHR-PH))                          RAYP054

      CL=CTHL*CTH+STHL*STH*SIN (PID2-(PHL-PH))                          RAYP055

      CEA=ATAN2(CR-CL*CLR,CL*SLR)                                       RAYP056

      NEW=1                                                             RAYP057

      IF (IHOP.NE.0) NEW=0                                              RAYP058

      IF (PLT.EQ.2.) GO TO 10                                           RAYP059

      CALL PLOT (R(1)*SIN(CEA-ALPHA),R(1)*SIN(PID2-(CEA-ALPHA)),NEW)    RAYP060

      RETURN                                                            RAYP061

   10 SL=SQRT (1.-CL**2)                                                RAYP062

      TMP1=STH*SIN (PH-PHL)                                             RAYP063

      TMP2=(CTH-CTHL*CL)/STHL                                           RAYP064

      ALPH2=0.                                                          RAYP065

      IF (TMP1.NE.0..OR.TMP2.NE.0.) ALPH2=ATAN2(TMP1,TMP2)              RAYP066

      CALL PLOT (EARTHR*CEA,EARTHR*ASIN(SL*SIN (ALPH1-ALPH2)),NEW)      RAYP067

      RETURN                                                            RAYP068

C                                                                       RAYP069

C          DRAW AXES AND CALL FOR LABELING AND TERMINATION OF THIS PLOT RAYP070

      ENTRY ENDPLT                                                      RAYP071

      TICKX=0.01*(YT-YB)                                                RAYP072

      IF (PLT.EQ.2.) GO TO 25                                           RAYP073

      R1=EARTHR-TICKX                                                   RAYP074

      X=XL                                                              RAYP075

      Y=YB                                                              RAYP076

      CALL PLOT (X,Y,1)                                                 RAYP077

      NTIC=2                                                            RAYP078

      IF (TIC.NE.0.) NTIC=NTIC+2.*ALPHA/TIC                             RAYP079

      NLINE=MAX0 (1,100/NTIC)                                           RAYP080

      DO 20 I=1,NTIC                                                    RAYP081

      ANG=-ALPHA+(I-1)*TIC                                              RAYP082

      CALL PLOT (R1*SIN (ANG),R1*SIN (PID2-ANG),0)                      RAYP083

      CALL PLOT (X,Y,0)                                                 RAYP084

      DO 20 J=1,NLINE                                                   RAYP085

      ANG=ANG+TIC/NLINE                                                 RAYP086

      X=EARTHR*SIN (ANG)                                                RAYP087

      Y=EARTHR*SIN (PID2-ANG)                                           RAYP088

   20 CALL PLOT (X,Y,0)                                                 RAYP089

      CALL PLOT (XR,YB,0)                                               RAYP090

      GO TO 50                                                          RAYP091

   25 DTIC=TIC*EARTHR                                                   RAYP092

      LTIC=DTIC/FACTR                                                   RAYP093

      TICY=XL+0.01*(XR-XL)                                              RAYP094

      NTIC=YT/LTIC                                                      RAYP095

      TIC1=-LTIC*NTIC                                                   RAYP096

      CALL PLOT (XL,YB,1)                                               RAYP097

      NTIC=2*NTIC+1                                                     RAYP098

      DO 30 I=1,NTIC                                                    RAYP099

      Y=TIC1+(I-1)*LTIC                                                 RAYP100

      CALL PLOT (XL,Y,0)                                                RAYP101

      CALL PLOT (TICY,Y,0)                                              RAYP102

   30 CALL PLOT (XL,Y,0)                                                RAYP103

      CALL PLOT (XL,YT,0)                                               RAYP104

      CALL PLOT (XL,0.,1)                                               RAYP105

      NTIC=(XR-XL)/DTIC                                                 RAYP106

      DO 40 I=1,NTIC                                                    RAYP107

      X=I*DTIC                                                          RAYP108

      CALL PLOT (X,0.,0)                                                RAYP109

      CALL PLOT (X,TICKX,0)                                             RAYP110

   40 CALL PLOT (X,0.,0)                                                RAYP111

      CALL PLOT (XR,0.,0)                                               RAYP112

   50 CALL LABPLT                                                       RAYP113

      CALL PLTEND                                                       RAYP114

      RETURN                                                            RAYP115

         END                                                            RAYP116-

      SUBROUTINE PLOT (X,Y,NEW)                                         PLOT001

      COMMON /PLT/ XMIN0,XMAX0,YMIN0,YMAX0,RESET                        PLOT002

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY                            PLOT003

      DATA (INITAL=1),(MINX=0),(MINY=0),(MAXX=1023),(MAXY=1023),        PLOT004

     1 (MINX0=23),(MINY0=23),(MAXX0=1023),(MAXY0=1023)                  PLOT005

C                                                                       PLOT006

C          INITIALIZE LIBRARY PLOTTING ROUTINES                         PLOT007

      IF (INITAL.EQ.0) GO TO 1                                          PLOT008

      INITAL=0                                                          PLOT009

      CALL DDINIT (1,1H )                                               PLOT010

C                                                                       PLOT011

C          COMPUTE SCALE FACTORS                                        PLOT012

    1 IF (RESET.EQ.0.) GO TO 5                                          PLOT013

      RESET=0.                                                          PLOT014

      XSCALE=(MAXX0-MINX0)/(XMAX0-XMIN0)                                PLOT015

      YSCALE=(MAXY0-MINY0)/(YMAX0-YMIN0)                                PLOT016

      XMIN=XMIN0-(MINX0-MINX)/XSCALE                                    PLOT017

      YMIN=YMIN0-(MINY0-MINY)/YSCALE                                    PLOT018

      XMAX=XMAX0+(MAXX-MAXX0)/XSCALE                                    PLOT019

      YMAX=YMAX0+(MAXY-MAXY0)/YSCALE                                    PLOT020

C                                                                       PLOT021

C          START A NEW LINE                                             PLOT022

    5 IF (NEW.EQ.0) GO TO 10                                            PLOT023

      IX=MINX0+(X-XMIN0)*XSCALE                                         PLOT024

      IY=MINY0+(Y-YMIN0)*YSCALE                                         PLOT025

      IF (IX.GE.MINX.AND.IX.LE.MAXX.AND.IY.GE.MINY.AND.IY.LE.MAXY)      PLOT026

     1 CALL DDBP                                                        PLOT027

      GO TO 50                                                          PLOT028

C                                                                       PLOT029

C          HORIZONTAL DISPLACEMENT                                      PLOT030

   10 XS=X-XOLD                                                         PLOT031

      YS=Y-YOLD                                                         PLOT032

      IF (XS) 11,12,16                                                  PLOT033

C          NEGATIVE                                                     PLOT034

   11 X1=XMAX                                                           PLOT035

      X2=XMIN                                                           PLOT036

      GO TO 20                                                          PLOT037

C          ZERO                                                         PLOT038

   12 IF (YS) 13,50,14                                                  PLOT039

   13 S1=(YMAX-YOLD)/YS                                                 PLOT040

      S2=(YMIN-YOLD)/YS                                                 PLOT041

      GO TO 40                                                          PLOT042

   14 S1=(YMIN-YOLD)/YS                                                 PLOT043

      S2=(YMAX-YOLD)/YS                                                 PLOT044

      GO TO 40                                                          PLOT045

C          POSITIVE                                                     PLOT046

   16 X1=XMIN                                                           PLOT047

      X2=XMAX                                                           PLOT048

C                                                                       PLOT049

C          VERTICAL DISPLACEMENT                                        PLOT050

   20 IF (YS) 21,22,26                                                  PLOT051

C          NEGATIVE                                                     PLOT052

   21 Y1=YMAX                                                           PLOT053

      Y2=YMIN                                                           PLOT054

      GO TO 30                                                          PLOT055

C          ZERO                                                         PLOT056

   22 S1=(X1-XOLD)/XS                                                   PLOT057

      S2=(X2-XOLD)/XS                                                   PLOT058

      GO TO 40                                                          PLOT059

C          POSITIVE                                                     PLOT060

   26 Y1=YMIN                                                           PLOT061

      Y2=YMAX                                                           PLOT062

C                                                                       PLOT063

   30 S1=AMAX1((X1-XOLD)/XS,(Y1-YOLD)/YS)                               PLOT064

      S2=AMIN1((X2-XOLD)/XS,(Y2-YOLD)/YS)                               PLOT065

C                                                                       PLOT066

C          PLOT LINE -- CHECKING FOR BORDER CROSSINGS                   PLOT067

   40 S=SQRT(XS**2+YS**2)                                               PLOT068

      IF (S2.LT.0..OR.S*S1-S.GT.0.) GO TO 50                            PLOT069

      IF (S1.LT.0.) GO TO 42                                            PLOT070

C          PREVIOUS POINT OFF GRAPH                                     PLOT071

      IX=MINX0+(XOLD+XS*S1-XMIN0)*XSCALE+0.5                            PLOT072

      IY=MINY0+(YOLD+YS*S1-YMIN0)*YSCALE+0.5                            PLOT073

      CALL DDBP                                                         PLOT074

   42 IF (S*S2-S.GT.0.) GO TO 44                                        PLOT075

C          CURRENT POINT OFF GRAPH                                      PLOT076

      IX=MINX0+(XOLD+XS*S2-XMIN0)*XSCALE+0.5                            PLOT077

      IY=MINY0+(YOLD+YS*S2-YMIN0)*YSCALE+0.5                            PLOT078

      CALL DDVC                                                         PLOT079

      GO TO 50                                                          PLOT080

C          CURRENT POINT ON GRAPH                                       PLOT081

   44 IX=MINX0+(X-XMIN0)*XSCALE+0.5                                     PLOT082

      IY=MINY0+(Y-YMIN0)*YSCALE+0.5                                     PLOT083

      CALL DDVC                                                         PLOT084

C                                                                       PLOT085

C          EXIT ROUTINE                                                 PLOT086

   50 XOLD=X                                                            PLOT087

      YOLD=Y                                                            PLOT088

      RETURN                                                            PLOT089

C                                                                       PLOT090

C          TERMINATE THE CURRENT PLOT                                   PLOT091

      ENTRY PLTEND                                                      PLOT092

      CALL DDFR                                                         PLOT093

      RETURN                                                            PLOT094

         END                                                            PLOT 95-

      SUBROUTINE LABPLT                                                 LABP001

C          LABEL THE CURRENT PLOT                                       LABP002

      DIMENSION LABEL(9),TYPE(3)                                        LABP003

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY                            LABP004

      COMMON /CONST/ PI,PIT2,PID2,DEGS,DUM(4)                           LABP005

      COMMOÎ /FLG¯ NTYP,NEWWR,NEWWP,PENET,LINES,IHOP,HPUNCÈ             LABP006

      COMMON /WW/ ID(10),W0,W(400)                                      LABP007

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(AZ1,W(10)),(PLT,W(81)),       LABP008

     1 (FACTR,W(82)),(TIC,W(87))                                        LABP009

      LOGICAL NEWWR,NEWWP,PENET                                         LABP010

      REAL LTIC                                                         LABP011

      DATA (TYPE=8HEXTRAORD,8HNO FIELD,8HORDINARY)                      LABP012

      IOR=IT=0                                                          LABP013

      IS=2                                                              LABP014

      IX=0     $  IY=1023  $  CALL DDTAB  $  CALL DDTEXT (7,ID)         LABP015

      NDATE=IDATE(0)                                                    LABP016

      CALL DDTEXT (1,NDATE)                                             LABP017

      AZA=AZ1*DEGS                                                      LABP018

      DTIC=TIC*EARTHR                                                   LABP019

      ENCODE (72,1000,LABEL)  F,AZA,TYPE(NTYP),DTIC                     LABP020

 1000 FORMAT (3HF =,F7.3,6H, AZ =,F7.2,2H, ,A8,2H, ,F7.2,24H KM BETWEEN LABP021

     1TICK MARKS".)                                                     LABP022

      IX=0     $  IY=991   $  CALL DDTAB  $  CALL DDTEXT (9,LABEL)      LABP023

      IF (PLT.EQ.1.) RETURN                                             LABP024

      LTIC=DTIC/FACTR                                                   LABP025

      ENCODE (32,2000,LABEL)  LTIC                                      LABP026

 2000 FORMAT (F7.2,24H KM BETWEEN TICK MARKS".)                         LABP027

      IOR=1                                                             LABP028

      IX=0     $  IY=0     $  CALL DDTAB  $  CALL DDTEXT (4,LABEL)      LABP029

      IOR=0                                                             LABP030

      RETURN                                                            LABP031

         END                                                            LABP 32-

      SUBROUTINE RKAM                                                   RKAM001

C        NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS                RKAM002

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

      COMMON Y(20),T,STEP,DYDT(20)                                      RKAM004

      DIMENSION DELY(4,20),BET(4),XV(5),FV(4,20),YU(5,20)               RKAM005

      TYPE DOUBLE YU                                                    RKAM006

      IF (RSTART.EQ.0.) GO TO 1000                                      RKAM007

      LL=MM=1                                                           RKAM008

      IF (MODE.EQ.1) MM=4                                               RKAM009

      ALPHA=T                                                           RKAM010

      EPM=0.0                                                           RKAM011

      BET(1)=BET(2)=0.5                                                 RKAM012

      BET(3)=1.0                                                        RKAM013

      BET(4)=0.0                                                        RKAM014

      STEP=SPACE                                                        RKAM015

      R=19.0/270.0                                                      RKAM016

      XV(MM)=T                                                          RKAM017

      IF (E1MIN.LE.0.) E1MIN=E1MAX/55.                                  RKAM018

      IF (FACT.LE.0.) FACT=0.5                                          RKAM019

      CALL HAMLTN                                                       RKAM020

      DO 320 I=1,NN                                                     RKAM021

      FV(MM,I)=DYDT(I)                                                  RKAM022

  320 YU(MM,I)=Y(I)                                                     RKAM023

      RSTART=0.                                                         RKAM024

      GO TO 1001                                                        RKAM025

 1000 IF (MODE.NE.1) GO TO 2000                                         RKAM026

C                                                                       RKAM027

C          RUNGE-KUTTA                                                  RKAM028

 1001 DO 1034 K=1,4                                                     RKAM029

      DO 1350 I=1,NN                                                    RKAM030

      DELY(K,I)=STEP*FV(MM,I)                                           RKAM031

      Z=YU(MM,I)                                                        RKAM032

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

      T=BET(K)*STEP+XV(MM)                                              RKAM034

      CALL HAMLTN                                                       RKAM035

      DO 1034 I=1,NN                                                    RKAM036

 1034 FV(MM,I)=DYDT(I)                                                  RKAM037

      DO 1039 I=1,NN                                                    RKAM038

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

 1039 YU(MM+1,I)=YU(MM,I)+DEL                                           RKAM040

      MM=MM+1                                                           RKAM041

      XV(MM)=XV(MM-1)+STEP                                              RKAM042

      DO 1400 I=1,NN                                                    RKAM043

 1400 Y(I)=YU(MM,I)                                                     RKAM044

      T=XV(MM)                                                          RKAM045

      CALL HAMLTN                                                       RKAM046

      IF (MODE.EQ.1) GO TO 42                                           RKAM047

      DO 150 I=1,NN                                                     RKAM048

  150 FV(MM,I)=DYDT(I)                                                  RKAM049

      IF (MM.LE.3) GO TO 1001                                           RKAM050

C                                                                       RKAM051

C          ADAMS-MOULTON                                                RKAM052

 2000 DO 2048 I=1,NN                                                    RKAM053

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

      Y(I)=YU(4,I)+DEL                                                  RKAM055

 2048 DELY(1,I)=Y(I)                                                    RKAM056

      T=XV(4)+STEP                                                      RKAM057

      CALL HAMLTN                                                       RKAM058

      XV(5)=T                                                           RKAM059

      DO 2051 I=1,NN                                                    RKAM060

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

      YU(5,I)=YU(4,I)+DEL                                               RKAM062

 2051 Y(I)=YU(5,I)                                                      RKAM063

      CALL HAMLTN                                                       RKAM064

      IF (MODE.LE.2) GO TO 42                                           RKAM065

C                                                                       RKAM066

C          ERROR ANALYSIS                                               RKAM067

      SSE=0.0                                                           RKAM068

      DO 3033 I=1,NN                                                    RKAM069

      EPSIL=R*ABS(Y(I)-DELY(1,I))                                       RKAM070

      IF (MODE.EQ.3.AND.Y(I).NE.0.) EPSIL=EPSIL/ABS(Y(I))               RKAM071

      IF (SSE.LT.EPSIL) SSE=EPSIL                                       RKAM072

 3033 CONTINUE                                                          RKAM073

      IF (E1MAX.GT.SSE) GO TO 3035                                      RKAM074

      IF (ABS(STEP).LE.E2MIN) GO TO 42                                  RKAM075

      LL=MM=1                                                           RKAM076

      STEP=STEP*FACT                                                    RKAM077

      GO TO 1001                                                        RKAM078

 3035 IF (LL.LE.1.OR.SSE.GE.E1MIN.OR.E2MAX.LE.ABS(STEP)) GO TO 42       RKAM079

      LL=2                                                              RKAM080

      MM=3                                                              RKAM081

      XV(2)=XV(3)                                                       RKAM082

      XV(3)=XV(5)                                                       RKAM083

      DO 5363 I=1,NN                                                    RKAM084

      FV(2,I)=FV(3,I)                                                   RKAM085

      FV(3,I)=DYDT(I)                                                   RKAM086

      YU(2,I)=YU(3,I)                                                   RKAM087

 5363 YU(3,I)=YU(5,I)                                                   RKAM088

      STEP=2.0*STEP                                                     RKAM089

      GO TO 1001                                                        RKAM090

C                                                                       RKAM091

C          EXIT ROUTINE                                                 RKAM092

   42 LL=2                                                              RKAM093

      MM=4                                                              RKAM094

      DO 12 K=1,3                                                       RKAM095

      XV(K)=XV(K+1)                                                     RKAM096

      DO 12 I=1,NN                                                      RKAM097

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

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

      XV(4)=XV(5)                                                       RKAM100

      DO 52 I=1,NN                                                      RKAM101

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

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

      IF (MODE.LE.2) RETURN                                             RKAM104

      E=ABS(XV(4)-ALPHA)                                                RKAM105

      IF (E.LE.EPM) GO TO 2000                                          RKAM106

      EPM=E                                                             RKAM107

      RETURN                                                            RKAM108

         END                                                            RKAM109-

      SUBROUTINE HAMLTN                                                 HAML001

C********* CALCULATES HAMILTONS EQUATIONS FOR RAY TRACING               HAML002

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN                   HAML003

      COMMON /RIN/ MODRIN(3),COLL,FIELD,SPACE,KAY2,KAY2I,               HAML004

     1              H,HI,PHPT,PHPTI,PHPR,PHPRI,PHPTH,PHPTHI,PHPPH,PHPPHIHAML005

     2, PHPOM,PHPOMI,PHPKR,PHPKRI,PHPKTH,PHPKTI, PHPKPH,PHPKPI          HAML006

     3 ,KPHPK,KPHPKI,POLAR,POLARI,LPOLAR,LPOLRI                         HAML007

      COMMON R(20),T,STP,DRDT(20) /WW/ ID(10),W0,W(400)                 HAML008

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

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

     2 (DKPHDT,DRDT(6)),(F,W(6))                                        HAML011

      REAL KR,KTH,KPH,KPHPK,KPHPKI,LPOLAR,LPOLRI,LOGTEN,K,KAY2,KAY2I    HAML012

      OM=PIT2*1.E6*F                                                    HAML013

      STH=SIN(TH)                                                       HAML014

      CTH=SIN(PID2-TH)                                                  HAML015

      RSTH=R(1)*STH                                                     HAML016

      RCTH=R(1)*CTH                                                     HAML017

      CALL RINDEX                                                       HAML018

      DRDT=-PHPKR/(PHPOM*C)                                             HAML019

      DTHDT=-PHPKTH/(PHPOM*R(1)*C)                                      HAML020

      DPHDT=-PHPKPH/(PHPOM*RSTH*C)                                      HAML021

      DKRDT=PHPR/(PHPOM*C)+KTH*DTHDT+KPH*STH*DPHDT                      HAML022

      DKTHDT=(PHPTH/(PHPOM*C)-KTH*DRDT+KPH*RCTH*DPHDT)/R(1)             HAML023

      DKPHDT=(PHPPH/(PHPOM*C)-KPH*STH*DRDT-KPH*RCTH*DTHDT)/RSTH         HAML024

      NR=6                                                              HAML025

C********* PHASE PATH                                                   HAML026

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

      NR=NR+1                                                           HAML028

      DRDT(NR)=-     KPHPK/PHPOM/OM                                     HAML029

C********* ABSORPTION                                                   HAML030

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

      NR=NR+1                                                           HAML032

      DRDT(NR)= 10./LOGTEN*KPHPK*KAY2I/(KR*KR+KTH*KTH+KPH*KPH)/PHPOM/C  HAML033

C********* DOPPLER SHIFT                                                HAML034

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

      NR=NR+1                                                           HAML036

      DRDT(NR)=-PHPT/PHPOM/C/PIT2                                       HAML037

C********* GEOMETRICAL PATH LENGTH                                      HAML038

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

      NR=NR+1                                                           HAML040

      DRDT(NR)=-SQRT(PHPKR**2+PHPKTH**2+PHPKPH**2)/PHPOM   /C           HAML041

C********* OTHER CALCULATIONS                                           HAML042

   25 CONTINUE                                                          HAML043

      RETURN                                                            HAML044

         END                                                            HAML045-

