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 16-dec-86 c C SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE NITI002 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 /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*8 N2,N2I,LOGTEN,K,MAXSTP,INTYP,MAXERR,MU NITI021 COMPLEX*16 PNP,POLAR,LPOLAR NITI022 c c PRINT writes to unit 6 open (unit=6,file='data.out',status='old') c PUNCH writes to unit 7 open (unit=7,file='punch.out',status='old') c normal input is from unit 5 -- read by s/r readW open (unit=5,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 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=dLOG(10.d0) 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 IF (SKIP.EQ.0.) SKIP=MAXSTP NITI056 c 12 RAY=dSIGN(1.,RAY) NITI057 12 ray = Dsign ( Dabs(ray) , ray) NTYP=2.+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.) MODX(2)=' ' NITI069 IF (RAYSET.NE.0.) 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,7E10.3)) NITI073 write(6,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 write(6,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.) write(6,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=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 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 call timer(iticks) 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=dSIN (GAMMA) NITI124 CGAMMA=dSIN (PID2-GAMMA) NITI125 DO 40 I=1,NBETA NITI126 BETA=ELBEG+(I-1)*ELSTEP NITI127 EL=BETA*DEGS NITI128 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. 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 = 0 lines = 0 write(6,1000) ID,NDATE,MODX,MODY,MODZ,MODRIN,MFLD,KOLL NITI143 write(6,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(6,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=dSIGN (dSQRT (dABS (X))*F,X) NITI152 write(6,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 ccc OSEC=SECOND NITI163 ccc SECOND=KLOCK(0)*.001 NITI164 ccc DIFF=SECOND-OSEC NITI165 ccc PRINT 3500, DIFF call timer (jticks) itime = ( jticks - iticks )/100 print*,'freq,elevation ', f,el print*,' elapsed time in seconds ... ',itime 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.) write(7,5000) NITI175 5000 FORMAT (78X,1H-) NITI176 MODX(2)=MODSAV NITI177 GO TO 10 NITI178 END NITI179-