pe.f90

These routines model tropospheric radiowave propagation over variable terrain and calculates propaga

if( ef.lerr6 ) then
ierror = -6
return
else
tr.itp = min(tr.itp + 1, mxter)
tr.terx(tr.itp) = vnp.rmax*l.01
tr.tery(tr.itp) = tr.tery(tr.itp - 1)
end if
end if
! Test to see if the unit number for the scratch file is already attached to
! another file. If so, search for a unit number that is unattached.
inquire( iscr, OPENED = lopen )
do while (lopen)
iscr = iscr + 1
inquire( iscr, OPENED = lopen )
end do
open( iscr, status = 'SCRATCH')
write(iscr,*) tr.terx(l), tr.tery(l) ! Keep first point of
do i = 2, tr.itp-l ! terrain profile.
iml = i - 1
ipl = i + 1
xml = tr.terx(iml)
yml = tr.tery(iml)
xi = tr.terx(i)
yi = tr.tery(i)
xpl = tr.terx(ipl)
ypl = tr.tery(ipl)
dxl = amaxl( l.e-3, xi - xml)!
dx2 = amaxl( l.e-3, xpl - xi)!
sll = (yi - yml) / dxl
s12 = (ypl - yi) / dx2
scd = s12 - sll ! dx is taken to be I m
! If second derivative is large enough then keep this point.
if( abs(scd) .GT. l.e-3 ) write(iscr,*) xi, yi
end do
write(iscr,*) tr.terx(tr.itp), tr.tery(tr.itp) !Keep last point
rewind( iscr ) !in profile.
! Now the scratch file contains all the necessary points for the terrain
! profile. Go back and read them in the arrays TR.TERX), TR.TERY).
bugfix = 0.
lflag = eof(iscr)
tr.itp = 0
do while( .not. lflag)!
tr.itp = tr.itp + 1
read(iscr,*) tr.terx(tr.itp), tr.tery(tr.itp)
if( tr.terx(tr.itp) .ge. vnp.rmax ) exit
bugfix = 0.
lflag = eof(iscr)
end do
close(iscr)

! Determine minimum height of terrain profile. Then adjust entire terrain
! profile by this minimum height HMINTER such that this is the new 0 reference.
hminter = vnp.hmax
do i = 1, tr.itp
yi = tr.tery(i)
if ( yi .1t. hminter ) hminter = yi
end do
! Get maximum height of terrain, then adjust terrain elevations by height
! offset.
htermax = 0.
do i = 1, tr.itp
if( tr.tery(i) .gt. htermax ) htermax = tr.tery(i)
tr.tery(i) = tr.tery(i) - hminter
end do
! Return error code if VNP.HMAX does not exceed the maximum height of the
! terrain profile.
if( htermax .gt. vnp.hmax ) then
ierror = -8
return
end if
antref = sv.antht + tr.tery(l)
do i = 1, tr.itp-I
ipl = i + 1
yl = tr.tery(i)
xl = tr.terx(i)
y2= tr.tery(ipl)
x2 = tr.terx(ipl)
xdif = x2 - xl
ydif = y2 - yl
xdif = amaxl( xdif, l.e-5)
slope = ydif / xdif
slp(i) = slope
! Calculate angle made from each terrain point height to antenna height above
! reference (HMINTER). Determine maximum propagation angle so that direct ray
! angle will clear highest peak.
if( yl .gt. antref ) then
angle = atan( (yl-antref) / xl)
if( angle .gt. angu ) angu = angle
end if
end do
! Add 1/2 degree to the angle that clears the highest peak.
angu = angu + hdeg
end if
hmref = vnp.hmin - hminter
htlim = vnp.hmax-hminter
zlim = amaxl( htlim, antref )
! Initialize refractivity arrays.
call refinit( ef.lerrl2, vnp.rmax, rf, IERROR )

