	PROGRAM NITIAL                                                    

C		SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE

	CHARACTER*8 MFLD(2)			       	                 

	DOUBLE PRECISION PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN,W(400),HPUNCH

     1     ,N2,N2I,NTYP,R0,SP,CP,SDPH,CDPH,SL,CL,ALPHA,TH0,PH0,AZA,

     2     GAMMA,SGAMMA,CGAMMA,EL,CBETA,FN

	COMMON /LPLT/ NDATE

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

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

	CHARACTER*8 MODRIN(3),NDATE

	COMMON /RINC/ MODRIN

C	  ** SPACE IS LOGICAL, POLAR AND LPOLAR ARE COMPLEX **

	COMMON /RIN/ COLL,FIELD,SPACE,N2,N2I,PNP(10),POLAR,LPOLAR,

     1         SGN

	COMPLEX*16 PNP,POLAR,LPOLAR                                          

	LOGICAL SPACE,NEWWR,NEWWP,PENET                                   

	DOUBLE PRECISION STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT

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

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	CHARACTER*6 MODY

	COMMON /YYC/ MODY

	CHARACTER*6 MODZ

	COMMON /ZZC/ MODZ

	DOUBLE PRECISION Y(16),Z(4),R(20),T,STP,DRDT(20)

	COMMON /YY/ Y /ZZ/ Z                         

	CHARACTER*8 ID(10)

	COMMON R,T,STP,DRDT /WW/ ID,W0,W                 

	DOUBLE PRECISION RAY,EARTHR,XMTRH,TLAT,TLON,F,FBEG,FEND,FSTEP,

     1   AZ1,AZBEG,AZEND,AZSTEP,BETA,ELBEG,ELEND,ELSTEP,ONLY,HOP,

     2   MAXSTP,PLAT,PLON,INTYP,MAXERR,ERATIO,STEP1,STPMAX,STPMIN,

     3   FACTR,SKIP,RAYSET,PLT,PERT,MU

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

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

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

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

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

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

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

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

	CHARACTER*4 KOLL

	CHARACTER*6 MODSAV

C********* OPEN PLOT AND PUNCH FILES

C	OPEN(UNIT=5,FILE='PRN')

	OPEN(UNIT=4,FILE='PRINT.DAT',STATUS='NEW')

	OPEN(UNIT=2,FILE='PUNCH.DAT',STATUS='NEW')

	OPEN(UNIT=1,FILE='TEST.DAT',IOSTAT=NOK,STATUS='OLD')

	IF (NOK.NE.0) THEN

	  WRITE(*,*) 'FILE TEST.DAT NOT FOUND.....EXITING PROGRAM'

	  CALL EXIT

	ENDIF

C*********

C		GET DATE FROM USER

	PRINT 1

1	FORMAT (1X,'ENTER DATE AS MM/DD/YY:')

	READ (*,2) NDATE

2	FORMAT (A8)

C*********

C		SET UP PLOTTER

	WRITE (*,*) 'Turn Plotter ON.'

	PAUSE 'Press RETURN to Continue...'

	OPEN(UNIT=3,FILE='COM2',STATUS='OLD')

	WRITE(3,*) 'IN;IP;SC0,1023,0,1023;SP1;'

C*********

C	NDATE=IDATE(0)                                                    

C	SECOND=KLOCK(0)*.001                                              

	KOLL=4H  NO                                                       

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

C***********CONSTANTS                                                    

	PI=3.1415926536D0                                                 

	PIT2=2.0D0*PI                    					 

	PID2=PI/2.0D0                                         		 

	DEGS=180.0D0/PI	

	RAD=PI/180.0D0						 

	C=2.997925D5					 

	K=2.81785D-15*C**2/PI						 

	LOGTEN=DLOG(10.0D0)

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

	DO 5 NW=1,400

5	W(NW)=0.0D0

	PLON=0.0D0

	PLAT=PID2

	EARTHR=6370.0D0

	INTYP=3.0D0

	MAXERR=1.0D-4

	ERATIO=50.0D0

	STEP1=1.0D0

	STPMAX=100.0D0

	STPMIN=1.0D-8

	FACTR=0.5D0

	MAXSTP=1000.0D0

	HOP=1.0D0

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

10	CALL READ W

	F=0.0D0

	BETA=0.0D0

	AZ1=0.0D0

	IF (SKIP.EQ.0.0D0) SKIP=MAXSTP					 

12	RAY=DSIGN(1.0D0,RAY)					        	 

	NTYP=2.0D0+FIELD*RAY					         

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

13	MFLD(1)=8HEXTRAORD

	MFLD(2)=8HINARY

	GO TO 16

14	MFLD(1)=8HNO FIELD

	MFLD(2)=8H

	GO TO 16

15	MFLD(1)=8HORDINARY

	MFLD(2)=8H

16	MODSAV=MODX(2)

	IF (PERT.EQ.0.0D0) MODX(2)=6H

C********* PRINTS CHAPX,WAVE,DIPOLY EXPZ2 WITH 7 ELEMENTS EACH ON A

C********* VERTICAL SCREEN

C********* PRINTS INITIAL VALUES FOR THE W ARRAY PORTRAIT

C********* PRNTPOR SETS HP/LJ PRINTER TO PORTRAIT

C	CALL PRNTPOR

C	WRITE (5,*) CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),

C     1     CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),

C     1     CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),

C     1     CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),

C     1     CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),CHAR (),

	IF (RAYSET.NE.0.0D0) THEN

C	   PRINT 2000, ID,MODX(1),(W(NW),NW=101,107),      

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

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

	   WRITE (4,2000) ID,MODX(1),(W(NW),NW=101,107),      

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

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

	ENDIF

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

C********* PRINTS INITIAL VALUES FOR THE W ARRAY LANDSCAPE

C********* PRNTLAN SETS HP/LJ PRINTER TO LANDSCAPE

C	CALL PRNTLAN

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

	WRITE (4,1000) ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL              

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

     1  11H COLLISIONS/)					         

C	PRINT 1050

	WRITE (4,1050)

1050	FORMAT (' INITIAL VALUES FOR THE W ARRAY -- ALL ANGLES IN ',

     1    'RADIANS, ONLY NONZERO VALUES PRINTED'/)                        	         

	DO 17 NW=1,400

	   IF (W(NW).NE.0.0D0) THEN

C	       PRINT 1700, NW,W(NW)

	       WRITE (4,1700) NW,W(NW)

	   ENDIF			         

1700	   FORMAT (I4,1PE19.11,0P)

17	CONTINUE								 

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

	NEWWP=.TRUE.							 

	NEWWR=.TRUE.        						 

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

	N=6								 

	DO 20 NR=7,20							 

	  IF (W(50+NR).NE.0.0D0) N=N+1					 

20	CONTINUE								 

	MODE=INTYP							 

	STEP=STEP1							 

	E1MAX=MAXERR							 

	E1MIN=MAXERR/ERATIO                   				 

	E2MAX=STPMAX						         

	E2MIN=STPMIN							 

	FACT=FACTR							 

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

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

	R0=EARTHR+XMTRH							 

	SP=DSIN (PLAT)							 

	CP=DSIN (PID2-PLAT)						 

	SDPH=DSIN (TLON-PLON)						 

	CDPH=DSIN (PID2-(TLON-PLON))       				 

	SL=DSIN (TLAT)							 

	CL=DSIN (PID2-TLAT)						 

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

	TH0=DACOS (CDPH*CP*CL+SP*SL)					 

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

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

	NFREQ=1								 

	IF (FSTEP.NE.0.0D0) NFREQ=IDINT((FEND-FBEG)/FSTEP+1.5D0)			 

	NAZ=1

	IF (AZSTEP.NE.0.0D0) NAZ=IDINT((AZEND-AZBEG)/AZSTEP+1.5D0)

	NBETA=1

	IF (ELSTEP.NE.0.0D0) NBETA=IDINT((ELEND-ELBEG)/ELSTEP+1.5D0)			 

	DO 50 NF=1,NFREQ							 

	   F=FBEG+(NF-1)*FSTEP						 

	   DO 45 J=1,NAZ							 

	      AZ1=AZBEG+(J-1)*AZSTEP						 

	      AZA=AZ1*DEGS							 

	      GAMMA=PI-AZ1+ALPHA						 

	      SGAMMA=DSIN (GAMMA)						 

	      CGAMMA=DSIN (PID2-GAMMA)						 

	      DO 40 I=1,NBETA

		 BETA=ELBEG+(I-1)*ELSTEP

	         EL=BETA*DEGS

	         CBETA=DSIN (PID2-BETA)

	         R(1)=R0

	         R(2)=TH0								 

	         R(3)=PH0								 

	         R(4)=DSIN (BETA)							 

	         R(5)=CBETA*CGAMMA						 

	         R(6)=CBETA*SGAMMA					         

	         T=0.0D0								 

	         RSTART=1.							 

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

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

		 CALL RINDEX							 

	         IF (I.NE.1.AND.NPAGE.LT.3.AND.LINES.LE.17)

     1              GO TO 25

		 NPAGE=0

		 LINES=0							 

C

C************			PRINTING IN INCREMENTS OF 15 DEGREES

C

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

		 WRITE (4,1000)ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL

C	         PRINT 2400, F,AZA		

		 WRITE (4,2400) F,AZA				 

