pe.f90

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

! NP = Final number of refractivity profiles
! RLARGE = This is the maximum range at which the refractivity profile
! is "extrapolated". For range-dependent case, the last entered
! profile is then forced to be homogeneous to a range of RLARGE.
subroutine refinit( elerrl2, vrmax, rf, IERROR)
include 'tpem.inc'
common / parinit / rv2, refdum(mxlvls), htdum(mxlvls), profint(O:maxpts), ht(O:maxpts), is, lvlep
record / refractivity / rf
logical elerrl2
data hlarge/ l.e6 /
data rlarge / l.elO /
ierror = 0
! Test to see if last profile entered ( for range dependent case ) meets or
! exceeds VRMAX, otherwise, return error (unless error trapping is turned off
! - EF.LERR12 = .FALSE.).
if( rf.nprof .gt. 1 ) then
if(( rf.rngprof(rf.nprof) .1t. vrmax ) .and. ( elerrl2 )) then
ierror = -12
return
end if
end if
! Add extra level to tabulated profiles with extrapolated gradient. Test on
! HDIF greater than 0 for profiles that contain multiple height/M-unit values
! that are equal.
rf.lvlep = rf.lvlep + 1
do i = l,rf.nprof
hdif = 0.
lvlml = rf.lvlep
lvlm2 = rf.lvlep
do while( hdif .le. l.e-6
lvlml = lvlml - 1
lvlm2 = lvlml - 1
hdif = rf.hmsl(lvlml,i) - rf.hmsl(lvlm2,i)
end do
grad = (rf.refmsl(lvlml,i)-rf.refmsl(lvlm2,i)) / hdif
! If last gradient in refractivity profile is negative then return error.
if( grad .1t. 0 ) then
ierror = -14
return
end if
rf.hmsl(rf.lvlep, i) = hlarge
rf.refmsl( rf.lvlep, i ) = (hlarge-rf.hmsl(lvlml,i)) * grad + rf.refmsl( lvlml, i
end do
is = 1
rv2=rf.rngprof(is)
do i = 1, rf.lvlep
refdum(i) = rf.refmsl( i, is )

htduxn(i) = rf.hmsl( i, is
end do
np = rf.nprof + 1.
rf.rngprof(np) = riarge
do i = 1, rf.lvlep
npml = np - 1
rf.hmsl( i, np )=rf.hmsl( i, npml
rf.refmsl( i, np )=rf.refmsl( i, npml
end do
iviep = rf.lvlep
call remdup
end


! * SUBROUTINE TRACEA *
! Module Name: TRACEA
! Module Security Classification: UNCLASSIFIED
! Purpose: This routine performs a ray trace to determine the minimum angle
! required (based on the reflected ray) to obtain a PE solution for
! all heights up to ZLIM and all ranges beyond RLIM. THETAMAX
! is then determined from this angle. This is done only for smooth
! surface and automati! angle calculation. For terrian cases,
! THETAMAX has already been set to the larger of the critical angle
! (if a duct exists), the angle that clears the highest terrain
! peak, and the tangent angle determined from HMAX and RMAX.
! If PRANG is not equal to 0, then the user has overriden the
! default calculation and THETAMAX is then determined based on
! PRANG. However a ray trace must still be done in order to
! determine the initial launch angle such that the local angle of
! the ray remains less than PRANG. The initial launch angle is
! used in subroutine TRACEH.
! Version Number: 1.5
! INPUTS:
! Argument List: ACRIT, PRANG, TR structure
! Common: ANTREF, DMDH(), FTER, HTDUMO, JLS, LVLEP, RLIM,
! SLP(, THETAMAX, ZLIM
! OUTPUTS:
! Argument List: NONE
! Common: RPE, THETALAUNCH, THETAMAX
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT
! Routines called: NONE
! GLOSSARY: For common variables refer to main glossary
! Input Variables:
! ACRIT = Critical angle, i.e., angle above which no rays are trapped
! for ducting environment.
! PRANG = User-specified propagation angle in radians.
! TR = Terrain structure for external terrain data elements.
! TR.TERX() = range points of terrain profile in meters.
! TR.TERY() = height points of terrain profile in meters.
! Output Variables:
! For common variables refer to main glossary
! Local Variables:
! AO = Angle in radians of ray before trace step.
! Al = Angle in radians of ray after trace step.
! AMXCUR = Maximum of local angle along ray.
! AMXCURL = Last (or previous) AMXCUR.
! AS = Starting launch angle in radians.
! ASL = Last (or previous) starting launch angle.
! DEG15 = 15 degrees (in radians) - used as a maximum limit for
! THETAMAX and THETALAUNCH.
! GRAD = Gradient of current refractivity layer.
! HO = Height of ray in meters of ray before trace step.

! Hl = Height of ray in meters of ray after trace step.
! IDN = Integer with value of + or - 1. Used to increment or decrement
! initial launch angle.
! ISET = Flag to test whether or not to stop main loop.
! JL = Index of current refractivity layer ray tracing through.
! LOOP = Flag to test whether or not to stop nested loop of ray tracing
! individual ray.
! RAD = Radicand for ray trace formula
! RO = Range of ray in meters before trace step.
! R1 = Range of ray in meters after trace step.
! RREF = Range at which traced ray is reflected.
subroutine tracea( tr, prang, acrit )
include 'tpem.inc'
common / miscvar / fnorm, cnst, delp, thetamax, plcnst, qi, antref, rpe, hlim(mxrout), slp(mxter), fter, hmref, htlim
common / trvar / dmdh(mxlvls), zlim, jls, thetalaunch, rlim
common / parinit / rv2, refdum(mxlvls), htdum(mxlvls), profint(O:maxpts), ht(O:maxpts), is, lvlep
record / terrain / tr
complex qi
logical fter, loop
data degl5 / .2617994 / !15 degrees in radians
! All heights and ranges are in meters, gradients are in M-unit/meter * l.e-6
! and angles are in radians.
! Define in line ray trace functions:
radal( a, b ) = a**2 + 2. * grad * b !a=aO, b=hl-hO
rp( a, b ) = a + b / grad !a=rO, b=al-aO
ap( a, b ) = a + b * grad !a=aO, b=rl-rO
hp( a, b, c ) = a + ( b**2 - c**2 ) / 2. / grad !a=hO, b=al, c=aO
as = -thetamax
idn = -1
if( fter ) then
as = thetamax
if( prang .le. l.e-6 ) idn = 1
! For initial shallow slope and non-zero user-defined maximum propagation
! angle, determine the range at which ray is reflected (RREF). If this is
! less than the range of the 2nd terrain point, then treat as if this were a
! smooth surface problem. I.e., determine THETAMAX based on reflected ray,
! not direct ray.
if(( prang .gt. l.e-6 ) .and. ( slp(l) .le. l.e-6 )) then
rref = (antref - tr.tery(l)) / tan( thetamax )
if( rref .1t. tr.terx(2) ) then
idn = -1
as = -thetamax
end if
end if
end if
asl = 0.
amxcurl = 0.
iset = 0
do while( iset .eq. 0 )

! Decrease or increase angle by 1 mrad, depending on value of IDN.
! Initialize ray trace variables.
as = as + idn*l.e-3
hO = antref
ro = 0.
rpe = 0.
aO = as
jl = jls
amxcur = 0.
loop = .true.
! Perform ray trace until ray has reached ZLIM and/or RLIM where
! ZLIM = maximum of HMAX-HMINTER or ANTREF.
! RLIM = .9 * RMAX
do while( loop
grad = dmdh(jl)
if( a0 .it. 0. ) hl = htdum(jl)
if( aO .gt. 0. ) hl = htdum(jl + 1)
if( aO .eq. 0. ) then
if( grad .1t. 0. ) hl = htdum(jl)
if( grad .gt. 0. ) hl = htdum(jl+l)
end if
if( hl .gt. zlim ) hl = zlim
rad = radal( aO, hl-hO
if( rad .gt. 0 ) then
al = sign( 1., aO ) * sqrt( rad )
else
al = 0.
hl = hp( hO, al, aO
end if
rl = rp( rO, al-aO
if(( al .le. 0. ) .and. ( hl .le. htdum(jl) )) then
hl = htdum(jl)
rad = radal( a0, hl-hO )
al = -sqrt( rad )
rl = rp( rO, al-aO )
ji = jl - 1
if( jl .eq. 0 ) jl = 1
elseif(( al .ge. 0. ).and.( hl .ge. htdum(jl+l) ))then
hl = htdum(jl+l)
rad = radal( aO, hi-hO
al = sqrt( rad)
rl = rp( rO, al-aO
ji = ji + 1
if( jl .gt. ivlep ) jl = ivlep
end if
if( hl .gt. zlim ) then
hl = zlim
rad = radal( aO, hi-hO
al = sqrt( rad )
rl = rp( rO, al-aO )
end if
hO = hl
rO = rl
aO = al
! Set RPE to range at which reflected ray hits the ground.
if( hO .le. 1.e-4 ) then
aO = -a0

rpe = rO
end if
! If AO greater than 90 degrees then exit loop.
if( aO .ge. 1.57079 ) exit
amxcur = amaxl( amxcur, aO )
if(( hO .ge. zlim ) .and. ( aO .gt. 0.)) loop = .false.
if( rO .gt. rlim ) loop = .false.
end do
! Test to see if the current ray traced from launch angle AS meets criteria.
! If ray traced does not reach ZLIM AND is not within RLIM then the initial
! launch angle AS is increased by 1 mrad and ray trace is
! repeated. This is done only for smooth surface.
if(( rO .le. rlim ) .and. ( rpe .gt. 0. )) iset = 1
! If criteria is met then (if user specified a problem angle) make sure the
! local maximum angle is just within PRANG. If not then repeat ray trace
! until this occurs.
if( fter )then
if( prang .gt. l.e-6 ) then
iset = 1
if( amxcur .gt. thetamax ) iset = 0
if( as .le. acrit+l.e-3 ) iset = 1 !Don't let launch angle
!be less than critical
!angle.
else
if(( rO .le. rlim ) .and. ( hO .ge. zlim )) iset = 1
if( iset .eq. 1 ) thetamax = amaxl( abs(as), amxcur )
end if
else
if(( prang .gt. l.e-6 ) .and. (iset .eq. 1)) then
a = amaxl( abs(as), amxcur
if( a .1t. prang ) then
iset = 0
elseif( asl .ne. 0. ) then
as = asl
amxcur = amxcurl
end if
end if
end if
! Just as a safeguard - set absolute maximum of launch angle to 15 degrees.
if( as .le. -degl5 ) then
iset = 1
as = -degl5
amxcur = degl5
end if
asl = as
amxcurl = amxcur
end do
if( .not. fter ) thetamax = amaxl( abs(as), amxcur )
thetalaunch = abs(as)
end


! *************SUBROUTINE DIEINIT**************
! Module Name: DIEINIT
! Module Security Classification: UNCLASSIFIED
! Purpose: This routine calculates Conductivity and Permittivity
! as a function of frequency in MHz. All equations and coefc
! ficients were obtained by using a SUMMASKETCH digitizer with
! the CCIR volume 5 curves on page 74. The digitized data was
! then used with the TABLECURVE software to obtain the best fit
! equations and coefficients used in this subroutine. In some
! cases two sets of equations were required to obtain a decent
! fit across the 100 MHz - 100GHz range. These curves fit the
! digitized data to within 5%.
! Version Number: 1.5
! INPUTS:
! Argument List: SV structure, TR structure
! Common: NONE
! OUTPUTS: TR.DIELEC(,)
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT
! Routines called: NONE
! GLOSSARY:
! Input Variables:
! SV = System structure for external system data elements.
! SV.FREQ = Frequency in MHz.
! TR = Terrain structure for external terrain data elements.
! TR.IGR = Number of different ground types specified.
! TR.IGRND() = Type of ground composition for given terrain profile -
! can vary with range. Different ground types are:
! 0 = sea water, 1 = fresh water, 2 = wet ground,
! 3 = medium dry ground, 4 = very dry ground,
! 5 = user defined (in which case, values of relative
! permittivity and conductivity must be given).
! TR.RGRND() = Ranges at which the ground types apply.
! Output Variables:
! TR.DIELEC(,) = 2-dimensional array containing the relative
! permittivity and conductivity; DIELEC(l,i) and
! DIELEC(2,i), respectively. Only needs to be specified
! if using IGRND(i) = 5, otherwise, TPEM will
! calculate based on frequency and ground types 0-4.
! Local Variables:
! EPSILON = Relative permittivity.
! SIGMA = Conductivity.
subroutine dieinit( sv, tr )
include 'tpem.inc'
record / terrain / tr
record / systemvar / sv
dimension a(14), b(14), c(14), d(14), e(14), f(14)

data (a(i),i=1,14) / 1.4114535e-2, 3.8586749, 79.027635,
+ -0.65750351, 201.97103, 857.94335,
+ 915.31026, 0.8756665, 5.5990969e-3,
+ 215.87521, .17381269, 2.4625032e-2,
+ -4.9560275e-2, 2.2953743e-4 /
data (b(i),i=1,14) / -5.2122497e-8, -2.1179295e-5, -2.2083308e-5,
+ 5.5620223e-5, -2.5539582e-3, -8.9983662e-5,
+ -9.4530022e-6, 4.7236085e-5, 8.7798277e-5,
+ -7.6649237e-5, 1.2655183e-4, 1.8254018e-4,
+ 2.9876572e-5, -8.1212741e-7 /
data (c(i),i=1,14) / 5.8547829e-11, 9.1253873e-4, -3.5486605e-4,
+ 6.6113198e-4, 1.2197967e-2, 5.5275278e-2,
+ -4.0348211e-3, 2.6051966e-8, 6.2451017e-8,
+ -2.6151055e-3, -1.6790756e-9, -2.664754e-8,
+ -3.0561848e-10, 1.8045461e-9 /
data (d(i),i=1,14) / -7.6717423e-16, 6.5727504e-10, 2.7067836e-9,
+ 3.0140816e-10, 3.7853169e-5, 8.8247139e-8,
+ 4.892281e-8, -9.235936e-13, -7.1317207e-12,
+ 1.2565999e-8, 1.1037608e-14, 7.6508732e-12,
+ 1.1131828e-15, -1.960677e-12 /
data (e(i),i=1,14) / 2.9856318e-21, 1.5309921e-8, 8.210184e-9,
+ 1.4876952e-9, -1.728776e-6, 0.0,
+ 7.4342897e-7, 1.4560078e-17, 4.2515914e-16,
+ 1.9484482e-7, -2.9223433e-20, -7.4193268e-16,
+ 0.0, 1.2569594e-15 /
data (f(i),i=1,14) / 0., -1.9647664e-15, -1.0007669e-14, 0., 0.,
+ 0., 0., -1.1129348e-22, -1.240806e-20, 0.,
+ 0., 0., 0., -4.46811e-19 /
fl = sv.freq
f2 = fl * fl
f3 = fl * f2
f4 = fl * f3
f5 = fl * f4
f6 = fl * f5
f7 = fl * f6
f8 = fl * f7
f9 = fl * f8
do i = 1, tr.igr
select case ( tr.igrnd(i)
! Permittivity and conductivity for salt water.
case( 0 )
epsilon = 70.
sigma = 5.
m= 1
ml = m + 1
if( fl .gt. 2253.5895 ) epsilon = 1. / (a(m) +
+ b(m)*fl + c(m)*f2 + d(m)*f3 + e(m)*f4
if( fl .gt. 1106.207 ) then
sigma = a(ml) + c(ml)*fl + e(ml)*f2
sigma = sigma / ( 1.+ b(ml)*fl + d(ml)*f2 + f(ml)*f3
end if
! Permittivity and conductivity for fresh water.
case( 1 )
epsilon = 80.0
m=3
ml = m + 1
IF( fl .gt. 6165.776 ) THEN
epsilon = a(m) + c(m)*fl + e(m)*f2
epsilon = epsilon/(l. + b(m)*fl + d(m)*f2 + f(m)*f3 )

end if
IF( fl .gt. 5776.157) THEN
k =2
else
ml = ml + 1
k = -1
end if
sigma = a(ml) + c(ml)*fl + e(ml)*f2
sigma = (sigma / (1. + b(ml)*fl + d(ml)*f2))**k
! Permittivity and conductivity for wet ground.
case( 2 )
epsilon = 30.0
m= 6
IF( fl .ge. 4228.11 ) m = 7
if( fl .gt. 1312.054 ) then
epsilon = a(m) + c(m)*fl + e(m)*f2
epsilon = SQRT( epsilon / (1. + b(m)*fl + d(m)*f2) )
end if
IF( fl .gt. 15454.4) then
ml = 8
g = 3.3253339e-28
else
ml = 9
g = 1.3854354e-25
end if
sigma = a(ml) + b(ml)*fl + c(ml)*f2 + d(ml)*f3 + e(ml)*f4
sigma = sigma + f(ml)*f5 + g*f6
! Permittivity and conductivity for medium dry ground.
case( 3 )
epsilon = 15.0
IF( fl .gt. 4841.945) THEN
m = 10
epsilon = a(m) + c(m)*fl + e(m)*f2
epsilon = SQRT( epsilon / (1. + b(m)*fl + d(m)*f2) )
end if
ml = 12
IF( fl .gt. 4946.751) ml = 11
sigma = (a(ml) + b(ml)*fl + c(ml)*f2 + d(ml)*f3 + e(ml)*f4)**2
! Permittivity and conductivity for very dry ground.
case( 4 )
epsilon = 3.0
IF( fl .1t. 590.8924 ) then
sigma = 1.0e-4
else
IF( fl .gt. 7131.933) THEN
ml = 13
sigma = (a(ml) + b(ml)*fl + c(ml)*f2 + d(ml)*f3)**2
else