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-