2400		 FORMAT (18X,11HFREQUENCY =,F12.6,

     1             37H MHZ, AZIMUTH ANGLE OF TRANSMISSION =,F12.6,

     2             4H DEG)						 

25		 NPAGE=NPAGE+1

	         PRINT 2501, EL						 	 

		 WRITE (4,2500) EL

2500		 FORMAT (/31X,33HELEVATION ANGLE OF TRANSMISSION =,

     1             F12.6,4H DEG/)  

2501		 FORMAT (/1X,33HELEVATION ANGLE OF TRANSMISSION =,

     1             F12.6,8H DEGREES/)

	         IF (N2.GT.0.0D0) GO TO 30		                                 

	         CALL ELECTX							 

	         FN=DSIGN(DSQRT(DABS(X))*F,X)					 

C	         PRINT 2900,FN						         

		 WRITE (4,2900) FN

2900		 FORMAT ('0','TRANSMITTER IN EVANESCENT REGION, ',

     1   'TRANSMISSION IMPOSSIBLE'/,'0','PLASMA FREQUENCY = ',E17.10)

	         GO TO 44

30		 MU=DSQRT (DABS(N2/(R(4)**2+R(5)**2+R(6)**2)))

	         DO 34 NN=4,6

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

	         DO 35 NN=7,N

35		 R(NN)=0.0D0

		 CALL TRACE

C	         OSEC=SECOND							 

C	         SECOND=KLOCK(0)*.001						 

C	         DIFF=SECOND-OSEC							 

C		 PRINT 3500, DIFF							 

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

	         IF (PENET.AND.ONLY.NE.0.0D0.AND.IHOP.EQ.1)

     1             GO TO 44

40	      CONTINUE

44	      IF (PLT.NE.0.0D0) CALL ENDPLT

45	   CONTINUE

	   IF (PENET.AND.ONLY.NE.0.0D0.AND.IHOP.EQ.1.AND.

     1        NAZ.EQ.1.AND.NBETA.EQ.1) GO TO 55

50	CONTINUE

55	IF (RAYSET.NE.0.0D0) THEN

C	   PRINT 5000

	   WRITE (4,5000)

	ENDIF

5000	FORMAT (78X,1H-)

	MODX(2)=MODSAV

	GO TO 10

	END

C

C==========================================================================

C	END OF MAIN PROGRAM

C==========================================================================

C

C######################################################################

C

C UNIQUE "BLOCK DATA SETUP" FOR THIS CASE -- NEEDS TO BE MODIFIED

C  FOR EACH INDIVIDUAL LINKING WITH THE SPECIFIC REFRACTIVE INDEX

C  EQUATIONS AND CALLS



	BLOCK DATA SETUP

	CHARACTER*8 MODRIN(3)

	CHARACTER*6 MODX(2),MODY,MODZ

	COMMON /RINC/ MODRIN

	COMMON /XXC/ MODX

	COMMON /YYC/ MODY

	COMMON /ZZC/ MODZ

	COMPLEX*16 N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT,

     1    POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2,D2,

     2    KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,

     3    KPHPK

	COMMON /RIN/ COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,

     1     PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,SGN

	LOGICAL SPACE

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX

	DOUBLE PRECISION Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,

     1    PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP

	COMMON /YY/ Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,

     1    PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP

	DOUBLE PRECISION Z,PZPR,PZPTH,PZPPH

	COMMON /ZZ/ Z,PZPR,PZPTH,PZPPH

	DATA X,PXPR,PXPTH,PXPPH,PXPT,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,

     1       PYRPP,YTH,PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP,Z,PZPR,

     2       PZPTH,PZPPH /25*0.0D0/,

     3	     MODRIN /8HAPPLETON,8H-HARTREE,8H FORMULA/,COLL /1./,

     4       FIELD /1./,

     5       MODX /6H CHAPX,6H WAVE /,MODY /6HDIPOLY/,MODZ /6H EXPZ2/

	END

C

C  

C#####################################################################



	SUBROUTINE EXIT

C

C	GRACE EXIT OF PROGRAM

C

	WRITE(*,*) 'EXITING PROGRAM'

C********* PRNTPOR SET HP/LJ PRINTER TO PORTRAIT

C	CALL PRNTPOR

C

	CLOSE(1)

	CLOSE(2)

	WRITE (3,*) 'SP0;'

	CLOSE(3)

	CLOSE(4)

C	CLOSE(5)

	STOP

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE READ W

C                          READS W ARRAY                               

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

COL 18     DEGREES IN RADIANS                                            

COL 19     GREAT CIRCLE DISTANCE IN KM TO RADIANS                        

COL 20     NAUTICAL MILES TO KM                                          

COL 21     FEET TO KM							 

C								         

	DOUBLE PRECISION PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN,W(400)            

	COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,DUM                       

	CHARACTER*8 ID(10)

	COMMON /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,DUM(3)

	EQUIVALENCE (EARTHR,W(2))				         

	INTEGER DEG,FEET							 



	WRITE (*,*) 'INSERT Plotter Paper.'

	PAUSE 'Press RETURN to Continue...'

	READ (1,1000,END=30) ID

1000	FORMAT (10A8)

4	READ (1,1100,END=20) NW,W(NW),DEG,KM,NM,FEET				 

1100	FORMAT (I3,D14.7,5I1)

	IF (NW.EQ.0) GO TO 40

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

	WRITE(*,4000) NW

4000	FORMAT(1X,'THE SUBSCRIPT ',I3,' ON THE W-ARRAY INPUT IS OUT 

     1OF BOUNDS. ALLOWABLE VALUES ARE 1 THROUGH 400.')			 

	CALL EXIT							 

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

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

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

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

	GO TO 4								 

20	WRITE(*,*) 'INCOMPLETE W ARRAY READ...EXITING'

	CALL EXIT

30	WRITE(*,*) 'SUCCESSFUL READ....PROGRAM COMPLETED'

	CALL EXIT

40	RETURN

	END								 

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

 	SUBROUTINE WAVE

C    A PERTURBATION TO AN ALPHA-CHAPMAN ELECTRON DENSITY MODEL

	DOUBLE PRECISION PI,PIT2,PID2,W(400),DUM(5),H,SINW,COSW,CONS

	COMMON /CONST/ PI,PIT2,PID2,DUM

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX,R(6),TMP,EXPQ                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX

	CHARACTER*8 ID(10)

	COMMON R /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,Z0,SH,DELTA,VSUBX,LAMBDX,LAMBDZ,TP

	EQUIVALENCE (EARTHR,W(2)),(Z0,W(151)),(SH,W(152)),(DELTA,W(153)

     1   ),(VSUBX,W(154)),(LAMBDX,W(155)),(LAMBDZ,W(156)),(TP,W(157))

	ENTRY ELECT1

	IF (X.EQ.0.0D0.AND.PXPR.EQ.0.0D0.AND.PXPTH.EQ.0.0D0.AND.

     1	  PXPPH.EQ.0.0D0) RETURN

	IF (DELTA.EQ.0.0D0.OR.SH.EQ.0.0D0) RETURN

	H=R(1)-EARTHR

	EXPQ=DEXP(-((H-Z0)/SH)**2)	

	TMP=PIT2*(TP+(PID2-R(2))*EARTHR/LAMBDX+H/LAMBDZ)

	SINW=DSIN(TMP)

	COSW=DSIN(PID2-TMP)

	CONS=1.0D0+DELTA*EXPQ*COSW

	IF (H.NE.0.0D0) PXPR=PXPR*CONS-X*DELTA*EXPQ*(2.0D0/SH**2*(H-Z0)

     1              *COSW+PIT2/LAMBDZ*SINW)

	PXPTH=PXPTH*CONS+X*DELTA*PIT2*EARTHR/LAMBDX*SINW*EXPQ

	PXPPH=PXPPH*CONS

	PXPT=0.0D0

	IF (VSUBX.NE.0.0D0) PXPT=-PIT2*VSUBX/LAMBDX*X*DELTA*EXPQ*SINW

	X=X*CONS

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE EXPZ2

C		COLLISION FREQUENCY PROFILE FROM TWO EXPONENTIALS

	DOUBLE PRECISION PI,PIT2,PID2,W(400),DUM(5)

	COMMON /CONST/ PI,PIT2,PID2,DUM

	CHARACTER*6 MODZ

	COMMON /ZZC/ MODZ

	DOUBLE PRECISION Z,PZPR,PZPTH,PZPPH,R(6),H,EXP1,EXP2

	COMMON /ZZ/ Z,PZPR,PZPTH,PZPPH

	CHARACTER*8 ID(10)

	COMMON R /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,F,NU1,H1,A1,NU2,H2,A2

	EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(NU1,W(251)),(H1,W(252))

     1              ,(A1,W(253)),(NU2,W(254)),(H2,W(255)),(A2,W(256)) 

	ENTRY COLFRZ

	H=R(1)-EARTHR

	EXP1=NU1*DEXP(-A1*(H-H1))

	EXP2=NU2*DEXP(-A2*(H-H2))

	Z=(EXP1+EXP2)/(PIT2*F*1.0D6)

	PZPR=(-A1*EXP1-A2*EXP2)/(PIT2*F*1.0D6)

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE DIPOLY

	DOUBLE PRECISION PI,PIT2,PID2,W(400),DUM(5)

	COMMON /CONST/ PI, PIT2,PID2,DUM

	CHARACTER*6 MODY

	COMMON /YYC/ MODY

	DOUBLE PRECISION Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,

     1    PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP,R(6)

	COMMON /YY/ Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH

     1              ,PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP

	CHARACTER*8 ID(10)

	COMMON R /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,F,FH,SINTH,COSTH,TERM9,T1

	EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FH,W(201))

	ENTRY MAGY

	SINTH=DSIN(R(2))

	COSTH=DSIN(PID2-R(2))

	TERM9=DSQRT(DABS(1.0D0+3.0D0*COSTH**2))

	T1=FH*(EARTHR/R(1))**3/F

	Y=T1*TERM9

	YR=2.0D0*T1*COSTH

	YTH=T1*SINTH

	PYRPR=-3.0D0*YR/R(1)

	PYRPT=-2.0D0*YTH

	PYTPR=-3.0D0*YTH/R(1)

	PYTPT=.5D0*YR

	PYPR=-3.0D0*Y/R(1)

	PYPTH=-3.0D0*Y*SINTH*COSTH/TERM9**2

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE CHAPX