if( ierror .ne. 0 ) return
! Calculate gradients and other variables for use in ray tracing.
do i = 1, lvlep-I
rml = refdum(i)
rm2 = refdum(i+l)
hl = htdum(i)
h2 = htdum(i+l)
grad = ( rm2 - rml ) / ( h2 - hl )
if( abs( grad ) .lt. l.e-3 ) grad = sign( I., grad )*l.e-3
dmdh(i) = grad * 1.e-6 ! for ray trace formulas
end do
jls = 0
rmatht = 0.
do i = 1, lvlep-I
if((antref .1t. htdum(i+l)).and.(antref .ge. htdum(i))) then jls = i
rmatht = refdum(i) + (antref - htdum(i)) * dmdh(i)*l.e6
exit
end if
end do
rlim = .9 * vnp.rmax
prang = vnp.propang * radc
! Calculate the critical angle and perform ray trace to get the maximum
! propagation angle such that you get coverage at all heights and ranges
! >= 90% of maximum range. This is done for automatic angle calculation.
! Get minimum M-unit value of profile for all heights above transmitter height.
j = jls + 1
rmina = refdum(j)
do i = j, lvlep
if( refdum(i) .1t. rmina ) rmina = refdum(i)
end do
! Get minimum M-unit value of profile for all heights below transmitter heights.
rminb = refdum(jls)
do i = jls, 1, -1
if( refdum(i) .1t. rminb ) rminb = refdum(i)
end do
! Get critical angle if the transmitter is within or above a duct.
acrit = 0.
acritl = 0.
acrit2 = 0.
delml = rmatht - rmina
delm2 = rmatht - rminb
if( delml .gt. 0. ) acritl = sqrt( 2.e-6 * delml )
if( delm2 .gt. 0. ) acrit2 = sqrt( 2.e-6 * delm2 )
acrit = amaxl( acritl, acrit2 ) + l.e-4
thetamax = acrit
at = atan( (zlim-antref) / vnp.rmax )
thetamax = amaxl( angu, at, acrit )
! If user inputs non-zero propagation angle, use that value.
if( prang .ge. l.e-6 ) thetamax = prang
! Get THETAMAX based on shallowest reflected ray traced to reach maximum height

