pe.f90

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

ml = 14
g = 9.4623158e-23
h = -1.1787443e-26
s = 7.9254217e-31
t = -2.2088286e-35
sigma = a(ml) + b(ml)*fl + c(ml)*f2 + d(ml)*f3
sigma = sigma + e(ml)*f4 + f(ml)*f5 + g*f6
sigma = sigma + h*f7 + s*f8 + t*f9
end if
end if
! User-defined

case( 5 )
epsilon = tr.dielec(l,i)
sigma = tr.dielec(2,i)
case default
! Do nothing
end select
tr.dielec(l,i) = epsilon
tr.dielec(2,i) = sigma
end do
! Set dielectri! constants equal to last provided ground constants at le7 km.
tr.igr = tr.igr + 1
tr.rgrnd(tr.igr) = l.elO
tr.dielec(l,tr.igr) = epsilon
tr.dielec(2,tr.igr) = sigma
end

! ******************* SUBROUTINE GETFFTSZ ***************************
! Module Name: GETFFTSZ
! Module Security Classification: UNCLASSIFIED
! Purpose: Determines the FFT size needed for a given problem.
! Version Number: 1.5
! INPUTS:
! Argument List: ZLIM
! Common: FTER, THETAMAX, WL
! Parameter: MXNFFT
! OUTPUTS:
! Argument List: ZLIM
! Common: DELZ, LN, N, ZMAX
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT
! Routines called: NONE
! GLOSSARY: For common variables refer to main glossary. For parameters
! refer to TPEM.INC and FFTSIZ.INC
! Input Variables:
! ZLIM = Maximum height region where PE solution is valid = .75 * ZMAX.
! Output Variables:
! ZLIM = Calculates a new height limit equal to .75*ZMAX only if the
! maximum transform size needed is too large to do specified
! problem. Fixes transform size to maximum and adjusts ZMAX and
! ZLIM accordingly.
! Local Variables:
! STHETAMAX = Sine of THETAMAX
subroutine getfftsz( ZLIM )
include 'tpem. inc'
common / miscvar / fnorm, cnst, delp, thetamax, plcnst, qi,
+ antref, rpe, hlim(mxrout), slp(mxter), fter, hmref, htlim
common / pevar / wl, fko, delz, n, ln, zmax, n34, con, dz2, nml
logical fter
complex qi
sthetamax = sin( thetamax
delz= wl * .5 / sthetamax
! Set lower FFT limit to 2**9 for smooth surface cases, if terrain case then
! set lower FFT limit to 2**10.
ln = 9
if( fter) ln=10
N=2**LN

zmax=delz*float(n)
! Determine transform size needed to perform calculations to a height of ZLIM,
! up to the maximum FFT size allowed.
do while( .75*zmax .1t. zlim
in = in + 1
if( in .gt. mxnfft ) exit
n = 2**ln
zmax = delz * float(n)
end do
! If the transform size needed is too large then set ZMAX and ZLIM
! accordingly.
if( in .gt. mxnfft ) then
in = mxnfft
n = 2**ln
zmax = delz * float(n)
zlim = .75 * zmax
end if
end

! **************** SUBROUTINE XYINIT*************
! Module Name: XYINIT
! Module Security Classification: UNCLASSIFIED
! Purpose: Determines the initial PE starter field.
! Version Number: 1.5
! INPUTS:
! Argument List: SV structure, TR structure
! Common: DELP, FILT(, FKO, N, N34, RNG2, WL, ZMAX
! OUTPUTS:
! Argument List: NONE
! Common: RNG2, U()
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT
! Routines called: ANTPAT, GETALN
! GLOSSARY: For common variables see main glossary
! Input Variables:
! SV = System structure for external system data elements.
! SV.ANTHT = transmitting antenna height above local ground in meters.
! SV.IPAT = integer value indicating type of antenna pattern desired.
! IPAT = 0 ->omni
! IPAT = 1 -> gaussian
! IPAT = 2 ->sinc x
! IPAT = 3 -> csc**2 x
! IPAT = 4 -> generic height-finder
! SV.POLAR = 1-character string indicating polarization. H-horizontal,
! V-vertical.
! TR = Terrain structure for external terrain data elements.
! Output Variables:
! For common variables refer to main glossary
! Local Variables:
! ANTKO = Exponential term in calculation of DTERM and RTERM
! ATTN = Attenuation factor for filtering
! CRAD = Square root term in reflection coefficient calculation
! CTHETA = Cosine of angle PK
! DTERM = Field due to real point source at height SV.ANTHT, i.e.,
! DTERM=exp{-i*sv.antht*fko*sin(theta)}.
! DTHETA = Bin width in angle- (or p-) space
! FACD = Antenna pattern factor for direct ray angle.
! FACR = Antenna pattern factor for reflected ray angle.
! PK = Angle value for bin I, i.e., I*DTHETA
! REFCOEF = Complex reflection coefficient.
! RTERM = Field due to image point source at height -SV.ANTHT, i.e.,
! RTERM=exp{i*sv.antht*fko*sin(theta)).
! SGAIN = Normalization factor.
! SRAD = Sine term in reflection coefficient calculation
! STHETA = Sine of angle PK
! ZPK = Phase term for DTERM and RTERM

SUBROUTINE xyinit( sv, tr )

include 'tpem.inc'
common /arrays /u(O:maxpts), filt(O:maxn4), frsp(O:maxpts), envpr(O:maxpts), ulst(O:maxpts)
common /miscvar /fnorm, cnst, deip, thetamax, plcnst, qi,
+ antref, rpe, hlim(mxrout), slp(mxter), fter, hmref, htlim.
common /pevar / wi, fko, delz, n, in, zmax, n34, con, dz2, nml
common /impedance /alphav, rav(O:maxpts), rng, rng2, ci, c2, rk, din, c2m, ig, root
record /systemvar /sv
record /terrain / tr
logical fter
complex u, frsp, envpr, uist, qi, root
complex alphav, rav, rng, rng2, ci, c2, rk, din, c2m
complex refcoef, rterm, dterm, crad, srad
! Reflection coefficient is defaulted to -1 for horizontal polarization.
refcoef = cmplx( -1., 0.)
if( sv.polar .eq. IV' ) call getaln( tr )
sgain= sqrt( wl ) / zinax
dtheta =delp / fko
antko =fko * sv.antht
DO I=0,N
pk =float(i) * dtheta
zpk =pk * antko
! Get antenna pattern factors for the direct and reflected rays.
call antpat( sv.ipat, pk, FACD )
call antpat( sv.ipat, -pk, FACR )
! If vertical polarization, then determine reflection coefficient.
if( sv.polar .eq. I'V) then
ctheta = cos( pk )
stheta =sin( pk )
crad = csqrt( rng2 -ctheta**2 )
srad = rng2 * stheta
refcoef = (srad - crad) / (srad + crad)
end if
rterm = cmplx( cos( zpk ), sin( zpk ))
dterm = conjg( rterm)
u(i) = sgain *Cfacd * dterin + refcoef *facr *rterm
end do
! Filter upper 1/4 of the field.
do i = n34, n
attn = filt(i-n34)
u(i) = attn*u(i)
end do
END


! ********************** SUBROUTINE FFT *****************************
! Module Name: FFT
! Module Security Classification: UNCLASSIFIED
! Purpose: Performs fast Fourier sine transform on complex PE field array U.
! Version Number: 1.5
! INPUTS:
! Argument List: U()
! Common: LN, N
! Parameter: MAXPTS
! OUTPUTS:
! Argument List: U()
! Files Included: FFTSIZ.INC, TPEM.INC
! Calling Routines: PEINIT, PESTEP
! Routines Called: SINFFT
! GLOSSARY: For common variables refer to main glossary. For parameters,
! refer to TPEM.INC
! Input Variables:
! U() = Complex field to be transformed.
! Output Variables:
! U0) = Transform of complex field.
! Local Variables:
! X0) = Real part of field.
! Y() = Imaginary part of field.
subroutine fft( u )
include 'tpem.inc'
common / pevar / wl, fko, delz, n, ln, zmax, n34, con, dz2, nml
complex u(O:*)
dimension x(O:maxpts), y(O:maxpts)
do i = 0, n
x(i) = real( u(i) )
y(i) = imag( u(i) )
end do
call sinfft( ln, X )
call sinfft( ln, Y )
do i = 0, n
u(i) = cmplx( x(i), y(i) )
end do
end


SUBROUTINE SINFFT ( N, X )
! *********************************************************************************
!  PURPOSE: SINFFT replaces the real array X() by its finite discrete sine transform
!  METHOD : 
!           The algorithm is based on a mixed radix (8-4-2) real vector
!           fast Fourier synthesis routine published by Bergland:
! 
!  G.D. Bergland, 'A Radix-eight Fast Fourier Transform 
!  Subroutine for Real-valued Series,' IEEE Transactions on
!  Audio and Electro-acoustics', vol. AU-17, pp. 138-144, 1969
!  and sine and cosine transform algorithms for real series
!  published by Cooley, Lewis, and Welch: 
!  (J.W. COOLEY, P.A.W. LEWIS AND P.D. WELSH, 'The Fast Fourier
!  Transform Algorithm: Programming Considerations in the
!  Calculation of Sine, Cosine and Laplace Transforms',
!  J. SOUND VIB., vol. 12, pp. 315-337, 1970 ). *
!  
!  ARGUMENTS: *
!    INPUT:
!    N ...... transform size ( = 2**N ) *
!    X0 .... data array dimensioned 2**N in calling program *
!  -- OUTPUT -- *
!    X0 .... sine transform *
!    TABLES: array required size *
!        B 2**N *
!        JINDX 2**(N-l) *
!        COSTBL 2**(N-4) *
!        SINTBL 2**(N-4)
! 
!  Sub-programs called: - *
!
!  R8SYN ..... (radix 8 synthesis) *
! *************************************************************
!
include 'fftsiz.inc'
INTEGER*4 N

DIMENSION X(0:*)
INTEGER*4 NMAX2, NMAX16, NP, NPD2, NPD4

PARAMETER ( NMAX2 = MAXPTS/2
PARAMETER ( NMAX16 = MAXPTS/16
DIMENSION B(MAXPTS), JINDX (NMAX2)
DIMENSION SINES (MAXPTS)
DIMENSION COSTBL (NMAXI6), SINTBL (NMAXI6)


SAVE B, COSTBL, JINDX, SINES, SINTBL
SAVE NSAVE, N4, N8, NP, NPD2, NPD4, NPD16, NPMI

DOUBLE PRECISION ARG, DT, PI
DATA NSAVE / 0 /
DATA PI / 3.1415 92653 58979 32D0 /
! 
! - -----------------------------------------------------------------------
IF ( N .NE. NSAVE ) THEN
!  compute constants and construct tables
NSAVE = N
N8 = NSAVE / 3
N4 = NSAVE - 3 * N8 - 1
NP = 2**N
NPD2 = NP/ 2
NPD4 = NP / 4
NPD16 = NP / 16
NPMI = NP -1
!  build reciprical sine table
DT = PI /FLOAT ( NP)
DO 10 J = 1, NPMI
ARG = DT * J
SINES ( J ) = (0.5D0 / SIN ( ARG ))
10 CONTINUE
! construct bit reversed subscript table
J1 = 0
DO 30 J = 1, NPD2 - I
J2 = NPD2
20 CONTINUE
IF ( IAND ( Jl, J2 ) .NE. 0 ) THEN
Jl =IABS ( Jl -J2)
J2 = J2 / 2
GO TO 20
ENDIF
Jl = Jl + J2
JINDX ( J ) = Jl
30 CONTINUE
! 
!  form the trig tables for the radix-8 passes;
!  tables are stored in bit reversed order.
Jl = 0
DO 50 J = 1, NPD16 - 1
J2 = NPD16
40 CONTINUE
IF ( IAND ( Jl, J2 ) .NE. 0 ) THEN
Jl = IABS ( Jl - J2 )
J2 = J2 / 2
GO TO 40
ENDIF
Jl = Jl + J2
ARG = DT * FLOAT (Jl)
COSTBL ( J ) = COS ( ARG )
SINTBL ( J ) = -SIN ( ARG )
50 CONTINUE

ENDIF
! 
!  *** form the input Fourier coefficients ***
! 
!  sine transform
B ( 1 ) = -2. * X ( 1 )
B ( 2 ) = 2. * X ( NPMI )
J1 = 0
DO 110 J = 3, NPM1, 2

Jil = Ji + 1
J2 = JINDX ( Ji )
B (J ) =X (J2 - 1 )-X ( J2 + 1)
B ( J + 1 ) = X ( NP-J2 )
110 CONTINUE

!  * *
!  * Begin Fast Fourier Synthesis *
! 
! 
IF ( N8 .NE. 0 ) THEN
!  radix-8 iterations
INTT = 1
NT = NPDI6
DO 130 J = 1, N8
Jl = 1 + INTT
J2 = Jl + INTT
J3 = J2 + INTT
J4 = J3 + INTT
J5 = J4 + INTT
J6 = J5 + INTT
J7 = J6 + INTT
! **
CALL R8SYN (INTT, NT, COSTBL, SINTBL, B(1), B(Ji), B(J2),B(J3), B(J4), B(J5), B(J6), B(J7) )
! **
NT = NT / 8
INTT = 8 * INTT
130 CONTINUE
ENDIF
! radix-4 iteration
IF ( N4 .GT. 0 ) THEN
Jl = NPD4
J2 = 2*NPD4
J3 = 3*NPD4
DO 140 J = 1, NPD4
TO = B(J) + B(J + Jl)
Ti = B(J) - B(J + Jl)
T2 = 2. * B(J + J2)
T3 = 2. * B(J + J3)
B(J) = TO + T2
B(J + J2) = TO - T2
B(J + Jl) = TI + T3
B(J + J3) = TI - T3
140 CONTINUE

ELSE IF ( N4 .EQ. 0 ) THEN
! radix-2 iteration
K = NPD2
DO 150 J = 1, NPD2
K =K + 1
T =B(J) + B (K)
B(K) = B(J) - B (K)
B(J) = T
150 CONTINUE
ENDIF

! * *
! Form Transform *
! * *
! sine transform

Jl = NP
DO 160 J = 1, NPMl
X(J) = .25*(( B(J+I) + B(J1)) * SINES(J) - B(J+I) + B(Jl))
J1 = J1 - 1
160 CONTINUE
RETURN
END
! 
! 
SUBROUTINE R8SYN ( INTT, NT, COSTBL, SINTBL, BO, Bl, B2, B3, B4, B5, B6, B7 )
! 
! 
!  PURPOSE: Radix-8 synthesis subroutine used by mixed radix driver.
! 
! 
! 
DIMENSION COSTBL(*), SINTBL(*)
DIMENSION BO(*), BI(*), B2(*), B3(*), B4(*), B5(*), B6(*), B7(*)
! 
! 
!  /// Local variables //
! 
! 
DOUBLE PRECISION Cl, C2, C3, C4, C5, C6, C7
DOUBLE PRECISION SI, S2, S3, S4, S5, $6, S7
DOUBLE PRECISION CPI4, CPI8, R2, SP18
SAVE CPI4, CPI8, R2, SPI8
! 
DATA R2 / 1.41421 35623 7310D+0 /,