C		CHAPMAN LAYER WITH TILTS, RIPPLES, AND GRADIENTS

C		W(104)=0.5 FOR AN ALPHA-CHAPMAN LAYER

C		      =1.0 FOR A BETA-CHAPMAN LAYER

	DOUBLE PRECISION PI,PIT2,PID2,W(400),DUM(5)

	COMMON /CONST/ PI, PIT2,PID2,DUM

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX,R(6)

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX

	CHARACTER*8 ID(10)

	COMMON R /WW/ ID,W0,W

	EQUIVALENCE (THETA,R(2))

	DOUBLE PRECISION EARTHR,F,FC,HM,SH,ALPHA,A,B,C,E,PERT,THETA

	EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC,W(101)),(HM,W(102)),

     1   (SH,W(103)),(ALPHA,W(104)),(A,W(105)),(B,W(106)),(C,W(107)),

     2   (E,W(108)),(PERT,W(150))

	DOUBLE PRECISION EXZ,Z,H,D,THETA2,TEMP

	ENTRY ELECTX

	THETA2=THETA-PID2

	HMAX=HM+EARTHR*E*THETA2

	H=R(1)-EARTHR

	Z=(H-HMAX)/SH

	D=0.0D0

	IF (B.NE.0.0D0) D=PIT2/B

	TEMP=1.0D0+A*DSIN(D*THETA2)+C*THETA2

	EXZ=1.0D0-DEXP(-Z)

	X=(FC/F)**2*TEMP*DEXP(ALPHA*(EXZ-Z))

	PXPR=-ALPHA*X*EXZ/SH

	PXPTH=X*(D*A*DSIN(PID2-D*THETA2)+C)/TEMP-PXPR*EARTHR*E

	PXPPH=0.0D0

	IF (PERT.NE.0.0D0) CALL ELECT1

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE AHWFWC

C		CALCULATES THE REFRACTIVE INDEX AND ITS GRADIENT USING THE

C		APPLETON-HARTREE FORMULA WITH FIELD, WITH COLLISIONS

	DOUBLE PRECISION PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN,W(400)

	COMMON /CONST/ PI,PIT2,PID2,DEGS,RADIAN,K,C,LOGTEN

	CHARACTER*8 MODRIN(3)

	COMMON /RINC/ MODRIN

	COMPLEX*16 N2,PNPR,PNPTH,PNPPH,PNPVR,PNPVTH,PNPVPH,NNP,PNPT,

     1    POLAR,LPOLAR,I,U,RAD,D,PNPPS,PNPX,PNPY,PNPZ,UX,UX2,D2,

     2    KAY2,H,PHPT,PHPR,PHPTH,PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,

     3    KPHPK,GAM

C      **** ASSIGNMENT OF COMPLEX VARIABLES AFFECTS THE /RIN/ COMMON

C        *****    BLOCK STATEMENT (30 BYTES) BELOW	  *****

C      ********* ALLOCATIONS ARE TO TWO BYTES *********

C	 ***** SPACE IS LOGICAL  ALL ELSE ABOVE ******

	COMMON /RIN/ COLL,FIELD,SPACE,KAY2,H,PHPT,PHPR,PHPTH,

     1     PHPPH,PHPOM,PHPKR,PHPKTH,PHPKPH,KPHPK,POLAR,LPOLAR,SGN

	LOGICAL SPACE

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX

	CHARACTER*6 MODY

	COMMON /YYC/ MODY

	DOUBLE PRECISION Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,

     1    PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP

	COMMON /YY/ Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,

     1    PYTPR,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP

	CHARACTER*6 MODZ

	COMMON /ZZC/ MODZ

	DOUBLE PRECISION Z,PZPR,PZPTH,PZPPH

	COMMON /ZZ/ Z,PZPR,PZPTH,PZPPH

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R,TH,PH,KR,KTH,KPH,K2

	COMMON R,TH,PH,KR,KTH,KPH /WW/ ID,W0,W

	DOUBLE PRECISION STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT

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

	DOUBLE PRECISION RAY,F,OM,C2,OM2,VR,VTH,VPH,V2,VDOTY,YLV,YL2,

     1     YT2,YT4,PPSPR,PPSPTH,PPSPPH,SCALE,ABSLIM

	EQUIVALENCE (RAY,W(1)),(F,W(6))

	DATA I /(0.0D0,1.0D0)/,ABSLIM /1.D-5/

	ENTRY RINDEX

	OM=PIT2*1.D6*F

	C2=C*C

	K2=KR*KR+KTH*KTH+KPH*KPH

	OM2=OM*OM

	VR=C/OM*KR

	VTH=C/OM*KTH

	VPH=C/OM*KPH

	CALL ELECTX

	CALL MAGY

	V2=VR**2+VTH**2+VPH**2

	VDOTY=VR*YR+VTH*YTH+VPH*YPH

	YLV=VDOTY/V2

	YL2=VDOTY**2/V2

	YT2=Y**2-YL2

	YT4=YT2*YT2

	CALL COLFRZ

	U=DCMPLX(1.0D0,-Z)

	UX=U-X

	UX2=UX*UX

	RAD=RAY*CDSQRT(YT4+4.0D0*YL2*UX2)

	D=2.0D0*U*UX-YT2+RAD

	D2=D*D

	N2=1.0D0-2.0D0*X*UX/D

	PNPPS=2.*X*UX*(-1.+(YT2-2.*UX2)/RAD)/D2

	PPSPR=YL2/Y*PYPR-(VR*PYRPR+VTH*PYTPR+VPH*PYPPR)*YLV

	PPSPTH=YL2/Y*PYPTH-(VR*PYRPT+VTH*PYTPT+VPH*PYPPT)*YLV

	PPSPPH=YL2/Y*PYPPH-(VR*PYRPP+VTH*PYTPP+VPH*PYPPP)*YLV

	PNPX=-(2.0D0*U*UX2-YT2*(U-2.0D0*X)+(YT4*(U-2.0D0*X)+4.0D0*YL2*

     1     UX*UX2)/RAD)/D2

	PNPY=2.0D0*X*UX*(-YT2+(YT4+2.0D0*YL2*UX2)/RAD)/(D2*Y)

	PNPZ=I*X*(-2.0D0*UX2-YT2+YT4/RAD)/D2

	PNPR=PNPX*PXPR+PNPY*PYPR+PNPZ*PZPR+PNPPS*PPSPR

	PNPTH=PNPX*PXPTH+PNPY*PYPTH+PNPZ*PZPTH+PNPPS*PPSPTH

	PNPPH=PNPX*PXPPH+PNPY*PYPPH+PNPZ*PZPPH+PNPPS*PPSPPH

	PNPVR=PNPPS*(VR*YL2/V2-YLV*YR)

	PNPVTH=PNPPS*(VTH*YL2/V2-YLV*YTH)

	PNPVPH=PNPPS*(VPH*YL2/V2-YLV*YPH)

	NNP=N2-(2.0D0*X*PNPX+Y*PNPY+Z*PNPZ)

	PNPT=PNPX*PXPT

	SPACE=DREAL(N2).EQ.1.0D0.AND.DABS(DIMAG(N2)).LT.ABSLIM

	POLAR=-I*DSQRT(DABS(V2))*(-YT2+RAD)/(2.0D0*VDOTY*UX)

	GAM=(-YT2+RAD)/(2.0D0*UX)

	LPOLAR=I*X*DSQRT(DABS(YT2))/(UX*(U+GAM))

	KAY2=OM2/C2*N2

	IF (RSTART.EQ.0.) GO TO 1

	SCALE=DSQRT(DABS(DREAL(KAY2)/K2))

	KR=SCALE*KR

	KTH=SCALE*KTH

	KPH=SCALE*KPH

1	CONTINUE

C********* CALCULATES A HAMILTONIAN H

	H=.5D0*(C2*K2/OM2-N2)

C********* AND ITS PARTIAL DERIVATIVES WITH RESPECT TO

