SUBROUTINE DIFFRACT(IELEVAT, IDIST, FREQ, HTAMSL, DLOSS, CLUTTER) C***************************************************************** C Subroutine to determine diffraction loss C Inputs: IELEVAT Array of Elevations along path C IDIST Index to Distance of point under study C Outputs: DLOSS Diffraction Loss at point on the path C CLUTTER Clutter Loss at point on the path C C 5/17/1999 C C**************************************************************** real ielevat(*) dd = 0.1 EK = 1.333 ! Earth curvature (4/3 earth) HR = 9.1 ! Rcvr Ant ht (m), for 30 feet HRAMSL = IELEVAT(IDIST) + HR DIST = float(idist)*dd FRESMIN = HR + 1.0 DLOSS = 0.0 TDLOSS = 0.0 RDLOSS = 0.0 ADJ = 0.0 C************************************************************ C Determine clearance ratio C************************************************************ OBSTACLE = 1.0 ! Re-initialize for next azimuth RATIO = 1.0 clearmin = 1000 ! Minimum clearance (hilltop) SLOPE = (HRAMSL - HTAMSL)/ DIST DO 110 L = 1, IDIST - 1 XD = L * dd EARTHB = .0785 * XD * (DIST - XD) / EK ! Earth bulge term CEL = IELEVAT(L) + EARTHB ! Elev, incl curved earth HLOS = HTAMSL + SLOPE * XD CLEARANS = HLOS - CEL FRESNEL = 548. * SQRT(XD * (DIST - XD) / (DIST * FREQ)) IF (CLEARANS .GT. 0.0 ) go to 109 IF (FRESNEL .LT. FRESMIN) FRESNEL = FRESMIN 109 continue rat = clearans / fresnel if(rat .GT. ratio) go to 110 ratio = rat CLen = Dist - XD if(clearans .LE. clearmin) then clearmin = clearans obstacle = L obselev = CEL end if 110 CONTINUE IF (ratio .GT. .8) obstacle = idist - 1 ObDist = obstacle * dd RLength = Dist - Obdist ! in km PLength = RLength / dd ! in point counts 120 CONTINUE C---------------------------------------- C Find smooth earth Horizon distance Elmin = HtAMSL - 30. Txsite = IELEVAT(1) Init = 1 Last = IDist DO 121 J = Init, Last D = J Y = IELEVAT(J) ! Elevation at dist IF (Y .LT. Elmin) Elmin = Y 121 CONTINUE Ref = Elmin HrZ = 4.122 * SQRT(Hr) Ht = HTAMSL - Ref IF (Ht .LT. 30) Ht = 30 HORIZON = 4.122 * SQRT(Ht) Fract = dist/horizon !Distance in terms of HORIZON c ************************************ c Beyond Horizon Adjustment C (unknown reason -- but TASO data in smooth-earth C situations indicated such need) c ************************************ adj = 0. blend = 100/SQRT(freq) IF (Ratio .GT. 0.2) go to 125 freqfact = 4.5 - 25.57/(Log10(freq) + 5.53) beyond = dist - horizon IF (beyond+blend .LT. 0) go to 125 Horzfact = 224 /(horizon + 70) - .6 IF (Horzfact .LT. .2) Horzfact = .2 IF (Horzfact .GT. 1.8) Horzfact = 1.8 beyond = Horzfact * beyond beyfact = 11.7 - 222.3/(beyond + 19) IF (beyond .GT. 20) go to 124 blend = 100.0 / SQRT(freq) X1 = -blend ! blending near horizon X2 = 20 Y2 = 6.0 bslope = Y2 /(X2 - X1) beyfact = bslope * (beyond - X1) 124 Adj = beyfact * freqfact IF (Adj .LT. 0.0) Adj = 0.0 125 continue C ************************************************** C Smooth-earth & Knife-edge diffraction losses C ************************************************** EDGE = -1.377 * RATIO **2 + 11.31 * RATIO - 6. IF (RATIO .LT. -.5) EDGE = 50.4/(1.6 - RATIO) - 36.0 IF (EDGE .GT. 0.) EDGE = 0. SMOOTH = 38.68 * RATIO - 21.66 IF (ratio .LT. -0.8) smooth = 3030/(11.2-ratio) - 305.7 IF (SMOOTH .GT. EDGE) SMOOTH = EDGE IF (RATIO .GE. .57) go to 300 ! free space condition ! no diffraction loss C ********************************** C Varying terrain (non isolated) C *********************************** INIT = OBSTACLE - 100 !limit to 10 km LAST = OBSTACLE + 100 IF (INIT .LT. 1) INIT = 1 IF (LAST .GE. IDIST) LAST = IDIST - 1 SUMX = 0. SUMX2 = 0. SUMY = 0. SUMXY = 0. N = LAST - INIT + 1 DO 130 L = INIT, LAST X = L Y = IELEVAT(L) SUMX = SUMX + X SUMX2 = SUMX2 + X * X SUMY = SUMY + Y SUMXY = SUMXY + X * Y 130 CONTINUE DEN = N * SUMX2 - SUMX * SUMX TSLOPE = (N * SUMXY - SUMX * SUMY) / DEN TB = (SUMY * SUMX2 - SUMX * SUMXY) / DEN c Determine fit through MEDIAN (translate fitted line) c ---------------------------------------- MPTS = (N+1) / 2 BOUNDUP = 1000 BOUNDLOW = -1000 TEST = 0 GO TO 145 !start at mean value 141 BOUNDUP = TEST TEST = (BOUNDUP + BOUNDLOW) / 2 IF (BOUNDLOW .LE. -1000) TEST = BOUNDUP - 10 GO TO 145 142 BOUNDLOW = TEST TEST = (BOUNDUP + BOUNDLOW) / 2 IF (BOUNDUP .GE. 1000) TEST = BOUNDLOW + 10 145 COUNTL = 0 COUNTH = 0 DO 150 L = INIT, LAST X = L Y = IELEVAT(L) FIT = TSLOPE * X + TB + test DIFF = Y - FIT IF (DIFF .GT. TEST) COUNTL = COUNTL + 1 IF (DIFF .GE. TEST) COUNTH = COUNTH + 1 150 continue IF (boundup - boundlow .LE. 1) go to 151 IF (COUNTH .LT. MPTS) GO TO 141 IF (COUNTL .GT. MPTS) GO TO 142 151 TRANS = TEST C Determine Delta H C -------------------- SUMD = 0. SUML2 = 0. SUMH2 = 0. TMIN = 0. TMAX = 0. DO 160 L = INIT, LAST X = L Y = IELEVAT(L) FIT = TSLOPE * X + TB + Trans DIFF = Y - FIT IF (DIFF .LT. TMIN) TMIN = DIFF IF (DIFF .GT. TMAX) TMAX = DIFF SUMD = SUMD + DIFF IF (DIFF .GE. 0.) then ! Separate the upper part SUMH2 = SUMH2 + DIFF * DIFF ! from the lower part ELSE SUML2 = SUML2 + DIFF * DIFF END IF 160 CONTINUE SUMD2 = SUMH2 + SUML2 HILL = 1.6 * TMAX AVG = SUMD / N VARIAN = SUMD2 / N - AVG * AVG if(varian .LT. 0.0) varian = 0.0 SIGMA = SQRT(VARIAN) delta = 2 * 1.282 * SIGMA ! 90% Probability (two tails) DELTAH = delta IF (SUMD2 .LE. .1) go to 165 ! terrain is flat IF (delta .GT. HILL) THEN DELTAH = HILL ELSE DELTAH = Hill/2 + delta * SUML2/SUMD2 END IF 165 CONTINUE TBlock = -DeltaH c ---------------------------------- C Receiver very near Obstacle IF (RLength .GT. 5) go to 170 diff = Obselev - Ielevat(IDist) IF (diff .LE. 0) go to 170 proj = .35 + 3/(RLength + .4) Cliff = diff * proj IF (DeltaH .GT. Cliff) DeltaH = (DeltaH + Cliff) / 2 170 continue C --------------------------------- C Short paths (no secondary obstacles & line-of-sight) IF (ratio .LT. 0.0) go to 175 IF (Fract .LT. .3) then DELTAH = DELTAH + 100*(1.125 + .0422/(Fract-.3375)) END IF 175 continue C -------------------- Flim = 3413.8/(DeltaH + 266.7) - 5.8 IF (CLen .GE. Flim) go to 180 IF (Ratio .GE. 0) go to 180 Plow = Dist/2 IF (Dist .GT. 2*Flim) Plow = Flim F5 = 548. * SQRT(Plow * (DIST - Plow)/(DIST * FREQ)) PrimLim = -DeltaH / F5 IF (Ratio .LT. PrimLim) Ratio = PrimLim EDGE = -1.377 * RATIO **2 + 11.31 * RATIO - 6. IF (RATIO .LT. -.5) EDGE = 50.4/(1.6 - RATIO) - 36.0 IF (EDGE .GT. 0.) EDGE = 0. SMOOTH = 38.68 * RATIO - 21.66 IF (ratio .LT. -0.8) smooth = 3030/(11.2-ratio) - 305.7 IF (SMOOTH .GT. EDGE) SMOOTH = EDGE 180 Continue C -------------------------- ROUND = 75. / (DELTAH + 75.) ! Roundness factor IF (ROUND .LT. 0.) ROUND = 0. IF (ROUND .GT. 1.) ROUND = 1. RANGE = EDGE - SMOOTH DLOSS = EDGE - Round * RANGE c-------------------------------------- c Adjust roundness near horizon Pround = Round If (Ratio .LT. -1.1) go to 190 Part = fract IF (Part .GT. .7) then Power = .17 + .166 /(fract - .5) IF (Power .GT. 1 ) Power = 1 Diff = 1 - Power Rats = 1.33 - .266/(Ratio + 1.3) IF (Rats .LT. 0) Rats = 0 IF (Rats .GT. 1) Rats = 1 Power = 1 - Rats*Diff power = power * log10(round) Pround = 10**power DLOSS = EDGE - PRound * RANGE END IF 190 Continue c ----------------- ptlim = 5.0 * Round ! set dist to look for 2nd obs. ptlim = ptlim/dd ! in point counts IF (ptlim .LT. 10.) ptlim = 10.0 C======================================= C Beyond line-of-sight TRound = PRound Round2 = PRound IF (Round .GE. .8) go to 260 IF (ratio .GT. 0.0) go to 260 TAdj = 0. TOBS = 1.0 TRATIO = .57 clearmin = 1000 SLOPE = (Obselev - HTAMSL)/ Obdist Lastpt = obstacle - ptlim DO 210 L = 10, Lastpt XD = L * dd EARTHB = .0785 * XD * (DIST - XD) / EK ! Earth bulge term CEL = IELEVAT(L) + EARTHB ! Elev, incl curved earth HLOS = HTAMSL + SLOPE * XD CLEARANS = HLOS - CEL FRESNEL = 548. * SQRT(XD * (Obdist - XD) / (Obdist * FREQ)) IF (CLEARANS .GT. 0.0 ) go to 205 IF (FRESNEL .LT. FRESMIN) FRESNEL = FRESMIN 205 continue rat = clearans / fresnel if(rat .GT. Tratio) go to 210 Tratio = rat if(clearans .LE. clearmin) then clearmin = clearans Tobs = L end if 210 CONTINUE ARatio = TRatio - .3 IF (ARatio .GT. 0.0) then TFactor = 4.95 /(ARatio+4.5) - 1.1 else TFactor = 2 + 11/(ARatio - 5.5) end if TAdj = TFactor * (1.0 - PRound) TRound = PRound + TAdj IF (TROUND .LT. 0.) TROUND = 0. IF (TROUND .GT. 1.) TROUND = 1. 250 Continue C ! end of transmitter side c====================== c Adjusting roundness for long beyond line-of-sight paths c (loss between primary obstacle & receiver) Round2 = TRound IF (RLength .LE. 10) Go to 260 ROBS = Obstacle + 1.0 RRATIO = .57 clearmin = 1000 ! Minimum clearance (hilltop) SLOPE = (HRAMSL - obselev)/ Rlength Initpt = Obstacle + ptlim Lastpt = Idist - 5 DO 252 L = Initpt, Lastpt XD = L * dd AXD = XD - ObDist EARTHB = .0785 * XD * (DIST - XD) / EK ! Earth bulge term CEL = IELEVAT(L) + EARTHB ! Elev, incl curved earth HLOS = obselev + SLOPE * AXD CLEARANS = HLOS - CEL FRESNEL = 548.* SQRT(AXD* (Rlength-AXD)/(Rlength * FREQ)) IF (CLEARANS .GT. 0.0 ) go to 251 IF (FRESNEL .LT. FRESMIN) FRESNEL = FRESMIN 251 continue rat = clearans / fresnel IF (rat .GT. RRatio) go to 252 Rratio = rat Cnum = clearans CLen = Rlength - AXD IF (clearans .LE. clearmin) then clearmin = clearans Robs = L end if 252 CONTINUE IF (CLen .GE. 5) go to 255 IF (RRatio .GE. 0) go to 255 Plow = Rlength/2 IF (Rlength .GT. 10) Plow = 5 F5 = 548. * SQRT(Plow *(Rlength - Plow)/(Rlength * freq)) RxLim = Cnum/F5 IF (RRatio .LT. RxLim) RRatio = RxLim 255 Continue c ---- undo 2nd obstacle gain TX side ---- un = 1 IF (Round .GT. .5) go to 258 !smooth terrain IF (Fract .LT. .8) go to 259 IF (RRatio .GT. .1) go to 259 IF (RLength .LT. 15.) go to 259 un = 1. - 10. * RRatio IF (un .GT. 1.0) un = 1.0 258 IF (TAdj .GT. 0) go to 259 TRound = TRound - un * TAdj 259 Continue c ======================== ARatio = RRatio - .5 IF (ARatio .GT. 0.0) then RFactor = 4.95 /(ARatio+4.5) - 1.1 else RFactor = 2 + 11/(ARatio - 5.5) end if RAdj = RFactor * (1.0 - Round) Round2 = TRound + Radj ! adjusted for both sides IF (ROUND2 .LT. 0.) ROUND2 = 0. IF (ROUND2 .GT. 1.) ROUND2 = 1. 260 Continue c------------------------------- c terrain variations near receiver (limit smoothness) Second = Round IF (Round .LT. 0.4) go to 295 ! skip -- rough terrain IF (Dist/Horizon .LT. .8) go to 295 IF (Rlength .LT. 10) go to 285 261 Continue INIT = IDist - 200 IF (INIT .LT. obstacle + 50) INIT = obstacle + 50 LAST = IDist IF (Last-Init .LT. 50) go to 285 SUMX = 0. SUMX2 = 0. SUMY = 0. SUMXY = 0. N = LAST - INIT + 1 DO 265 L = INIT, LAST X = L Y = IELEVAT(L) SUMX = SUMX + X SUMX2 = SUMX2 + X * X SUMY = SUMY + Y SUMXY = SUMXY + X * Y 265 CONTINUE DEN = N * SUMX2 - SUMX * SUMX TSLOPE = (N * SUMXY - SUMX * SUMY) / DEN TB = (SUMY * SUMX2 - SUMX * SUMXY) / DEN c Determine fit through MEDIAN (translate fitted line) c ---------------------------------------- MPTS = (N+1) / 2 BOUNDUP = 1000 BOUNDLOW = -1000 TEST = 0 GO TO 275 !start at mean value 271 BOUNDUP = TEST TEST = (BOUNDUP + BOUNDLOW) / 2 IF (BOUNDLOW .LE. -1000) TEST = BOUNDUP - 10 GO TO 275 272 BOUNDLOW = TEST TEST = (BOUNDUP + BOUNDLOW) / 2 IF (BOUNDUP .GE. 1000) TEST = BOUNDLOW + 10 275 COUNTL = 0 COUNTH = 0 DO 276 L = INIT, LAST X = L Y = IELEVAT(L) FIT = TSLOPE * X + TB + test DIFF = Y - FIT IF (DIFF .GT. TEST) COUNTL = COUNTL + 1 IF (DIFF .GE. TEST) COUNTH = COUNTH + 1 276 continue IF (boundup - boundlow .LE. 1) go to 277 IF (COUNTH .LT. MPTS) GO TO 271 IF (COUNTL .GT. MPTS) GO TO 272 277 TRANS = TEST C Determine Delta H C ------------------ SUMD = 0. SUMd2 = 0. DO 282 L = INIT, LAST X = L Y = IELEVAT(L) FIT = TSLOPE * X + TB + TRANS DIFF = Y - FIT SUMD = SUMD + DIFF SUMd2 = SUMd2 + DIFF * DIFF 282 CONTINUE AVG = SUMD / N VARIAN = SUMD2 / N - AVG * AVG IF (varian .LT. 0.0) varian = 0.0 SIGMA = SQRT(VARIAN) delta2 = 2 * 1.282 * SIGMA Second = 75. / (DELTA2 + 75.) ! Roundness near receiver 285 Continue Ref = ObDist IF (ObDist .GT. Horizon) Ref = Horizon X = dist/Ref Fact = .4 * X + .6 IF (Fact .LT. 1) Fact = 1 Slim = Fact * Second IF (Slim .GT. 1.) Slim = 1. IF (Round/Second .LT. .9) go to 295 c similar terrain throughout radial c --------------------------------- ARatio = RRatio - Ratio/4 IF (ARatio .GT. 1 ) ARatio = 1 IF (ARatio .GT. 0.0) then RFactor = 4.95 /(ARatio+4.5) - 1.1 else RFactor = 2 + 11/(ARatio - 5.5) end if RAdj = RFactor * (1.0 - Round) Round2 = round + RAdj IF (Round2 .GT. Slim) Round2 = Slim c------------------------------ 295 continue DLOSS = edge + Round2 * (smooth - edge) 300 Continue DLOSS = DLOSS - Adj C ************************************ c Clutter loss c ************************************ 350 continue C Clutter loss as function of frequency (smoothed out near LOS) CMod = 1.0 ! clutter modifier (from land-usage ?) C = 4.5 * LOG10(FREQ) - 4.0 C = CMod * C Span = .28 - .144 * (Round - 1.2) Para = C / Span**2 IF (Ratio .GT. .4) Para = 6.25 * (C - 1) CLUTTER = Para * (RATIO - .4)**2 - C IF (CLUTTER .GT. 0.) CLUTTER = 0. C Adjustment for steep angle of arrival CSlope = SQRT(freq)/350 steep = 1 + CSlope * (dist - Horizon) IF (steep .LT. 0) steep = 0 IF (steep .GT. 1) steep = 1 CLUTTER = CLUTTER * steep C +++++++++++++++++++++++++++++++++++ C Troposcattering Tropo100 = 20.34 - .077 * Dist FreqAdj = 23.978 - 58026.76 / (Freq + 2320) IF (Freq .GT. 1000) FreqAdj = 24.5 - 7200/(Freq+3000) TropoFd = Tropo100 - FreqAdj FS_field = 106.9 - 20 * LOG10(Dist) Scatter = TropoFd - FS_field !loss ref to free space DiffL = Scatter - DLOSS Combine = 150/(20 - DiffL) - 5 IF (DiffL .LT. -10) Combine = 0 IF (DiffL .GT. 10) Combine = DiffL DLOSS = DLOSS + Combine RETURN END