C Code comes from "Prediction of Tropospheric Radio Transmission C Loss Over Irregular Terrain, A Computer Method-1968", Longley and C Rice. C C I started with the OCR code from the report and fixed the many C errors. Typing from scratch might have been faster. Corrections C were also made as given in errata. The original, uncorrected C lines remain, commented out with "C!!!". C C Mike Markowski, mike.ab3ap@gmail.com C Dec 2022 PROGRAM COMTE C PROGRAM TO DETERMINE PARAMETERS AND WRITE OUTPUT C COMMON /M/F,D,NS,A,DH,DHS,S,E,POL,KM COMMON /MP/ H1E,H2E,H1G,H2G,DLS1,DLS2,DL1,DL2,DL,DLS,TE1,TE2,TE,KL COMMON /MLDS/ AG,AD,AS,ACR,AED,MD,AH50,AH5,D5,MS,AES,DX,H5 COMMON /ML/ D0,D1,D01,D02,A0,A1,K1,K2,AL,ALS,A0G C!!! DIMENSION ANS(3),DKM(6),DELH(6),SD(6),SA(6) DIMENSION ANS(3),DKM(6),SD(6),SA(6) REAL NS,MD,MDO,MS,MSS,MDS,K1,K2,K3,K4,LBF DATA (ANS=290.,290.,312.) DATA (DKM=5.,10.,20.,30.,50.,80.) C!!! DATA (DELH=105.,165.,234.,315.,575.) C C CALCULATION OF INPUT PARAMETERS C S=.005 $ E=15. DO 500 IX=1,3 NS=ANS(IX) A=6370./(1.-.04665*EXPF(.005577*NS)) WRITE ( 2,56)1 56 FORMAT (R1) IF (IX .EQ. 1) WRITE ( 2,57) IF (IX .EQ. 2) WRITE ( 2,58) IF (IX .EQ. 3) WRITE ( 2,59) 57 FORMAT (2X,*COLORADO PLAINS NS=290.*//) 58 FORMAT (2X,*COLORADO MOUNTAINS NS=290.*//) 59 FORMAT (2X,*OHIO NS=312.*//) DO 400 I=1,9 KK=0 DO 300 IZ=1,6 IF (IX .EQ. 2 .AND. IZ .EQ. 6) GO TO 300 IF (IX .EQ. 3 .AND. (IZ .EQ. 1 .OR. IZ .EQ. 6)) GO TO 300 D=DKM(IZ) DH=DHS=90. IF (IX .EQ. 2) DH=DHS=650. F=100. IF (I .GT. 6) F=50. IF (I .EQ. 9) F=20. POL=+1. IF (I .GT. 3 .AND. I .LT. 7) POL=-1. H1G=H1E=4. IF (IX .NE. 3 .AND. I .EQ. 9) H1G=H1E=3.3 IF (IX .EQ. 3 .AND. I .GT. 6) H1G=H1E=4.24 IF (IX .EQ. 3 .AND. I .EQ. 9) H1G=H1E=3.68 H2G=H2E=3. IF (I .EQ. 2 .OR. I .EQ. 5) H2G=H2E=6. IF (I .EQ. 3 .OR. I .EQ. 6) H2G=H2E=9. IF (IX .EQ. 3 .AND. I .EQ. 7) H2G=H2E=1. IF (IX .NE. 3 .AND. I .EQ. 7) H2G=H2E=.55 IF (IX .NE. 3 .AND. I .EQ. 8) H2G=H2E=1.7 IF (IX .NE. 3 .AND. I .EQ. 9) H2G=H2E=1.3 DLS1=SQRTF(.002*A*H1E) DLS2=SQRTF(.002*A*H2E) DLS=DLS1+DLS2 C!!! DL1=DLS1*EXPF(-.07*SQRTF(DH/H1E)) C!!! DL2=DLS2*EXPF(-.07*SQRTF(DH/H2E)) DL1=DLS1*EXPF(-.07*SQRTF(DH/MAX(5,H1E))) DL2=DLS2*EXPF(-.07*SQRTF(DH/MAX(5,H2E))) DL=DL1+DL2 TE1=(.00065/DLS1)*((DLS1/DL1-1.)*DH-3.077*H1E) TE2=(.00065/DLS2)*((DLS2/DL2-1.)*DH-3.077*H2E) TE=MAX1F((TE1+TE2),(-DL/A)) KK=KK+1 IF (D .GT. DLS) GO TO 40 CALL LOS SD(KK)=D SA(KK)=ACR GO TO 300 40 CALL DIFF SD(KK)=D SA(KK)=ACR IF (IX .EQ. 1 .AND. KK .NE. 6) GO TO 300 IF (IX .EQ. 2 .AND. KK .NE. 5) GO TO 300 IF (IX .EQ. 3 .AND. KK .NE. 4) GO TO 300 AE=A0G-K1*D0-K2*ALOG10(D0) C C WRITE OUTPUT C WRITE ( 2,60) F,DH,H1G,H2G,TE,DX 60 FORMAT (4X,*F=*F6.1,* DH=*F6.2,* H1G=*F6.2,* H2G=*F6.2, C* TE=*F10.6,* DX=*F8.2) WRITE ( 2,61) AE,K1,K2,DLS,ALS 61 FORMAT (4X,*AE=*F8.2,* K1=*F10.5,* K2=*F10.5,* DLS=*F8.2, C* ALS=*F8.2) ADX=AED+MD*DX WRITE( 2,62) AED,MD,AES,MS,ADX 62 FORMAT (4X,*AED=*F8.2,* MD=*F10.5,* AES=*F8.2,* MS=*F10.5, C* ADX=*F8.2) WRITE ( 2,63) (SD(JZ), JZ=1,KK) WRITE ( 2,64) (SA(JZ), JZ=1,KK) 63 FORMAT (4X,'D'6F10.2) 64 FORMAT (4X,'A'6F10.2) WRITE ( 2,56) 300 CONTINUE 400 CONTINUE 500 CONTINUE CALL EXIT END SUBROUTINE DIFF C SUBROUTINE TO COMPUTE DIFFRACTION ATTENUATION C COMMON /NR/ JZ,W,SW3,SW4,SA3,SA4,SAFO COMMON /MAR14/D3,D4,T5 COMMON /M/F,D,NS,A,DH,DHS,S,E,POL,KM COMMON /MP/ H1E,H2E,H1G,H2G,DLS1,DLS2,DL1,DL2,DL,DL5,TE1,TE2,TE,KL COMMON /MLDS/ AG,AD,AS,ACR,AED,MD,AH50,AH5,D5,MS,AES,DX,H5 REAL NS,MD,MDO,MS,MSS,MDS,K1,K2,K3,K4 FNA(C)=6.02+9.11*C-1.27*C*C FNB(C)=12.953+20.*ALOG10(C) FNC(C)=416.4*F**.3333333333*(1.607-C) FND(C)=(.36278/(C*F)**.3333333333)*1/((E-1.)**2+X*X)**.25 FNE(C)=C*SQRTF(E*E+X*X) RDL=DL KK=0 10 KK=KK+1 D3=DL+.5*(A*A/F)**.3333333333 IF (D3 .LT. DLS) D3=DLS D4=D3+(A*A/F)**.3333333333 T3=TE+D3/A T4=TE+D4/A C C CALCULATION OF KNIFE EDGE DIFFRACTION C V13=1.2915*T3*SQRTF(F*DL1*(D3-DL)/(D3-DL2)) V23=1.2915*T3*SQRTF(F*DL2*(D3-DL)/(D3-DL1)) V14=1.2915*T4*SQRTF(F*DL1*(D4-DL)/(D4-DL2)) V24=1.2915*T4*SQRTF(F*DL2*(D4-DL)/(D4-DL1)) AV13=FNA(V13) IF(V13 .GT. 2.4) AV13=FNB(V13) AV23=FNA(V23) IF(V23.GT.2.4) AV23=FNB(V23) AV14=FNA(V14) IF(V14.GT.2.4) AV14=FNB(V14) AV24=FNA(V24) IF(V24.GT.2.4) AV24=FNB(V24) AK3=AV13+AV23 AK4=AV14+AV24 C C CALCULATION OF ROUNDED EARTH DIFFRACTION C A1=DL1*DL1/(.002*H1E) A2=DL2*DL2/(.002*H2E) A3=(D3-DL)/T3 A4=(D4-DL)/T4 X=18000.*S/F K1=FND(A1) K2=FND(A2) K3=FND(A3) K4=FND(A4) IF (POL .EQ. -1.) GO TO 15 K1=FNE(K1) K2=FNE(K2) K3=FNE(K3) K4=FNE(K4) 15 B1=FNC(K1) B2=FNC(K2) B3=FNC(K3) B4=FNC(K4) X1=B1*DL1/A1**.6666666666 X2=B2*DL2/A2**.66666666666 X3=B3*(D3-DL)/A3**.666666666+X1+X2 X4=B4*(D4-DL)/A4**.666666666+X1+X2 XL1=450./ABSF(ALOG10(K1)**3) XL2=450./ABSF(ALOG10(K2)**3) IF(X1.GT.0..AND.X1.LE. 200. .AND.K1 .GE. 0..AND. K1.LE. .00001) C16,17 16 T=40.*ALOG10(X1)-117. T1=-117. T2=MIN1F((ABSF(T)),(ABSF(T1))) FX1=T IF (T2 .EQ. ABSF(T1)) FX1=T1 17 IF(X2.GT.0..AND. X2.LE.200. .AND. K2 .GE.0..AND. K2.LE. .00001) C18,19 18 T=40.*ALOG10(X2)-117. T1=-117. T2=MIN1F((ABSF(T)),(ABSF(T1))) FX2=T IF (T2 .EQ. ABSF(T1)) FX2=T1 19 IF(X1 .GT.0. .AND.X1 .LE.200. .AND.K1 .GT. .00001 .AND. K1 .LT. 1. C) 21,22 21 FX1=40.*ALOG10(X1)-117. IF (X1 .LE. XL1) FX1=20.*ALOG10(K1)+2.5*1. E-5*X1*X1/K1-15. C!!!22 IF(X2 .GT.0. .AND.X2 .LE.200 .AND.K2 .GT. .00001 .AND. K1 .LT. 1. 22 IF(X2 .GT.0. .AND.X2 .LE.200 .AND.K2 .GT. .00001 .AND. K2 .LT. 1. C) 23,24 23 FX2=40.*ALOG10(X2)-117. IF (X2 .LE. XL2) FX2=20.*ALOG10(K2)+2.5*1. E-5*X2*X2/K2-15. 24 W1=.0134*X1*EXPF(-.005*X1) $ W2=.0134*X2*EXPF(-.005*X2) IF(X1.GT.200. .AND. X1.LE.2000.) C FX1=W1*(40.*ALOG10(X1)-117.)+(1.-W1)*(.05751*X1-10.*ALOG10(X1)) IF(X2 .GT. 200 .AND. X2 .LE. 2000.) CFX2=W2*(40.*ALOG10(X2)-117.)+(1.-W2)*(.05751*X2-10.*ALOG10(X2)) IF(X1.GT. 2000.) FX1=.05751*X1-10.*ALOG10(X1) IF(X2 .GT. 2000.) FX2=.05751*X2-10.*ALOG10(X2) GX3=.05751*X3-10.*ALOG10(X3) GX4=.05751*X4-10.*ALOG10(X4) AR3=GX3-FX1-FX2-20. AR4=GX4-FX1-FX2-20. C C COMBINATION OF ROUNDED EARTH AND KNIFE EDGE DIFFRACTION C 28 DHD3=DH*(1.-.8 *EXPF(-.02*D3)) DHD4=DH*(1.-.8 *EXPF(-.02*D4)) P13=SQRTF((H1E*H2E)/(H1G*H2G))+(A*TE+DL)/D3 P14=SQRTF((H1E*H2E)/(H1G*H2G))+(A*TE+DL)/D4 DOL3=MIN1F(1000.,(DHD3*F/299.7925)) DOL4=MIN1F(1000.,(DHD4*F/299.7925)) 63 W3=1./(1.+.1 *SQRTF(DOL3*P13)) W4=1./(1.+.1 *SQRTF(DOL4*P14)) 31 A3=(1.-W3)*AK3+W3*AR3 A4=(1.-W4)*AK4+W4*AR4 MD=(A4-A3)/(D4-D3) AED=A4-MD*D4 DHDLS=DH*(1.-.8*EXPF(-.02*DLS)) SHDLS=.78*DHDLS*EXPF(-.5*(DHDLS**.25)) AFO=5.*ALOG10(1.+H1G*H2G*F*SHDLS*.00001) AFO=MIN1F(AF0,15.) 91 AED=AED+AFO IF (KM .EQ. 2) GO TO 40 IF (KK.EQ. 2) GO TO 20 SAFO=AFO SW3=W3 $ SW4=W4 $ SA3=A3 $ SA4=A4 29 AD=AED+MD*D 32 TD=(TE+D/A)*D SDL1=DL1 $ SDL2=DL2 $ SDL=DL $ STE1=TE1 $ STE2=TE2 STE=TE C C CALCULATION OF SCATTER ATTENUATION C CALL SCATT IF (H5 .LE. 10.) AES=AH5 -MS*D5 IF (H5 .LE. 10.) AS=AES+MS*D IF(H5.LE.10.) GO TO 30 MDS=MD AEDS=AED DH=0. GO TO 10 20 ADO=AED MDO=MD MD=MOS AED=AEDS DX1=(AH50-MS*D5-ADO)/(MD0-MS) DX2=(RDL+.25*(A*A/F)**.3333333333*ALOG10(F)) DXO=DX1*(3.-.2*H5)+DX2*(.2*H5-2.) AXO=ADO+MDO*DXO ASX=AXO+(AH5-AH50) AES=ASX-MS*DXO AS=AES+MS*D 30 DX=(AES-AED)/(MD-MS) DXN=DL+.25*(A*A/F)**.3333333333*ALOG10(F) C!!! IF (DXN .GT. DX) AES=AED+(MD-MS)*DX IF (DXN .GT. DX) AES=AED+(MD-MS)*DXN IF (DXN .GT. DX) AS=AES+MS*D IF (DXN .GT. DX) DX=DXN ACR=AD IF (D .GT. DX) ACR=AS DL1=SDL1 $ DL2=SDL2 $ DL=SDL $ TE1=STE1 $ TE2=STE2 TE=STE DH=DHS 41 CONTINUE 40 RETURN END SUBROUTINE SCATT C SUBROUTINE TO COMPUTE SCATTER PARAMETERS C COMMON /M/F,D,NS,A,DH,DHS,S,E,POL,KM COMMON /MP/ H1E,H2E,H1G,H2G,DLS1,DLS2,DL1,DL2,DL,DLS,TE1,TE2,TE,KL COMMON /MAR14/D3,D4,T5 COMMON /MLDS/ AG,AD,AS,ACR,AED,MD,AH50,AH5,D5,MS,AES,DX,H5 REAL NS,MD,MDO,MS,MSS,MDS,K1,K2,K3,K4 KK=0 10 KK=KK+1 D5=DL+ 200. D6=DL+ 400. 11 T5=TE+D5/A T6=TE+D6/A H5=MIN1F(((1./H1E +1./H2E)/(T5*F*ABSF(.007-.058*T5))), (15.)) H6=MIN1F(((1./H1E+1./H2E)/(T6*F*ABSF(.007-.058*T6))),(15.)) S5=H5+10.*ALOG10(F*T5**4)-.1*(NS-301.)*EXPF(-T5*D5/40.) S6=H6+10.*ALOG10(F*T6**4)-.1*(NS-301.)*EXPF(-T6*D6/40.) IF(T5*D5 .LE.10.) AH5=S5+103.4+.332*T5*D5-10.*ALOG10(T5*D5) IF(T6*D6 .LE.10.)AH6=S6+103.4+.332*T6*D6-10.*ALOG10(T6*D6) IF(T5*D5 .GT. 10. .AND.T5*D5.LE. 70.) AH5=S5+97.1+.212*T5*D5-2.5* CALOG10(T5*D5) IF(T6*D6 .GT. 10. .AND.T6*D6 .LE. 70.) AH6=S6+97.1+.212*T6*D6-2.5 C*ALOG10(T6*D6) IF(T5*D5 .GT. 70.) AH5=S5+86.8+.157*T5*D5+5.*ALOG10(T5*D5) IF(T6*D6 .GT. 70.) AH6=S6+86.8+.157*T6*D6+5.*ALOG10(T6*D6) MS=(AH6-AH5)/(D6-D5) IF (KK .EQ. 2) GO TO 25 IF (H5 .EQ. 10) GO TO 30 IF (KK .EQ. 1) GO TO 20 25 MS=MSS AH50=AH5 AH5=AH5S D5=D5S GO TO 30 20 DH=0. C!!! DL1=DLS1*EXPF(-.07*SQRTF(DH/H1E)) C!!! DL2=DLS2*EXPF(-.07*SQRTF(DH/H2E)) DL1=DLS1*EXPF(-.07*SQRTF(DH/MAX(5,H1E))) DL2=DLS2*EXPF(-.07*SQRTF(DH/MAX(5,H2E))) DL=DL1+DL2 C!!! TE1=(.00065/DLS1)*((DLS1/DL1-1.)*DH-3.077*H1G) C!!! TE2=(.00065/DLS2)*((DLS2/DL2-1.)*DH-3.077*H2G) TE1=(.00065/DLS1)*((DLS1/DL1-1.)*DH-3.077*H1E) TE2=(.00065/DLS2)*((DLS2/DL2-1.)*DH-3.077*H2E) TE=MAX1F((TE1+TE2),(-DL/A)) T=TE+D/A AH5S=AH5 MSS=MS D5S=D5 GO TO 10 30 CONTINUE RETURN END SUBROUTINE LOS C SUBROUTINE TO COMPUTE LINE OF SIGHT ATTENUATION C COMMON /M/F,D,NS,A,DH,DHS,S,E,POL,KM COMMON /NR/ JZ,W,SW3,SW4,SA3,SA4,SAFO COMMON /MP/ H1E,H2E,H1G,H2G,DLS1,DLS2,DL1,DL2,DL,DLS,TE1,TE2,TE,KL COMMON /MLDS/ AG,AD,AS,ACR,AED,MD,AH50,AH5,D5,MS,AES,DX,H5 COMMON /ML/ D0,D1,D01,D02,A0,A1,K1,K2,AL,ALS,A0G REAL NS,MD,MDO,MS,MSS,MDS,K1,K2,K3,K4,M KM=2 CALL DIFF KM=0 C C CALCULATION OF TWO RAY THEORY C D01=.00004*H1E*H2E*F D02=MIN1F((-AED/MD),(DL-2.)) IF (AED .GE. 0.) DO=MIN1F(D01,(.5*DL)) IF (AED .LT. 0.) DO=.5*DL IF (AED .LT. 0. .AND. D02 .GE. D0) D0=D02 D1=D0+.25*(DL-D0) IF (D1 .LE. D0) D1=D0+.25*(DLS-D0) J=0 $ DS=D IF (J.EQ. 0) D=D0 22 DIV=1. $ PSI=ATANF((H1E+H2E)/(1000.*D)) 2 DHD=DH*(1.-.8 *EXPF(-.02*D)) SH=.78*DHD*EXPF(-.5*DHD**.25) SP=SINF(PSI) X=18000.*S/F P2=(SQRTF((E-COSF(PSI)*COSF(PSI))**2+X*X)+E-COSF(PSI)*COSF(PSI))/2 C. P=SQRTF(P2) Q=X/(2.*P) IF (POL .EQ. 1.)B=(E*E+X*X)/(P2+Q*Q) IF(POL .EQ. -1.)B=1./(P2+Q*Q) IF (POL .EQ. 1.)M=2.*(P*E+Q*X)/(P2+Q*Q) IF (POL .EQ. -1.)M=2.*P/(P2+Q*Q) R2=(1.+B*SP*SP-M*SP)/(1.+B*SP*SP+M*SP) RE=SQRTF(SP) SQEXF=SQRTF(R2)*EXPF(-.0209584473*F*SH*SP)*DIV IF (SQEXF .GT. .5 .AND. SQEXF .GT. RE) RE=SQEXF C=ATANF(Q/(P+SP))-ATANF(Q/(P-SP)) IF (POL .EQ. -1.)GO TO 40 Y1=(X*SP+Q)/(E*SP+P) Y2=(X*SP-Q)/(E*SP-P) IF (E*SP .GE. P) C=ATANF(Y1)-ATANF(Y2)+3.141592654 IF (E*SP .LT. P .AND. P*SP .GT. .5) C=ATANF(Y1)+ATANF(Y2) IF (E*SP .LT. P .AND. P*SP .LE. .5) C=ATANF(Y1)-ATANF(Y2) 40 IF (J .EQ. 0) CA0=-10.*ALOG10(1.+RE*RE-2.*RE*COSF(.000041917*F*H1E*H2E/D0-C)) IF (J .EQ. 1) GO TO 3 D=D1 $ J=1 GO TO 22 3 CONTINUE D=DS IF (J .EQ. 1) CA1=-10.*ALOG10(1.+RE*RE-2.*RE*COSF(.000041917*F*H1E*H2E/D1-C)) C C COMBINATION OF TWO RAY THEORY AND DIFFRACTION C ALS=AED+MD*DLS AL=AED+MD*DL DED0=AED+MD*D0 DED1=AED+MD*D1 SA0=A0 $ SA1=A1 C!!! W=(MIN1F(H1E,H2E)/MAX1F(H1E,H2E))/(1.+F*DH*.0001) W=1./(1.+F*DH*.0001) A0=MIN1F((W*A0+(1.-W)*DED0),DED0) A1=MIN1F((W*A1+(1.-W)*DED1),DED1) 10 K2=((ALS-A0)*(D1-D0)-(A1-A0)*(DLS-D0))/((D1-D0)*ALOG10(DLS/D0)- C(DLS-D0)*ALOG10(D1/D0)) K2=MAX1F(K2,0.) K1=((ALS-A0)-K2*ALOG10(DLS/D0))/(DLS-D0) IF (K1 .GE. 0.) GO TO 50 K1=0. K2=(ALS-A0)/(ALOG10(DLS/D0)) 50 AG=A0+K1*(D-D0)+K2*ALOG10(D/D0) IF (AG .LT. 0.) AG=0. 51 A0G=A0 53 ACR=AG RETURN END