C********* TIME, R, THETA, PHI, OMEGA, KR, KTHETA AND KPHI.

	PHPT=-PNPT

	PHPR=-PNPR

	PHPTH=-PNPTH

	PHPPH=-PNPPH

	PHPOM=-NNP/OM

	PHPKR=C2/OM2*KR-C/OM*PNPVR

	PHPKTH=C2/OM2*KTH-C/OM*PNPVTH

	PHPKPH=C2/OM2*KPH-C/OM*PNPVPH

	KPHPK=N2

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE HAMLTN

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

	DOUBLE PRECISION PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN,W(400),

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

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

     3     KPHPK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,KAY2,KAY2I

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

	CHARACTER*8 MODRIN(3)

	COMMON /RINC/ MODRIN

C	 ****** SPACE IS LOGICAL ******

	COMMON /RIN/ COLL,FIELD,SPACE,KAY2,KAY2I,

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

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

     3     KPHPK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN

	LOGICAL SPACE

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R(20),T,STP,DRDT(20)

	COMMON R,T,STP,DRDT /WW/ ID,W0,W

	DOUBLE PRECISION TH,PH,KR,KTH,KPH,DTHDT,DPHDT,DKRDT,DKTHDT,

     1    DKPHDT,F,OM,STH,CTH,RSTH,RCTH

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

     1    (KPH,R(6)),(DTHDT,DRDT(2)),(DPHDT,DRDT(3)),(DKRDT,DRDT(4)),

     2    (DKTHDT,DRDT(5)),(DKPHDT,DRDT(6)),(F,W(6))

	OM=PIT2*1.D6*F

	STH=DSIN(TH)

	CTH=DSIN(PID2-TH)

	RSTH=R(1)*STH

	RCTH=R(1)*CTH

	CALL RINDEX

C	CHANGED ALL DRDT'S TO DRDT(1)	**************

	DRDT(1)=-PHPKR/(PHPOM*C)

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

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

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

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

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

	NR=6

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

	IF (W(57).EQ.0.0D0) GO TO 10

	NR=NR+1

	DRDT(NR)=-KPHPK/PHPOM/OM

C********* ABSORPTION

10	IF (W(58).EQ.0.0D0) GO TO 15

	NR=NR+1

	DRDT(NR)=10.0D0/LOGTEN*KPHPK*KAY2I/(KR*KR+KTH*KTH+KPH*KPH)/

     1     PHPOM/C

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

15	IF (W(59).EQ.0.0D0) GO TO 20

	NR=NR+1

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

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

20	IF (W(60).EQ.0.0D0) GO TO 25

	NR=NR+1

	DRDT(NR)=-DSQRT(DABS(PHPKR**2+PHPKTH**2+PHPKPH**2))/PHPOM/C

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

25	CONTINUE

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE RKAM

C	   NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS

	DOUBLE PRECISION SPACE,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT

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

	DOUBLE PRECISION Y(20),T,STEP,DYDT(20)

	COMMON Y,T,STEP,DYDT

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

	DOUBLE PRECISION YU,ALPHA,E,EPM,XV,FV,DELY,Z,DEL,SSE,EPSIL

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

	LL=1

	MM=1

	IF (MODE.EQ.1.0D0) MM=4

	ALPHA=T

	EPM=0.0D0

	BET(1)=0.5

	BET(2)=0.5

	BET(3)=1.0

	BET(4)=0.0

	STEP=SPACE

	R=19.0/270.0

	XV(MM)=T

	IF (E1MIN.LE.0.0D0) E1MIN=E1MAX/55.0D0

	IF (FACT.LE.0.0D0) FACT=0.5D0

	CALL HAMLTN

	DO 320 I=1,NN

	   FV(MM,I)=DYDT(I)

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

	RSTART=0.

	GO TO 1001

1000	IF (MODE.NE.1.0D0) GO TO 2000

C

C		RUNGE-KUTTA

1001	DO 1034 K=1,4

	   DO 1350 I=1,NN

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

	      Z=YU(MM,I)

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

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

	   CALL HAMLTN

	   DO 1034 I=1,NN

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

	   DO 1039 I=1,NN

	      DEL=(DELY(1,I)+2.0D0*DELY(2,I)+2.0D0*DELY(3,I)+

     1           DELY(4,I))/6.0D0

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

	   MM=MM+1

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

	   DO 1400 I=1,NN

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

	   T=XV(MM)

	   CALL HAMLTN

	   IF (MODE.EQ.1.0D0) GO TO 42

	   DO 150 I=1,NN

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

	   IF (MM.LE.3) GO TO 1001

C

C		ADAMS-MOULTON

2000	   DO 2048 I=1,NN

	      DEL=STEP*(55.0D0*FV(4,I)-59.0D0*FV(3,I)+37.0D0*FV(2,I)-

     1           9.0D0*FV(1,I))/24.0D0

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

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

	   T=XV(4)+STEP

	   CALL HAMLTN

	   XV(5)=T

	   DO 2051 I=1,NN

	      DEL=STEP*(9.0D0*DYDT(I)+19.0D0*FV(4,I)-5.0D0*FV(3,I)+

     1           FV(2,I))/24.0D0

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

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

	   CALL HAMLTN

	   IF (MODE.LE.2.0D0) GO TO 42

C

C		ERROR ANALYSIS

	   SSE=0.0D0

	   DO 3033 I=1,NN

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

	      IF (MODE.EQ.3.0D0.AND.Y(I).NE.0.0D0)

     1          EPSIL=EPSIL/DABS(Y(I))

	      IF (SSE.LT.EPSIL) SSE=EPSIL

3033	   CONTINUE

	   IF (E1MAX.GT.SSE) GO TO 3035

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

	   LL=1

	   MM=1

	   STEP=STEP*FACT

	   GO TO 1001

3035	   IF (LL.LE.1.OR.SSE.GE.E1MIN.OR.E2MAX.LE.DABS(STEP))

     1        GO TO 42

	   LL=2

	   MM=3

	   XV(2)=XV(3)

	   XV(3)=XV(5)

	   DO 5363 I=1,NN

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

	      FV(3,I)=DYDT(I)

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

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

	   STEP=2.0D0*STEP

	   GO TO 1001

C

C		EXIT ROUTINE

42	   LL=2

	   MM=4

	   DO 12 K=1,3

	      XV(K)=XV(K+1)

	      DO 12 I=1,NN

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

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

	   XV(4)=XV(5)

	   DO 52 I=1,NN

	      FV(4,I)=DYDT(I)

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

	   IF (MODE.LE.2.0D0) RETURN

	   E=DABS(XV(4)-ALPHA)

	   IF (E.LE.EPM) GO TO 2000

	   EPM=E

	   RETURN

	   END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

   	SUBROUTINE PRINTR(NWHY,CARD)

C		PRINTS OUTPUT AND PUNCHES RAYSETS WHEN REQUESTED

	DOUBLE PRECISION CARD

	CHARACTER*8 NWHY

	DIMENSION G(3,3),G1(3,3),RPRINT(20),NPR(20)

	CHARACTER*1 TYPE(3)

	CHARACTER*6 HEADR1(20),HEADR2(20),UNITS(20),HEAD1(20),HEAD2(20)

     1              ,UNIT(20)

	DOUBLE PRECISION PI,PIT2,PID2,DEGS,RAD,W(400),HPUNCH,N2,N2I,

     1         POLAR(2),LPOLAR(2),NTYP,DUM(3),G,G1,DENM

	COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,DUM

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

	CHARACTER*8 MODRIN(3)

	COMMON /RINC/ MODRIN

C	 **** SPACE IS LOGICAL ****

	COMMON /RIN/ COLL,FIELD,SPACE,N2,N2I,PNP(10),

     1     POLAR,LPOLAR,SGN

	COMPLEX*16 PNP

	LOGICAL SPACE,NEWWR,NEWWP,PENET

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R(20),T

	COMMON R,T /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,XMTRH,TLAT,TLON,F,AZ1,BETA,RCVRH,

     1     HOP,PLAT,PLON,RAYSET,THETA,PHI,SPL,CPL,SP,CP,SL,CL,R0,

     2     XR,YR,ZR,CTHR,STHR,PHIR,ALPH,PDEV,ABSORB,DOPP,TLOND,TLATD,

     3     AZ,EL,V,H,STH,CTH,XP,YP,ZP,EPS,ETA,ZETA,RCE2,RCE,RANGE

     4     ELL,SR,RPRINT,ANGA,AZDEV,ANA,SANA,SPHI,CPHI,AZA,TDEV

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

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

     2             (HOP,W(22)),(PLAT,W(24)),(PLON,W(25)),

     3             (RAYSET,W(72)),(THETA,R(2)),(PHI,R(3))

	DATA TYPE /1HX,1HN,1HO/,

     1  HEADR1(7) /6H  PHAS/,HEADR2(7) /6HE PATH/,UNITS(7) /6H   KM /,

     2  HEADR1(8) /6H  ABSO/,HEADR2(8) /6HRPTION/,UNITS(8) /6H   DB /,

     3  HEADR1(9) /6H   DOP/,HEADR2(9) /6HPLER  /,UNITS(9) /6H  C/S /,

     4  HEADR1(10) /6H PATH /,HEADR2(10) /6HLENGTH/,

     5  UNITS(10) /6H   KM /

	CALL RINDEX

	IF (.NOT.NEWWP) GO TO 10

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

	NEWWP=.FALSE.

	SPL=DSIN(PLON-TLON)

	CPL=DSIN(PID2-(PLON-TLON))

	SP=DSIN(PLAT)

	CP=DSIN(PID2-PLAT)

	SL=DSIN(TLAT)

	CL=DSIN(PID2-TLAT)

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

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

	G(1,2)=SPL*SP

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

	G(2,1)=-SPL*CL

	G(2,2)=CPL

	G(2,3)=SL*SPL

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

	G(3,2)=CP*SPL

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

	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)*

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

     2  G(3,2)*G(2,3)

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

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

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

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

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

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

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

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

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

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

	R0=EARTHR+XMTRH

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

	XR=R0*G(1,1)

	YR=R0*G(2,1)

	ZR=R0*G(3,1)

	CTHR=G(3,1)

	STHR=DSIN(DACOS(CTHR))

	PHIR=DATAN2(YR,XR)

	ALPH=DATAN2(G(3,2),G(3,3))