! and still be within 90% of maximum range (for smooth surface). For terrain case
! the direct ray angle is used.
call tracea( tr, prang, acrit )
! Add buffer for filter region.
thetamax = thetamax / .75
! Put lower limit on THETAMAX depending on frequency ( in MHz ):
! for 5000 <= SV.FREQ <= 9000, THETAMAX >= .5 deg
! for 4100 <= SV.FREQ < 5000, THETAMAX >= .6 deg
! for 2900 <= SV.FREQ < 4100, THETAMAX >= .7 deg
! for 2500 <= SV.FREQ < 2900, THETAMAX >= .8 deg
! for 1500 <= SV.FREQ < 2500, THETAMAX >= .9 deg
! for 600 < SV.FREQ < 1500, THETAMAX >= 1 deg
! for 400 < SV.FREQ <= 600, THETAMAX >= 2 deg
! for 200 < SV.FREQ <= 400, THETAMAX >= 3 deg
! for SV.FREQ <= 200, THETAMAX >= 4 deg
if( sv.freq .le. 9000. ) thetamax = amaxl(thetamax, 8.72665e-3)
if( sv.freq .1t. 5000. ) thetamax = amaxl(thetamax, 1.047197e-2)
if( sv.freq .1t. 4100. ) thetamax = amaxl(thetamax, 1.22173e-2)
if( sv.freq .1t. 2900. ) thetamax = amaxl(thetamax, 1.396263e-2)
if( sv.freq .1t. 2500. ) thetamax = amaxl(thetamax, 1.570796e-2)
if( sv.freq .1t. 1500. ) thetamax = amaxl(thetamax, 1.745329e-2)
if( sv.freq .le. 600. thetamax = amaxl(thetamax, 3.4906585e-2)
if( sv.freq .le. 400. thetamax = amaxl(thetamax, 5.2359877e-2)
if( sv.freq .le. 200. thetamax = amaxl(thetamax, 6.981317e-2)
if(( sv.polar .eq. 'V' ) .and. ( prang .le. l.e-6 ) thetamax = 2. * thetamax
! Get FFT size based on THETAMAX
call getfftsz( ZLIM )
! Maximize THETAMAX within determined FFT size for terrain cases and if
! using automatic angle calculation.
if(( fter ) .and. (prang .le. l.e-6)) then
! Use 74% of ZMAX instead of 75% to leave some slop and ensure the FFT size is
! not surpassed.
if( .74*zmax .gt. zlim ) then
thetafrac = thetalaunch / thetamax
zmax = zlim / .74
sthetamax = float(n) * wl * .5 / zmax
! Put upper limits on THETAMAX depending on frequency.
if( sv.freq .gt. 1000. ) then
sthetamax = aminl( sthetamax, sdeglO
else
sthetamax = aminl( sthetamax, sdegl5
end if
delz= wl * .5 / sthetamax
thetamax = asin( sthetamax )
zmax = float(n) * delz
thetalaunch = thetafrac * thetamax
end if
end if
! Determine horizon range based on transmitter height and 0 receiver height
! by RHOR = 3572. * sqrt( 1.3333 * antref )

rhor = 4124.5387 * sqrt( sv.antht)
dr = 2. * fko * delz**2
rkm = vnp.rmax * 1.e-3
! Determine range step.
if( fter ) then
dr = aminl( dr, 700.)
if( rkm ge. 5. ) rllim = 75.
if( rkm ge. 10. rilim = 90.
if( rkm ge. 15. rllim = 100.
if( rkm ge. 20. rllim = 110.
if( rkm ge. 30. rllim = 175.
if( rkm ge. 50. rllim = 200.
if( rkm ge. 75. rllim = 250.
if( rkm ge. 100. ) rllim = 300.
dr = amaxl( dr, rllim )
else
dr = aminl( dr, 1000. )
dr = amaxl( dr, 30. )
if( vnp.rmax .ge. rhor ) dr = amaxl( 300., dr
end if
dr2 = .5 * dr
! Path loss constant.
plcnst=20.*alogl0(2.*fko)
! Initialize variables for free-space propagator phase calculations.
delp = pi/zmax
FNorm = 2. / N
cnst = delp / fko
nml = n - 1
dz2 = 2. delz
! Initialize variables and set-up filter array.
no4 = n/4
n34 = 3.* no4
cn75 = 4.* pi / N
do i = 0, no4
fj= cn75 * float(i)
filt(i) = .5 + .5 * cos(fj)
end do
! Initialize dielectric ground constants for vertical polarization.
ig = 1
if( tr.igr .eq. 0 ) then
tr.igr = 1
tr.rgrnd(l) = 0.
tr.igrnd(l) = 0
end if
if( sv.polar .eq. 'V' ) call dieinit( sv, tr )
! Initialize starter field.
call xyinit( sv, tr )
! Transform to z-space.
call fft( u )
! Initialize C1 and C2 for start of PE calculations - only for vertical
! polarization. NOTE: THIS IS ONLY FOR SMOOTH SURFACE.

if( sv.polar .eq. 'V' ) then
clc = cmplx( 0., 0.)
c2c = cmplx( 0., 0.)
do i = 0, n
nmi = n - i
ui = u(i)
unmi = u(nmi)
if(( i .eq. 0 ) .or. (i .eq. n )) then
ui = .5 * ui
unmi = .5 * unmi
end if
iv = mod( i, 2
cx = rav(i)
if( iv .eq. 1 ) cx = -rav(i)
clc = ui * rav(i)
c2c = unmi * cx
cl = cl + clc
c2 = c2 + c2c
end do
cl = cl * rk
c2 = c2 * rk
end if
ylast = 0.
if( fter ) ylast = tr.tery(1)
ycurm = 0.
rout = 0.
ycur = 0.
! Define mesh array in height
do i=0,n
ht(i)= float(i)*delz
end do
! If smooth surface, trace THETALAUNCH ray and store all heights at each
! output range step in array HLIM().
call traceh( vnp.nrout )
! Determine the free-space propagator (p-space) arrays.
call phasel
! If smooth surface and range-independent case then initialize all refractivity
! and z-space propagator arrays now.
if( rf.nprof .eq. 1 ) call profref( hminter, 0 )
if(( .not. fter ) .and. (rf.nprof .eq. 1 )) then
call intprof
call phase2
end if
end


! ********************** SUBROUTINE ANTPAT ******************************
! Module Name: ANTPAT
! Module Security Classification: UNCLASSIFIED
! Purpose: Determines the antenna pattern factor for angle passed to routine.
! Version Number: 1.5
! INPUTS:
! Argument List: IPTRN, SANG
! Common: AFAC, BW, ELV, PELEV, SBW, UMAX
! OUTPUTS:
! Argument List: PATFAC
! Files Included: NONE
! Calling Routines: XYINIT
! Routines called: NONE
! GLOSSARY: For common variables refer to main glossary.
! Input Variables:
! IPTRN = Type of antenna pattern.
! SANG = Sine of angle for which antenna pattern is sought.
! Output Variables:
! PATFAC = Antenna pattern factor.
! Local Variables:
! ARG = Angle argument used for SINX/X and generic height-finder antenna
! pattern
! DIRANG = Sine of direct ray angle = abs(SANG)
! PR = Sine of angle U relative to sine of elevation angle (i.e.,
! sine(u) - sine(elv)
! U = Angle for which antenna pattern is sought.
! UDIF = Angle U relative to the antenna elevation angle (i.e., U-ELV)
subroutine antpat( iptrn, sang, PATFAC )
common / pattern / pelev, afac, bw, elv, umax, sbw
! In the following pattern definitions, "u" refers to the angle for which
! the antenna pattern is sought, and "uO" refers to the elevation angle.
! IPTRN = 0 gives Omnidirectional antenna pattern factor : f(u) = 1
patfac = 1.
if( iptrn .gt. 1 ) then
u = asin( sang
udif = u - ely
end if
! IPTRN = 1 gives Gaussian antenna pattern based on
! f(p-pO) = exp(-w**2 * ( p-p )**2 ) / 4, where p = sin(u) and
! pO = sin(uO)
if( iptrn .eq. 1 ) then
pr = sang - pelev

patfac = exp(-pr * pr * afac)
! IPTRN = 2 gives sin(x)/x pattern based on
! f(u-uO) = sin(x) / x where x = afac * sin(u-u0) for Iu-uOI <= umax
! f(u-uO) = .03 for Iu-uOI > umax
! IPTRN = 4 gives height-finder pattern which is a special case of sin(x)/x
elseif(( iptrn .eq. 2 ) .or. ( iptrn .eq. 4 )) then
if( iptrn .eq. 4 ) then
dirang = abs( sang )
if( dirang .gt. elv ) udif = u - dirang
end if
if( abs(udif) .le. l.e-6 ) then
patfac = 1.
elseif( abs( udif ) .gt. umax ) then
patfac = .03
else
arg = afac * sin( udif )
patfac = aminl( 1., amaxl( .03, sin( arg ) / arg ))
end if
! IPTRN = 3 gives csc-sq pattern based on
! f(u) = 1 for u-uO <= bw
! f(u) = sin(bw) / sin(u-uO) for u-u0 > bw
! f(u) = maximum of .03 or [l+(u-uO)/bw] for u-uO < 0
elseif( iptrn .eq. 3 ) then
if( udif .gt. bw ) then
patfac = sbw / sin( udif )
elseif( udif .1t. 0 ) then
patfac = aminl( 1., amaxl( .03, (1. + udif/bw) ) )
end if
end if
end


! *************SUBROUTINE REFINIT***************
! Module Name: REFINIT
! Module Security Classification: UNCLASSIFIED
! Purpose: Initializes refractivity arrays used for subsequent PE
! calculations.
! Version Number: 1.5
! INPUTS:
! Argument List: ELERR12, RF structure, VRMAX,
! Common: NONE
! OUTPUTS:
! Argument List: IERROR
! Common: HTDUM(), IS, LVLEP, REFDUM(), RV2
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT
! Routines called: REMDUP
! GLOSSARY: For common variables refer to main glossary
! Input Variables:
! ELERR12 = Element of user-provided error flag structure EF that will
! trap on certain errors if set to .TRUE. Refer to user's
! manual.
! RF = Refractivity structure for external environmental data elements.
! RF.LVLEP = Number of levels in refractivity profile (for range
! dependent case all profiles must have same number of
! levels).
! RF.REFMSL(,) = 2-dimensional array containing refractivity with
! respect to mean sea level of each profile. Array
! format must be REFMSL(I,J) = M-unit value at Ith
! level of Jth profile. J = 1 for range-independent
! cases.
! RF.HMSL(,) = 2-dimensional array containing heights in meters with
! respect to mean sea level of each profile. Array format
! must be HMSL(I,J) = height of Ith level of Jth profile.
! J = 1 for range-independent cases.
! RF.RNGPROF() = Ranges of each profile in meters, i.e., RNGPROF(I) =
! range of Ith profile. RNGPROF(1) should always be
! equal to 0.
! RF.NPROF = Number of profiles. Equals 1 for range-independent cases.
! VRMAX = Maximum range in meters.
! Output Variables:
! IERROR = -12 -> Range of last refractivity profile entered (for range
! dependent case) is less than VRNAX. (This is returned from
! subroutine REFINIT). Will only return this error if error flag
! ELERR12 is set to .TRUE.).
! Local Variables:
! GRAD = Gradient of current refractivity level.
! HDIF = Height difference between last two differing height levels in
! each refractivity profile.
! HLARGE = This is the maximum height at which the refractivity profile
! is extrapolated.
! LVLMI = Last user-specified level in refractivity profile

! LVLM2 = Second-to-last user-specified level in refractivity profile,
! (i.e, LVLMl-1).