c ............. file name RAYTRACE_DP.FOR
PROGRAM NITIAL NITI001
implicit real*8 (a-h,o-z)
c
c Jones - Stephenson raytracing program
c----- last modified 23-mar-88
c
C SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE NITI002
real*4 time0,time1,time2,timef,secnds
real*4 sthl,sphl,splast,stlast,slat,slon,sdist,sbear
character*8 MFLD(2) NITI003
character*8 id,ndate,modx,mody,modz,modrin,koll,modsav
COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,C,LOGTEN NITI004
common /trac / ground,perige,there,mindis,newray,smt
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
c apogee comes from TRACE; range comes from PRINT
common / leo / apogee, range
dimension ww(100),nww(100)
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
logical ground,perige,there,mindis
REAL*8 N2,N2I,LOGTEN,K,MAXSTP,INTYP,MAXERR,MU NITI021
COMPLEX*16 PNP,POLAR,LPOLAR NITI022
c
c PRINT writes to unit 2
open (unit=2,file='data.out',status='old')
c the final ground range (ONE-HOP ONLY) is written to JONES.OUT
open (unit=4,file='jones.out',status='old')
c PUNCH writes to unit 7
open (unit=7,file='punch.out',status='old')
c normal input is from unit 1 -- read by s/r readW
open (unit=1,file='ray.dat',status='old')
c
ccc NDATE=IDATE(0) NITI025
ccc SECOND=KLOCK(0)*.001 NITI026
KOLL=' NO' NITI027
IF (COLL.NE.0.) KOLL='WITH' NITI028
C********* CONSTANTS NITI029
c PI=3.1415926536 NITI030
pi = datan(1.d0) * 4.d0
PIT2=2.d0*PI NITI031
PID2=PI/2.d0 NITI032
DEGS=180.d0/PI NITI033
RAD=PI/180.d0 NITI034
C=2.997925d5 NITI035
K=2.81785d-15*C**2/PI NITI036
LOGTEN=dLOG(10.d0) NITI037
C********* INITIALIZE SOME VARIABLES IN THE W ARRAY NITI038
DO 5 NW=1,400 NITI039
5 W(NW)=0.d0 NITI040
PLON=0.d0 NITI041
PLAT=PID2 NITI042
EARTHR=6370.d0 NITI043
INTYP=3. NITI044
MAXERR=1.d-4 NITI045
ERATIO=50.d0 NITI046
STEP1=1.d0 NITI047
STPMAX=100.d0 NITI048
STPMIN=1.d-8 NITI049
FACTR=0.5d0 NITI050
MAXSTP=1000.d0 NITI051
HOP=1.d0 NITI052
C********* READ W ARRAY AND PRINT NON-ZERO VALUES NITI053
10 CALL READ W NITI054
c for the no-field case, put he "geomc" pole at the geogc pole
if(ray.eq.0.) w(24) = pid2
if(ray.eq.0.) w(25) = 0.d0
IF (SKIP.EQ.0.) SKIP=MAXSTP NITI056
c 12 RAY=dSIGN(1.,RAY) NITI057
12 ray = Dsign ( Dabs(ray) , ray)
NTYP=2.d0+FIELD*RAY NITI058
GO TO (13,14,15), NTYP NITI059
13 MFLD(1)='EXTRAORD' NITI060
MFLD(2)='INARY' NITI061
GO TO 16 NITI062
14 MFLD(1)='NO FIELD' NITI063
MFLD(2)=' ' NITI064
GO TO 16 NITI065
15 MFLD(1)='ORDINARY' NITI066
MFLD(2)=' ' NITI067
16 MODSAV=MODX(2) NITI068
IF (PERT.EQ.0.d0) MODX(2)=' ' NITI069
IF (RAYSET.NE.0.d0) write(7,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,7d10.3)) NITI073
write(2,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
c
c first call to ELECTX to read in the array P for DUD_3D
c other progs require a sensible value of r(1) [at least]
r(1) = w(2)
call electx
c
write(2,1050) NITI077
1050 FORMAT (85H INITIAL VALUES FOR THE W ARRAY -- ALL ANGLES IN RADIANNITI078
1S, ONLY NONZERO VALUES PRINTED/) NITI079
nnw = 0
DO 17 NW=1,400 NITI080
c IF (W(NW).NE.0.d0) write(2,1700) NW,W(NW) NITI081
if(w(nw).eq.0.d0) go to 17
nnw = nnw + 1
if(nnw.gt.100) stop 'NNW too large in MAIN'
ww(nnw) = w(nw)
nww(nnw) = nw
c1700 FORMAT (I4,E19.11) NITI082
17 CONTINUE NITI083
write(2,1700) (nww(i),ww(i),i=1,nnw)
1700 format(8(i4,d11.3))
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.d0) 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=dSIN (PLAT) NITI102
CP=dSIN (PID2-PLAT) NITI103
SDPH=dSIN (TLON-PLON) NITI104
CDPH=dSIN (PID2-(TLON-PLON)) NITI105
SL=dSIN (TLAT) NITI106
CL=dSIN (PID2-TLAT) NITI107
ALPHA=dATAN2(-SDPH*CP,-CDPH*CP*SL+SP*CL) NITI108
TH0=dACOS (CDPH*CP*CL+SP*SL) NITI109
PH0=dATAN2(SDPH*CL,CDPH*SP*CL-CP*SL) NITI110
C********* LOOP ON FREQUENCY, AZIMUTH ANGLE, AND ELEVATION ANGLE NITI111
NFREQ=1 NITI112
print*,' fstep,fbeg,fend ',fstep,fbeg,fend
IF (FSTEP.NE.0.d0) NFREQ=(FEND-FBEG)/FSTEP+1.5d0 NITI113
NAZ=1 NITI114
IF (AZSTEP.NE.0.d0) NAZ=(AZEND-AZBEG)/AZSTEP+1.5d0 NITI115
NBETA=1 NITI116
IF (ELSTEP.NE.0.d0) NBETA=(ELEND-ELBEG)/ELSTEP+1.5d0 NITI117
time0 = 0
time0 = secnds(0.0)
if(w(399).eq.0.d0) go to 2211
c trace ray at VI and equiv vert freq (SECOND PASS)
nf = 1
naz = 1
nbeta = 1
2211 continue
DO 50 NF=1,NFREQ NITI118
print*,' fbeg,fstep,f ',fbeg,fstep,f
F=FBEG+(NF-1)*FSTEP NITI119
if(w(399).eq.1.d0) f = fv
DO 45 J=1,NAZ NITI120
AZ1=AZBEG+(J-1)*AZSTEP NITI121
AZA=AZ1*DEGS NITI122
GAMMA=PI-AZ1+ALPHA NITI123
SGAMMA=dSIN (GAMMA) NITI124
CGAMMA=dSIN (PID2-GAMMA) NITI125
DO 40 I=1,NBETA NITI126
BETA=ELBEG+(I-1)*ELSTEP NITI127
if(w(399).eq.1.d0) beta = pid2
EL=BETA*DEGS NITI128
time1 = 0
time1 = secnds(0.0)
CBETA=dSIN (PID2-BETA) NITI129
R(1)=R0 NITI130
R(2)=TH0 NITI131
R(3)=PH0 NITI132
R(4)=dSIN (BETA) NITI133
R(5)=CBETA*CGAMMA NITI134
R(6)=CBETA*SGAMMA NITI135
T=0.d0 NITI136
RSTART=1.d0 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 = 0
lines = 0
write(2,1000) ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL NITI143
write(2,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
write(2,2500) EL NITI148
2500 FORMAT (/31X,33HELEVATION ANGLE OF TRANSMISSION =,F12.6,4H DEG/) NITI149
c
c check to see if ray has gone out of the supplied grid
c 64 is justa number chosen at random
if(x.ge.64.d0) print*,' Ray has gone outside the supplied grid '
if(x.ge.64.d0) write(2,4323)
4323 format(/3x,' Ray has gone outside the supplied grid'/)
if(x.ge.64.d0) go to 40
c
IF (N2.GT.0.d0) GO TO 30 NITI150
CALL ELECTX NITI151
FN=dSIGN (dSQRT (dABS (X))*F,X) NITI152
write(2,2900) FN NITI153
2900 FORMAT (58H0TRANSMITTER IN EVANESCENT REGION, TRANSMISSION IMPOSSINITI154
1BLE/20H0PLASMA FREQUENCY = ,E17.10) NITI155
GO TO 44 NITI156
30 MU=dSQRT (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
print*,' IHOP = ', ihop
ccc OSEC=SECOND NITI163
ccc SECOND=KLOCK(0)*.001 NITI164
ccc DIFF=SECOND-OSEC NITI165
ccc PRINT 3500, DIFF NITI166
3500 FORMAT (36X,26HTHIS RAY CALCULATION TOOK ,F8.3,4H SEC) NITI167
time2 = 0.
time2 = secnds(time1)
print*,' Computation time for this ray in secs ',time2
write(2,5555) time2
c last position of ray
thlast = r(2) * degs
phlast = r(3) * degs
print*,' freq, elevation, range, group path, apogee ',
a f,el,range,t,apogee
print*,' Last coords of ray ',thlast,phlast
IF (PENET.AND.ONLY.NE.0..AND.IHOP.EQ.1) GO TO 44 NITI168
print*,' freq,elevation ', f,el
c
c calculate the equivalent vertical frequency
term = ((earthr * dcos (beta))/(earthr + apogee) )**2
term = dsqrt(1.d0 - term)
fv = f * term
c
c calculate the height of the equivalent mirror (no tilt)
hmirr = (earthr*dcos(beta))/dcos(beta + range/2./earthr)
a - earthr
if(hmirr.le.0.d0) hmirr = 0.d0
print*,' elev, equiv freq & mirror height ',el,fv,hmirr
write(2,2222) el,f,apogee,range,fv,hmirr
2222 format(5x,'Elev,f,apogee,range,fv,hmirr',2f6.2,2f6.1,f6.2,f6.1)
c
c find geographic location of last ground reflection point
if(.not.there) go to 200
ccc print*,' the height should be 6370 - it is ',r(1)
stlast = sngl (90.d0 - r(2) * degs)
splast = sngl (r(3) * degs)
sthl = stlast
sphl = splast
c if the position is geomagnetic, change it to geographic
if(ray.ne.0.d0) call ggm (1,sphl,sthl,splast,stlast)
write(2,201) sthl,sphl
write(6,201) sthl,sphl
201 format(/5x,'Ray last hit the ground at geogc lat & lon',
a 2f8.2)
c find the bearing from the TX to this point
slat = sngl(tlat) * degs
slon = sngl(tlon) * degs
call distbear ( slat,slon,sthl,sphl,sdist,sbear)
write(2,202) sdist,sbear
write(6,202) sdist,sbear
202 format(/5x,'Range & bearing from TX (SKYLOC) to this point',
a 2f8.1)
write(2,203) aza
write(6,203) aza
203 format(/5x,'Initial azimuth was',f8.1)
ccc print 4444, sbear,sdist,sthl,sphl
write(4,4444) sbear,sdist,sthl,sphl
4444 format(f10.1,f10.0,2f10.1)
200 if(.not.there) range = -1.d0
zero = 0.d0
if(.not.there) write(4,4444) zero,range,zero,zero
c
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.d0) write(7,5000) NITI175
5000 FORMAT (78X,1H-) NITI176
timef = 0
timef = secnds(time0)
print*,' Computation time for all rays in secs ',timef
write(2,5555) timef
5555 format(//5x,'Elapsed time for raytracing in secs ',f6.2/)
MODX(2)=MODSAV NITI177
GO TO 10 NITI178
END NITI179-