C*********

	NR=6

	NP=0

	DO 7 NN=7,20

	   IF (W(NN+50).EQ.0.0D0) GO TO 7

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

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

	   NR=NR+1

	   IF (W(NN+50).NE.2.0D0) GO TO 7

C********* DEPENDENT VARABLE NUMBER NN IS BEING INTEGRATED AND PRINTED.

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

C********* PRINTED

	   NP=NP+1

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

	   NPR(NP)=NR

	   HEAD1(NP)=HEADR1(NN)

	   HEAD2(NP)=HEADR2(NN)

	   UNIT(NP)=UNITS(NN)

7	CONTINUE

	NP1=MIN0(NP,3)

	PDEV=0.0D0

	ABSORB=0.0D0

	DOPP=0.0D0

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

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

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

	WRITE (4,1100) (HEAD1(NN),HEAD2(NN),NN=1,NP1)

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

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

     2   26HPOLARIZATION    GROUP PATH,5A6,A5)

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

	WRITE (4,1150) (UNIT(NN),NN=1,NP1)

1150	FORMAT (13X,2(8X,2HKM),2X,2(6X,3HDEG,5X,3HDEG),6X,

     1   12HREAL    IMAG,7X,2HKM,4X,3(4X,A6,2X))

	IF (RAYSET.EQ.0.0D0) GO TO 12

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

	TLOND=TLON*DEGS

	IF (TLOND.LT.0.0D0) TLOND=TLOND+360.0D0

	TLATD=TLAT*DEGS

	IF (TLATD.LT.0.0D0) TLATD=TLATD+360.0D0

	AZ=AZ1*DEGS

	EL=BETA*DEGS

	NHOP=IDINT(HOP)

	WRITE(2,1200) ID(1),TYPE(NTYP),XMTRH,TLATD,TLOND,RCVRH,F,AZ,EL,

     1    POLAR(1),POLAR(2),NHOP

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

     1     'TITLE',0P)

C********* 

12	V=0.0D0

	IF (N2.NE.0.0D0) V=(R(4)**2+R(5)**2+R(6)**2)/N2-1.0D0

	H=R(1)-EARTHR

	STH=DSIN(THETA)

	CTH=DSIN(PID2-THETA)

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

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

	YP=R(1)*STH*DSIN(PHI)-YR

	ZP=R(1)*CTH-ZR

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

C********* AND ROTATED

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

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

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

	RCE2=ETA**2+ZETA**2

	RCE=DSQRT(DABS(RCE2))

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

	RANGE=EARTHR*DATAN2(RCE,EARTHR+EPS+XMTRH)

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

	ELL=DATAN2(R(4),DSQRT(DABS(R(5)**2+R(6)**2)))*DEGS

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

	SR=DSQRT(DABS(RCE2+EPS**2))

	IF (NP.LE.1) GO TO 16

	DO 15 I=1,NP

	   NN=NPR(I)

15	RPRINT(I)=R(NN)

16	IF (SR.GE.1.0D-6) GO TO 20

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

C********* TRANSMITTER

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

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

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

	GO TO 40

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

20	EL=DATAN2(EPS,RCE)*DEGS

	IF (RCE.GE.1.0D-6) GO TO 30

C********* NEARLY DIRECTLY ABOVE OR BELOW TRANSMITTER. CANNOT CALCULTATE

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

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

	WRITE (4,2500) V,NWHY,H,RANGE,EL,ELL,POLAR,T,(RPRINT(NN),

     1     NN=1,NP1)

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

     1   4F12.4)

	GO TO 40

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

30	ANGA=DATAN2(ETA,ZETA)

	AZDEV=180.0D0-DMOD(540.0D0-(AZ1-ANGA)*DEGS,360.0D0)

	IF (R(5).NE.0.0D0.OR.R(6).NE.0.0D0) GO TO 34

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

C********* CALCULATED

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

C     1   NN=1,NP1)

	WRITE (4,3000) V,NWHY,H,RANGE,AZDEV,EL,ELL,POLAR,T,

     1     (RPRINT(NN),NN=1,NP1)

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

     1   F8.3,4F12.4)

	GO TO 40

34	ANA=ANGA-ALPH

	SANA=DSIN(ANA)

	SPHI=SANA*STHR/STH

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

     1   *CTHR

	AZA=180.0D0-DMOD(540.0D0-(DATAN2(SPHI,CPHI)-DATAN2(R(6),R(5)))*

     1     DEGS,360.0D0)

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

C     1   ,NN=1,NP1)

	WRITE (4,3500) V,NWHY,H,RANGE,AZDEV,AZA,EL,ELL,POLAR,T,

     1    (RPRINT(NN),NN=1,NP1)

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

     1   4F12.4)

C********* 

40	LINES=LINES+1

	IF (NP.LE.3) GO TO 45

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

C********* VARIABLES

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

	WRITE (4,4000) (RPRINT(NN),NN=4,NP)

4000	FORMAT (99X,3F12.4)

	LINES=LINES+1

45	IF (CARD.EQ.0.0D0) RETURN

C

C********* PUNCH A RAY SET

	IF (AZDEV.LT.-90.0D0) AZDEV=AZDEV+360.0D0

	IF (AZA.LT.-90.0D0) AZA=AZA+360.0D0

	TDEV=T-SR

	NR=6

	IF (W(57).EQ.0.0D0) GO TO 47

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

	NR=NR+1

	PDEV=R(NR)-SR

47	IF (W(58).EQ.0.0D0) GO TO 48

C********* ABSORPTION

	NR=NR+1

	ABSORB=R(NR)

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

48	IF (W(59).NE.0.0D0) DOPP=R(NR+1)

	WRITE(2,4500) HPUNCH,RANGE,AZDEV,AZA,ELL,SR,TDEV,PDEV,ABSORB

     1   ,DOPP,POLAR(1),POLAR(2),IHOP,NWHY

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

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE POL CAR

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

	COMMON R(6) /COORD/ S

	DOUBLE PRECISION PI,PIT2,PID2,W(400),R,DUM(5),SINA,COSA,SINP,

     1     COSP,X0,X,R0,TEMP,S

	COMMON /CONST/ PI,PIT2,PID2,DUM

C

C           CONVERTS SPHERICAL COORDINATES TO CARTESIAN

	IF (R(5).EQ.0.0D0.AND.R(6).EQ.0.0D0) GO TO 1

	VERT=0.

	SINA=DSIN(R(2))

	COSA=DSIN(PID2-R(2))

	SINP=DSIN(R(3))

	COSP=DSIN(PID2-R(3))

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

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

	X0(3)=R(1)*COSA

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

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

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

	RETURN

C		VERTICAL INCIDENCE

1	VERT=1.

	R0(1)=R(1)

	R0(2)=R(2)

	R0(3)=R(3)

	R0(4)=DSIGN(1.0D0,R(4))

	RETURN

C

C		STEPS THE RAY A DISTANCE S, AND THEN

C		CONVERTS CARTESIAN COORDINATES TO SPHERICAL COORDINATES

	ENTRY CAR POL

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

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

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

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

	TEMP=DSQRT(DABS(X(1)**2+X(2)**2))

	R(1)=DSQRT(DABS(X(1)**2+X(2)**2+X(3)**2))

	R(2)=DATAN2(TEMP,X(3))

	R(3)=DATAN2(X(2),X(1))

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

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

     1   *TEMP)

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

	RETURN

C                    VERTICAL INCIDENCE

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

	R(2)=R0(2)

	R(3)=R0(3)

	R(4)=R0(4)

	R(5)=0.

	R(6)=0.

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE BACK UP(HS)

	DOUBLE PRECISION HS,W(400)

	DOUBLE PRECISION STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT

	CHARACTER*8 STR8

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

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

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R(20),T,STP,DRDT(20)

	COMMON R,T,STP,DRDT /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,INTYP,STEP1,SMT,H

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

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

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

C 	STR8='BACK UP0'

C	CALL PRINTR (STR8,0.0D0)

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

	HOME=DRDT(1)*(R(1)-EARTHR-HS).GE.0.0D0

	IF (HS.GT.0.0D0.AND..NOT.HOME.OR.HS.EQ.0.0D0.AND.

     1     DRDT(1).GT.0.0D0) GO TO 30

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

	DO 10 I=1,10

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

	 STEP=DSIGN(DMIN1(DABS(STP),DABS(STEP)),STEP)

	 IF (DABS(R(1)-EARTHR-HS).LT..5D-4.AND.DABS(STEP).LT.1.0D0)

     1       GO TO 60

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

C 	 STR8='BACK UP1'

C	 CALL PRINTR(STR8,0.0D0)

	 MODE=1.0D0

	 RSTART=1.

	 CALL RKAM

10	RSTART=1.

C

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

	ENTRY GRAZE (HS)

C************

	THERE=.FALSE.

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

C 	STR8='GRAZE 0 '

C	CALL PRINTR (STR8,0.0D0)

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

	DO 20 I=1,10

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

	  STEP=DSIGN(DMIN1(DABS(STP),DABS(STEP)),STEP)

	  IF (DABS(R(4)).LE.1.D-6.AND.DABS(STEP).LT.1.0D0) GO TO 60

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

C 	  STR8='GRAZE 1 '

C         CALL PRINTR (STR8,0.0D0)

	  MODE=1.0D0

	  RSTART=1.

	  CALL RKAM

	  RSTART=1.

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

	  IF (R(5).EQ.0.0D0.AND.R(6).EQ.0.0D0) GO TO 60

20	CONTINUE

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

C************ PROBABLY MEANS THAT THE RAY INTERSECTS THE HEIGHT HS

30	CONTINUE

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

C 	STR8='BACK UP2'

C	CALL PRINTR (STR8,0.0D0)

	MODE=1.0D0

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

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

	H=R(1)-EARTHR

	STEP=(-R(4)-DSQRT(DABS(R(4)**2-2.*(H-HS)*R(4)/DRDT(1)*

     1     DRDT(4))))/DRDT(4)

	RSTART=1.

	CALL RKAM

	RSTART=1.

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

	DO 40 I=1,10

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

	   STEP=DSIGN(DMIN1(DABS(STP),DABS(STEP)),STEP)

	   IF (DABS(R(1)-EARTHR-HS).LT..5D-4.AND.DABS(STEP).LT.1.0D0)

     1        GO TO 50

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

C 	   STR8='BACK UP3'

C	   CALL PRINTR(STR8,0.0D0)

	   MODE=1.0D0

	   RSTART=1.

	   CALL RKAM

40	RSTART=1.

50	THERE=.TRUE.

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

60	MODE=INTYP

	STEP=STEP1

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE REACH

C           CALCULATES THE STRAIGHT LINE RAY PATH BETWEEN THE EARTH

C           AND THE IONOSPHERE OR BETWEEN IONOSPHERIC LAYERS

	CHARACTER*8 STR8

	DOUBLE PRECISION W(400),N2,N2I,SMT

	DOUBLE PRECISION STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT

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

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

	COMMON /COORD/ S

	CHARACTER*8 MODRIN(3)

	COMMON /RINC/ MODRIN

C	 **** SPACE IS LOGICAL ****

	COMMON /RIN/ COLL,FIELD,SPACE,N2,N2I,PNP(10),

     1               POLAR,LPOLAR,SGN

	COMPLEX*16 PNP,POLAR,LPOLAR

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R(20),T,STP,DRDT(20)

	COMMON R,T,STP,DRDT /WW/ ID,W0,W

	DOUBLE PRECISION EARTHR,XMTRH,RCVRH,S,SSTEP,H,V,UP,RADG,

     1                   DISTG,SG,RADR,DISTR,SR,S1

	EQUIVALENCE (EARTHR,W(2)),(XMTRH,W(3)),(RCVRH,W(20))

	LOGICAL CROSS,CROSSG,CROSSR,SPACE,GROUND,PERIGE,THERE,MINDIS,

     1        NEWRAY,RSPACE

	DATA NSTEP /500/

	CALL HAMLTN

	H=R(1)-EARTHR

 	STR8='EXIT ION'

	IF (.NOT.NEWRAY.AND..NOT.RSPACE) CALL PRINTR(STR8,0.0D0)

	NEWRAY=.FALSE.

	V=DSQRT(DABS(R(4)**2+R(5)**2+R(6)**2))

C********** NORMILIZE THE WAVE NORMAL DIRECTION TO ONE

	R(4)=R(4)/V

	R(5)=R(5)/V

	R(6)=R(6)/V

C********** NEGATIVE OF DISTANCE ALONG RAY TO CLOSEST APPROACH TO CENTER

C********** OF EARTH

	UP=R(1)*R(4)

	RADG=EARTHR**2-R(1)**2*(R(5)**2+R(6)**2)

	DISTG=DSQRT(DABS(DMAX1(0.0D0,RADG)))

C********** DISTANCE ALONG RAY TO FIRST INTERSECTION WITH OR CLOSEST

C********** APPROACH TO THE EARTH

	SG=-UP-DISTG

C********** CROSSG IS TRUE IF THE RAY WILL INTERSECT OR TOUCH THE EARTH

	CROSSG=UP.LT.0.0D0.AND.RADG.GE.0.0D0

	RADR=(EARTHR+RCVRH)**2-R(1)**2*(R(5)**2+R(6)**2)

	DISTR=DSQRT(DABS(DMAX1(0.0D0,RADR)))

C********** DISTANCE ALONG RAY TO THE FIRST INTERSECTION WITH OR CLOSEST

C********** APPROACH TO THE RECEIVER HEIGHT

	SR=DISTR-UP

	IF (UP.LT.0.0D0.AND.DISTR.LT.-UP.AND.R(1).NE.EARTHR+RCVRH)

     1       SR=-DISTR-UP

C********** CROSSR IS TRUE IF THE RAY WILL INTERSECT WITH OR MAKE A

C********** CLOSEST APPROACH TO THE RECEIVER HEIGHT

	 CROSSR=R(4).LT.0.0D0.OR.R(1).LT.(EARTHR+RCVRH)

	 CROSS=CROSSG.OR.CROSSR

C********** MAXIMUM DISTANCE IN WHICH TO LOOK FOR THE IONOSPHERE

	 S1=DMIN1(SR,SG)

	 IF (.NOT.CROSSG) S1=SR

	 IF (UP.GE.0.0D0) GO TO 15

	 CROSS=.TRUE.

C********** IF RAY IS GOING DOWN, S1 IS AT MOST THE DISTANCE TO A PERIGEE

	 S1=DMIN1(S1,-UP)

C********** CONVERT THE POSITION AND DIRECTION OF THE RAY TO CARTESIAN

C********** COORDINATES

15	CALL POL CAR

	SSTEP=100.0D0

	S=SSTEP

	DO 20 I=1,NSTEP

	 IF ((S-SSTEP).GT.S1.AND.CROSS) GO TO 25

C********** CONVERT POSITION AND DIRECTION TO SPHERICAL POLAR COORDINATES

C********** AT A DISTANCE S ALONG THE RAY

	 CALL CAR POL

	 CALL ELECTX

C********** FREE SPACE

	 IF (X.EQ.0.0D0) GO TO 20

	 CALL RINDEX

C********** EFFECTIVELY FREE SPACE

	 IF (SPACE) GO TO 20

	 IF (SSTEP.LT.0.5D-4) GO TO 25

C********** RAY IN THE IONOSPHERE. STEP BACK OUT

	 S=S-SSTEP

C********** DECREASE STEP SIZE

	 SSTEP=SSTEP/10.0D0

20	 S=S+SSTEP

	 PRINT 2000, NSTEP

2000	 FORMAT (9H EXCEEDED,I5,26H STEPS IN SUBROUTINE REACH)

	 CALL EXIT

25	IF (CROSS) S=DMIN1(S,S1)

C********** CONVERT POSITION AND DIRECTION TO SPHERICAL POLAR COORDINATES

C********** AT A DISTANCE S ALONG THE RAY

	CALL CAR POL

C********** AVOID THE RAY BEING SLIGHTLY UNDERGROUND

	R(1)=DMAX1(R(1),EARTHR)

C********** ONE STEP INTEGRATION

	IF (N.LT.7) GO TO 31

	DO 30 NN=7,N

30	R(NN)=R(NN)+S*DRDT(NN)

31	T=T+S

	CALL RINDEX

C********** AT A PERIGEE

	PERIGE=S.EQ.(-UP)

C********** CORRECT MINOR ERRORS

	IF (PERIGE) R(4)=0.0D0

C********** KEEP CONSISTENCY AFTER CORRECTING MINOR ERRORS

	DRDT(1)=R(4)

C********** ON THE GROUND

	GROUND=S.EQ.SG.AND.CROSSG

C********** AT THE RECEIVER HEIGHT

	THERE=S.EQ.SR.AND.CROSSR.AND..NOT.PERIGE

C********** AT A CLOSEST APPROACH TO THE RECEIVER HEIGHT

	MINDIS=PERIGE.AND.S.EQ.SR.AND.CROSSR

	RSPACE=SPACE

	V=DSQRT(DABS(N2/(R(4)**2+R(5)**2+R(6)**2)))

C********** RENORMALIZE THE WAVE NORMAL DIRECTION TO = SQRT(REAL(N**2))

	R(4)=R(4)*V

	R(5)=R(5)*V

	R(6)=R(6)*V

	RSTART=1.

 	STR8='ENTR ION'

	IF (.NOT.SPACE) CALL PRINTR (STR8,0.0D0)

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE LABPLT

C	LABEL THE CURRENT PLOT

	CHARACTER*8 NDATE,ID(10),TYPE(3)

	COMMON /LPLT/ NDATE

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

	COMMON /CONST/ PI,PIT2,PID2,DEGS,DUM

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

	COMMON /WW/ ID,W0,W

	DOUBLE PRECISION W(400),EARTHR,F,AZ1,PLT,FACTR,TIC,PI,PIT2,

     1     PID2,DEGS,DUM(4),LTIC,DTIC,HPUNCH,NTYP

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

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

	LOGICAL NEWWR,NEWWP,PENET

	DATA TYPE /8HEXTRAORD,8HNO FIELD,8HORDINARY/

	WRITE (3,*) 'PU;PA0,1023;LB',(ID(N),N=1,8),CHAR(3)

	WRITE (3,*) 'PU;PA900,1023;LB',NDATE,CHAR(3)

	AZA=AZ1*DEGS

	DTIC=TIC*EARTHR

	WRITE (3,1000) F,AZA,TYPE(NTYP),DTIC,CHAR(3)

1000	FORMAT(16HPU;PA0,991;LBF =,F7.3,6H, AZ =,F7.2,2H, ,A8,2H, ,

     1     F7.2,22H KM BETWEEN TICK MARKS,A1)

	IF(PLT.EQ.1) RETURN

	LTIC=DTIC/FACTR

	WRITE (3,2000) LTIC,CHAR(3)

2000	FORMAT(17HPU;PA0,0;DI0,1;LB,F7.2,22H KM BETWEEN TICK MARKS,

     1     A1)

	RETURN

	END	

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE PLOT (X,Y,NEW)

	DOUBLE PRECISION X,Y

	DOUBLE PRECISION XMIN0,XMAX0,YMIN0,YMAX0

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

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

	DATA MINX,MINY /2*0/,MAXX,MAXY /2*1023/,

     1	  MINX0,MINY0 /2*23/,MAXX0,MAXY0 /2*1023/

C

C		COMPUTE SCALE FACTORS

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

	RESET=0.

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

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

	XMIN=XMIN0-(MINX0-MINX)/XSCALE

	YMIN=YMIN0-(MINY0-MINY)/YSCALE

	XMAX=XMAX0-(MAXX-MAXX0)/XSCALE

	YMAX=YMAX0-(MAXY-MAXY0)/YSCALE

C

C		START A NEW LINE

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

	IX=MINX0+(X-XMIN0)*XSCALE

	IY=MINY0+(Y-YMIN0)*YSCALE

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

     1     WRITE (3,*) 'PU;PA',IX,IY,'PD;'

	GO TO 50

C

C		HORIZONTAL DISPLACEMENT

10	XS=X-XOLD

	YS=Y-YOLD

	IF (XS) 11,12,16

C		NEGATIVE

11	X1=XMAX

	X2=XMIN

	GO TO 20

C		ZERO

12	IF (YS) 13,50,14

13	S1=(YMAX-YOLD)/YS

	S2=(YMIN-YOLD)/YS

	GO TO 40

14	S1=(YMIN-YOLD)/YS

	S2=(YMAX-YOLD)/YS

	GO TO 40

C		POSITIVE

16	X1=XMIN

	X2=XMAX

C

C		VERTICAL DISPLACEMENT

20	IF (YS) 21,22,26

C		NEGATIVE

21	Y1=YMAX

	Y2=YMIN

	GO TO 30

C		ZERO

22	S1=(X1-XOLD)/XS

	S2=(X2-XOLD)/XS

	GO TO 40

C		POSITIVE

26	Y1=YMIN

	Y2=YMAX

C

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

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

C

C		PLOT LINE -- CHECKING FOR BORDER CROSSINGS

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

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

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

C		PREVIOUS POINT OFF GRAPH

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

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

	WRITE (3,*) 'PU;PA',IX,IY,'PD;'

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

C		CURRENT POINT OFF GRAPH

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

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

	WRITE (3,*) 'PA',IX,IY,';'

C		CURRENT POINT ON GRAPH

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

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

	WRITE (3,*) 'PA',IX,IY,';'

C		EXIT ROUTINE

50	XOLD=X

	YOLD=Y

	RETURN

C

C		TERMINATE THE CURRENT PLOT

	ENTRY PLTEND

	WRITE (*,*) 'Wait Until Plotting is Complete.'

	CLOSE (3)

	OPEN(UNIT=3,FILE='COM2',STATUS='OLD')

	WRITE(3,*) 'IN;IP;SC0,1023,0,1023;SP1;'

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE RAYPLT

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

C		=2. PLOTS PROJECTION OF RAYPATH ON GROUND

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

	DOUBLE PRECISION PI,PIT2,PID2,W(400),HPUNCH,NTYP,DUM(5)

	COMMON /CONST/ PI,PIT2,PID2,DUM

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

	CHARACTER*8 ID(10)

	COMMON R /WW/ ID,W0,W

	DOUBLE PRECISION TH,PH,EARTHR,PLAT,PLON,PLT,FACTR,LLAT,LLON,

     1    RLAT,RLON,TIC,HB,LTIC,SW,CW,SLM,CLM,SRM,CRM,CDPHI,PHL,CTHL,

     2    STHL,PHR,CTHR,STHR,CLR,SLR,R0,ALPHA,XR,XL,YB,YT,ALPH1,STH,

     3    CTH,CR,CL,CEA,SL,TMP1,TMP2,ALPH2,TICKX,R1,X,Y,ANG,DTIC,

     4    TICY,TIC1,R(6)

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

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

     1    W(81)),(FACTR,W(82)),(LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),

     2    (RLON,W(86)),(TIC,W(87)),(HB,W(88))

	LOGICAL NEWWR,NEWWP,PENET

	IF (.NOT. NEWWR) GO TO 5

C

C	     NEW W ARRAY -- REINITIALIZE

	NEWWR=.FALSE.

	RESET=1.

CONVERT COORDINATES OF VERTICAL PLANE FROM GEOGRAPHIC TO GEOMAGNETIC

	SW=DSIN(PLAT)

	CW=DSIN(PID2-PLAT)

	SLM=DSIN(LLAT)

	CLM=DSIN(PID2-LLAT)

	SRM=DSIN(RLAT)

	CRM=DSIN(PID2-RLAT)

	CDPHI=DSIN(PID2-(LLON-PLON))

	PHL=DATAN2(DSIN(LLON-PLON)*CLM,CDPHI*SW*CLM-CW*SLM)

	CTHL=CDPHI*CW*CLM+SW*SLM

	STHL=DSIN(DACOS(CTHL))

	CDPHI=DSIN(PID2-(RLON-PLON))

	PHR=DATAN2(DSIN(RLON-PLON)*CRM,CDPHI*SW*CRM-CW*SRM)

	CTHR=CDPHI*CW*CRM+SW*SRM

	STHR=DSIN(DACOS(CTHR))

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

	SLR=DSQRT(DABS(1.0D0-CLR**2))

	IF (PLT.EQ.2.0D0) GO TO 3

	FACTR=1.0D0

	R0=EARTHR+HB

	ALPHA=.5D0*DACOS(CLR)

	XR=R0*DSIN(ALPHA)

	XL=-XR

	YB=R0*DSIN(PID2-ALPHA)

	YT=YB+2.0D0*XR

	GO TO 5

3	IF (FACTR.EQ.0.0D0) FACTR=1.0D0

	ALPH1=DATAN2(STHR*DSIN(PHR-PHL),(CTHR-CTHL*CLR)/STHL)

	XL=0.0D0

	XR=EARTHR*DACOS(CLR)

	YT=0.5D0*XR/FACTR

	YB=-YT

5	STH=DSIN(TH)

	CTH=DSIN(PID2-TH)

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

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

	CEA=DATAN2(CR-CL*CLR,CL*SLR)

	NEW=1

	IF (IHOP.NE.0) NEW=0

	IF (PLT.EQ.2.0D0) GO TO 10

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

	RETURN

10	SL=DSQRT(DABS(1.0D0-CL**2))

	TMP1=STH*DSIN(PH-PHL)

	TMP2=(CTH-CTHL*CL)/STHL

	ALPH2=0.0D0

	IF (TMP1.NE.0.0D0.OR.TMP2.NE.0.0D0) ALPH2=DATAN2(TMP1,TMP2)

	CALL PLOT (EARTHR*CEA,EARTHR*DASIN(SL*DSIN(ALPH1-ALPH2)),NEW)

	RETURN

C

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

	ENTRY ENDPLT

	TICKX=0.01D0*(YT-YB)

	IF (PLT.EQ.2.0D0) GO TO 25

	R1=EARTHR-TICKX

	X=XL

	Y=YB

	CALL PLOT (X,Y,1)

	NTIC=2

	IF (TIC.NE.0.0D0) NTIC=NTIC+IDINT(2.0D0*ALPHA/TIC)

	NLINE=MAX0(1,100/NTIC)

	DO 20 I=1,NTIC

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

	   CALL PLOT (R1*DSIN(ANG),R1*DSIN(PID2-ANG),0)

	   CALL PLOT (X,Y,0)

	   DO 20 J=1,NLINE

	      ANG=ANG+TIC/NLINE

	      X=EARTHR*DSIN(ANG)

	      Y=EARTHR*DSIN(PID2-ANG)

20	CALL PLOT (X,Y,0)

	CALL PLOT (XR,YB,0)

	GO TO 50

25	DTIC=TIC*EARTHR

	LTIC=DTIC/FACTR

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

	NTIC=IDINT(YT/LTIC)

	TIC1=-LTIC*NTIC

	CALL PLOT (XL,YB,1)

	NTIC=2*NTIC+1

	DO 30 I=1,NTIC

	   Y=TIC1+(I-1)*LTIC

	   CALL PLOT (XL,Y,0)

	   CALL PLOT (TICY,Y,0)

30	   CALL PLOT (XL,Y,0)

	CALL PLOT (XL,YT,0)

	CALL PLOT (XL,0.0D0,1)

	NTIC=IDINT((XR-XL)/DTIC)

	DO 40 I=1,NTIC

	   X=I*DTIC

	   CALL PLOT (X,0.0D0,0)

	   CALL PLOT (X,TICKX,0)

40	   CALL PLOT (X,0.0D0,0)

	CALL PLOT (XR,0.0D0,0)

50	CALL LABPLT

	CALL PLTEND

	RETURN

	END

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

	SUBROUTINE TRACE                                                  

C         CALCULATES THE RAY PATH					 

	CHARACTER*8 STR8

	DOUBLE PRECISION W(400),HPUNCH

	DIMENSION ROLD(20),DROLD(20)					 

	DOUBLE PRECISION STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,NTYP

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

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

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

	CHARACTER*8 MODRIN(3)

	COMMON /RINC/ MODRIN

C	 **** SPACE IS LOGICAL ****

	COMMON /RIN/ COLL,FIELD,SPACE,N2,PNP(10),

     1     POLAR,LPOLAR,SGN

	COMPLEX*16 N2,PNP,POLAR,LPOLAR

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

     1    NEWWR,NEWWP,PENET,NEWRAY,WAS					 

	CHARACTER*6 MODX(2)

	COMMON /XXC/ MODX

	DOUBLE PRECISION X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	COMMON /XX/ X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  

	CHARACTER*8 ID(10)

	DOUBLE PRECISION R(20),T,STP,DRDT(20)

	COMMON R,T,STP,DRDT /WW/ID,W0,W			 

	DOUBLE PRECISION MAXSTP,EARTHR,RCVRH,HOP,SKIP,RAYSET,PLT,

     1     HTMAX,TOLD,APHT,H,ROLD,DROLD,TMP,SMT

	EQUIVALENCE (EARTHR,W(2)),(RCVRH,W(20)),(HOP,W(22)),

     1   (MAXSTP,W(23)),(SKIP,W(71)),(RAYSET,W(72)),(PLT,W(81))                          

	NHOP=IDINT(HOP)								 

	MAX=IDINT(MAXSTP)							 

	NSKIP=IDINT(SKIP)							 

	RSTART=1.						  	 

	CALL HAMLTN

	HOME=DRDT(1)*(R(1)-EARTHR-RCVRH).GE.0.0D0				 

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

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

	IHOP=0								 

 	STR8='XMTR    '

	CALL PRINTR (STR8,0.0D0)					 

	IF (PLT.NE.0.0D0) CALL RAYPLT					 

	HTMAX=0.0D0							 

	NEWRAY=.TRUE.							 

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

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

10	IHOP=IHOP+1							 

	IF (IHOP.GT.NHOP) RETURN						 

	PENET=.FALSE.							 

	APHT=RCVRH							 

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

	DO 79 J=1,MAX							 

	  H=R(1)-EARTHR							 

	  IF (DABS(H-RCVRH).GT.DABS(APHT-RCVRH)) APHT=H			 

	  HTMAX=DMAX1(H,HTMAX)						 

	  IF (.NOT.SPACE) GO TO 12						 

	  CALL REACH							 

	  RSTART=1.							 

	  H=R(1)-EARTHR							 

	  IF (DABS(H-RCVRH).GT.DABS(APHT-RCVRH)) APHT=H			 

	  HTMAX=DMAX1(H,HTMAX)						 

	  IF (.NOT.SPACE) GO TO 12

	  STR8='PERIGEE '

	  IF (PERIGE) CALL PRINTR (STR8,0.0D0)				 

	  IF (THERE) GO TO 51						 

	  IF (MINDIS) GO TO 40						 

	  IF (GROUND) GO TO 60						 

	  IF (PLT.NE.0.0D0) CALL RAYPLT					 

	  IF (PERIGE) GO TO 79						 

12	  DO 13 L=1,N							 

	   ROLD(L)=R(L)							 

13	  DROLD(L)=DRDT(L)							 

	  TOLD=T								 

	  WAS=THERE							 

	  CALL RKAM							 

	  H=R(1)-EARTHR							 

	  THERE=.FALSE.						 	 

	  WASNT=.NOT.HOME							 

          HOME=DRDT(1)*(H-RCVRH).GE.0.0D0					 

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

	  SMT=0.0D0

	  IF (TMP.NE.0.0D0) SMT=0.5D0*(R(1)-ROLD(1)+0.5D0*TMP)**2/

     1       DABS(TMP)

	  IF (((H-RCVRH)*(ROLD(1)-EARTHR-RCVRH).LT.0.0D0.AND..NOT.WAS)

     1       .OR.(WAS.AND.DRDT(1)*DROLD(1).LT.0.0D0.AND.HOME)) GO TO 50

	  IF (HOME.AND.WASNT) GO TO 30

	  IF (H.LT.0.0D0.OR.DRDT(1).GT.0.0D0.AND.DROLD(1).LT.0.0D0.AND.

     1      SMT.GT.H) GO TO 20

 	  STR8='PERIGEE '

	  IF (DROLD(1).LT.0.0D0.AND.DRDT(1).GT.0.0D0)

     1       CALL PRINTR(STR8,0.0D0)

 	  STR8='APOGEE  '

	  IF (DROLD(1).GT.0.0D0.AND.DRDT(1).LT.0.0D0)

     1       CALL PRINTR(STR8,0.0D0)

 	  STR8='MAX LAT '

	  IF (DROLD(2)*DRDT(2).LT.0.0D0) CALL PRINTR(STR8,0.0D0)

	  STR8='MAX LONG'

	  IF (DROLD(3)*DRDT(3).LT.0.0D0) CALL PRINTR(STR8,0.0D0)

	  DO 14 I=4,6

	   STR8='WAVE REV'

	   IF (ROLD(I)*R(I).LT.0.0D0) CALL PRINTR(STR8,0.0D0)

14	  CONTINUE

	  GO TO 75

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

20	  CALL BACK UP(0.0D0)

	  GO TO 60

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

30	  CALL GRAZE (RCVRH)

	  IF (THERE) GO TO 51

40	  DRDT(1)=0.0D0

	  HPUNCH=R(1)-EARTHR

	  STR8='MIN DIST'

	  CALL PRINTR(STR8,RAYSET)

	  IF (PLT.NE.0.0D0) CALL RAYPLT

	  IF (IHOP.GE.NHOP) RETURN

	  IHOP=IHOP+1

	  STR8='MIN DIST'

	  CALL PRINTR(STR8,RAYSET)

	  GO TO 89

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

50	  CALL BACK UP(RCVRH)

	  THERE=.TRUE.

51	  R(1)=EARTHR+RCVRH

	  HTMAX=DMAX1(RCVRH,HTMAX)

	  HPUNCH=APHT

	  STR8='RCVR    '

	  CALL PRINTR(STR8,RAYSET)

	  IF (PLT.NE.0.0D0) CALL RAYPLT

	  IF (RCVRH.NE.0.0D0) GO TO 89

	  IF (IHOP.GE.NHOP) RETURN

	  IHOP=IHOP+1

	  APHT=RCVRH

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

60	  R(1)=EARTHR

	  IF (DABS(RCVRH).GT.DABS(APHT-RCVRH)) APHT=0.

	  R(4)=DABS(R(4))

	  DRDT(1)=DABS(DRDT(1))

	  RSTART=1.

	  HPUNCH=HTMAX

	  STR8='GRND REF'

	  CALL PRINTR(STR8,RAYSET)

	  HTMAX=0.0D0

	  IF (RCVRH.NE.0.0D0) GO TO 65

	  THERE=.TRUE.

	  HPUNCH=APHT

	  STR8='RCVR    '

	  CALL PRINTR (STR8,RAYSET)

	  GO TO 89

65	  H=0.0D0

	  THERE=.FALSE.

C**********

75	  IF (PLT.NE.0.0D0) CALL RAYPLT

	  IF (H.GT.HMAX.AND.H.GT.RCVRH.AND.DRDT(1).GT.0.0D0) GO TO 90

	  STR8='        '

	  IF (MOD(J,NSKIP).EQ.0) CALL PRINTR(STR8,0.0D0)

79	CONTINUE

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

	HPUNCH=H

	STR8='STEP MAX'

	CALL PRINTR(STR8,RAYSET)

	RETURN

C***********

89	HOME=.TRUE.

	GO TO 10

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

90	PENET=.TRUE.

	HPUNCH=H

	STR8='PENETRAT'

	CALL PRINTR(STR8,RAYSET)

	RETURN

	END

