Coverage Report: RefCoef.f90

Generated from GCOV analysis of Fortran source code

47.3%
Lines Executed
110 total lines
54.8%
Branches Executed
208 total branches
100.0%
Calls Executed
50 total calls
0
-
Source:RefCoef.f90
0
-
Graph:RefCoef.gcno
0
-
Data:RefCoef.gcda
0
-
Runs:73
1
-
!! Provides data and functions for working with reflection coefficients
2
-
3
-
MODULE RefCoef
4
-
!! Provides reflection coefficient data
5
-
6
-
USE FatalError
7
-
USE monotonicMod
8
-
9
-
IMPLICIT NONE
10
-
PUBLIC
11
-
SAVE
12
-
13
-
INTEGER, PARAMETER :: BRCFile = 31, TRCFile = 32, IRCFile = 12
14
-
INTEGER :: NBotPts, NTopPts
15
-
INTEGER :: NkTab
16
-
INTEGER, ALLOCATABLE :: iTab( : )
17
-
REAL (KIND=8), ALLOCATABLE :: xTab( : )
18
-
COMPLEX (KIND=8), ALLOCATABLE :: fTab( : ), gTab( : )
19
-
20
-
TYPE ReflectionCoef
21
-
REAL(KIND=8) :: theta, R, phi
22
-
END TYPE ReflectionCoef
23
-
24
-
TYPE(ReflectionCoef), ALLOCATABLE :: RBot( : ), RTop( : )
25
-
26
-
CONTAINS
27
-
28
71*
SUBROUTINE ReadReflectionCoefficient( FileRoot, BotRC, TopRC, PRTFile )
29
-
30
-
! Optionally read in reflection coefficient for Top or Bottom boundary
31
-
32
-
USE MathConstants
33
-
IMPLICIT NONE
34
-
INTEGER, INTENT( IN ) :: PRTFile ! unit number for print file
35
-
CHARACTER (LEN=1 ), INTENT( IN ) :: BotRC, TopRC ! flag set to 'F' if refl. coef. is to be read from a File
36
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
37
-
INTEGER :: itheta, ik, IOStat, iAllocStat
38
-
REAL ( KIND = 8 ) :: freq
39
-
CHARACTER (LEN=80) :: Title2
40
-
41
71
IF ( BotRC == 'F' ) THEN
42
1
WRITE( PRTFile, * ) '__________________________________________________________________________'
43
1
WRITE( PRTFile, * )
44
1
WRITE( PRTFile, * ) 'Using tabulated bottom reflection coef.'
45
1*
OPEN( FilE = TRIM( FileRoot ) // '.brc', UNIT = BRCFile, STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
46
1
IF ( IOStat /= 0 ) THEN
47
#####
WRITE( PRTFile, * ) 'BRCFile = ', TRIM( FileRoot ) // '.brc'
48
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Unable to open Bottom Reflection Coefficient file' )
49
-
END IF
50
-
51
1
READ( BRCFile, * ) NBotPts
52
1
WRITE( PRTFile, * ) 'Number of points in bottom reflection coefficient = ', NBotPts
53
-
54
1*
IF ( ALLOCATED( RBot ) ) DEALLOCATE( RBot )
55
1*
ALLOCATE( RBot( NBotPts ), Stat = IAllocStat )
56
1
IF ( IAllocStat /= 0 ) THEN
57
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Insufficient memory for bot. refl. coef.: reduce # points' )
58
-
END IF
59
-
60
4*
READ( BRCFile, * ) ( RBot( itheta ), itheta = 1, NBotPts )
61
4*
IF ( .NOT. monotonic( RBot( : )%theta, NBotPts ) ) THEN
62
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Bottom reflection coefficients must be monotonically increasing' )
63
-
END IF
64
-
65
1
CLOSE( BRCFile )
66
4*
RBot%phi = DegRad * RBot%phi ! convert to radians
67
-
68
-
ELSE ! should allocate something anyway, since variable is passed
69
70*
IF ( ALLOCATED( RBot ) ) DEALLOCATE( RBot )
70
70*
ALLOCATE( RBot( 1 ), Stat = IAllocStat )
71
-
END IF
72
-
73
-
! Optionally read in top reflection coefficient
74
-
75
71
IF ( TopRC == 'F' ) THEN
76
1
WRITE( PRTFile, * ) '__________________________________________________________________________'
77
1
WRITE( PRTFile, * )
78
1
WRITE( PRTFile, * ) 'Using tabulated top reflection coef.'
79
1*
OPEN( FILE = TRIM( FileRoot ) // '.trc', UNIT = TRCFile, STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
80
1
IF ( IOStat /= 0 ) THEN
81
#####
WRITE( PRTFile, * ) 'TRCFile = ', TRIM( FileRoot ) // '.trc'
82
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Unable to open Top Reflection Coefficient file' )
83
-
END IF
84
-
85
1
READ( TRCFile, * ) NTopPts
86
1
WRITE( PRTFile, * ) 'Number of points in top reflection coefficient = ', NTopPts
87
-
88
1*
IF ( ALLOCATED( RTop ) ) DEALLOCATE( RTop )
89
1*
ALLOCATE( RTop( NTopPts ), Stat = IAllocStat )
90
1
IF ( iAllocStat /= 0 ) THEN
91
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Insufficient memory for top refl. coef.: reduce # points' )
92
-
END IF
93
-
94
4*
READ( TRCFile, * ) ( RTop( itheta ), itheta = 1, NTopPts )
95
4*
IF ( .NOT. monotonic( RTop( : )%theta, NTopPts ) ) THEN
96
#####
CALL ERROUT( 'ReadReflectionCoefficient', 'Top reflection coefficients must be monotonically increasing' )
97
-
END IF
98
-
99
1
CLOSE( TRCFile )
100
4*
RTop%phi = DegRad * RTop%phi ! convert to radians
101
-
ELSE ! should allocate something anyway, since variable is passed
102
70*
IF ( ALLOCATED( RTop ) ) DEALLOCATE( RTop )
103
70*
ALLOCATE( RTop( 1 ), Stat = iAllocStat )
104
-
END IF
105
-
106
-
! Optionally read in internal reflection coefficient data
107
-
108
71
IF ( BotRC == 'P' ) THEN
109
#####
WRITE( PRTFile, * ) 'Reading precalculated refl. coeff. table'
110
#####
OPEN( FILE = TRIM( FileRoot ) // '.irc', UNIT = IRCFile, STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
111
#####
IF ( IOStat /= 0 ) CALL ERROUT( 'ReadReflectionCoefficient', &
112
#####
'Unable to open Internal Reflection Coefficient file' )
113
-
114
#####
READ( IRCFile, * ) Title2, freq
115
#####
READ( IRCFile, * ) NkTab
116
#####
WRITE( PRTFile, * )
117
#####
WRITE( PRTFile, * ) 'Number of points in internal reflection coefficient = ', NkTab
118
-
119
#####
IF ( ALLOCATED( xTab ) ) DEALLOCATE( xTab, fTab, gTab, iTab )
120
#####
ALLOCATE( xTab( NkTab ), fTab( NkTab ), gTab( NkTab ), iTab( NkTab ), Stat = iAllocStat )
121
#####
IF ( iAllocStat /= 0 ) CALL ERROUT( 'ReadReflectionCoefficient', 'Too many points in reflection coefficient' )
122
-
123
#####
READ( IRCFile, FMT = "( 5G15.7, I5 )" ) ( xTab( ik ), fTab( ik ), gTab( ik ), iTab( ik ), ik = 1, NkTab )
124
#####
CLOSE( IRCFile )
125
-
END IF
126
-
127
71
END SUBROUTINE ReadReflectionCoefficient
128
-
129
143339
SUBROUTINE InterpolateReflectionCoefficient( RInt, R, NPts, PRTFile )
130
-
131
-
! Given an angle RInt%ThetaInt, returns the magnitude and
132
-
! phase of the reflection coefficient (RInt%R, RInt%phi).
133
-
134
-
! Uses linear interpolation using the two nearest abscissas
135
-
! Assumes phi has been unwrapped so that it varies smoothly.
136
-
! I tried modifying it to allow a complex angle of incidence but
137
-
! stopped when I realized I needed to fuss with a complex atan2 routine
138
-
139
-
IMPLICIT NONE
140
-
INTEGER, INTENT( IN ) :: PRTFile, NPts ! unit number of print file, # pts in refl. coef.
141
-
TYPE(ReflectionCoef), INTENT( IN ) :: R( NPts ) ! Reflection coefficient table
142
-
TYPE(ReflectionCoef), INTENT( INOUT ) :: RInt ! interpolated value of refl. coef.
143
-
INTEGER :: iLeft, iRight, iMid
144
-
REAL ( KIND=8 ) :: alpha, Thetaintr
145
-
146
143339
iLeft = 1
147
143339
iRight = NPts
148
-
149
143339
thetaIntr = REAL( RInt%Theta ) ! This should be unnecessary? probably used when I was doing complex angles
150
-
151
-
! Three cases: ThetaInt left, in, or right of tabulated interval
152
-
153
143339*
IF ( thetaIntr < R( iLeft )%theta ) THEN
154
-
!iRight = 2
155
#####
RInt%R = 0.0 ! R( iLeft )%R
156
#####
RInt%phi = 0.0 ! R( iLeft )%phi
157
#####
WRITE( PRTFile, * ) 'Warning in InterpolateReflectionCoefficient : Refl. Coef. being set to 0 outside tabulated domain'
158
#####
WRITE( PRTFile, * ) 'angle = ', thetaintr, 'lower limit = ', R( iLeft)%theta
159
-
160
143339*
ELSE IF( thetaIntr > R( iRight )%theta ) THEN
161
-
!iLeft = NPts - 1
162
#####
RInt%R = 0.0 ! R( iRight )%R
163
#####
RInt%phi = 0.0 ! R( iRight )%phi
164
-
! WRITE( PRTFile, * ) 'Warning in InterpolateReflectionCoefficient : Refl. Coef. being set to 0 outside tabulated domain'
165
-
! WRITE( PRTFile, * ) 'angle = ', ThetaIntr, 'upper limit = ', R( iRight )%theta
166
-
167
-
ELSE
168
-
! Search for bracketing abscissas: Log2( NPts ) stabs required for a bracket
169
-
170
286678
DO WHILE ( iLeft /= iRight - 1 )
171
143339
iMid = ( iLeft + iRight ) / 2
172
143339*
IF ( R( iMid )%theta > thetaIntr ) THEN
173
27263
iRight = iMid
174
-
ELSE
175
116076
iLeft = iMid
176
-
END IF
177
-
END DO
178
-
179
-
! Linear interpolation for reflection coef
180
-
181
143339*
alpha = ( RInt%theta - R( iLeft )%theta ) / ( R( iRight )%theta - R( iLeft )%theta )
182
143339*
RInt%R = ( 1 - alpha ) * R( iLeft )%R + alpha * R( iRight )%R
183
143339*
RInt%phi = ( 1 - alpha ) * R( iLeft )%phi + alpha * R( iRight )%phi
184
-
185
-
END IF
186
-
187
143339
END SUBROUTINE InterpolateReflectionCoefficient
188
-
189
#####
SUBROUTINE InterpolateIRC( x, f, g, iPower, xTab, fTab, gTab, iTab, NkTab )
190
-
191
-
! Internal reflection coefficient interpolator.
192
-
! Returns f, g, iPower for given x using tabulated values.
193
-
! Uses polynomial interpolation to approximate the function between the tabulated values
194
-
195
-
USE PolyMod
196
-
IMPLICIT NONE
197
-
198
-
INTEGER, PARAMETER :: N = 3 ! order of the polynomial for interpolation
199
-
INTEGER, INTENT( IN ) :: NkTab, iTab( NkTab )
200
-
REAL ( KIND=8 ), INTENT( IN ) :: xTab( NkTab )
201
-
COMPLEX ( KIND=8 ), INTENT( IN ) :: fTab( NkTab ), gTab( NkTab ), x
202
-
INTEGER, INTENT( OUT ) :: iPower
203
-
COMPLEX ( KIND=8 ), INTENT( OUT ) :: f, g
204
-
INTEGER :: i, j, iDel, iLeft, iMid, iRight, NAct
205
-
REAL ( KIND=8 ) :: xReal
206
-
COMPLEX ( KIND=8 ) :: xT( N ), fT( N ), gT( N )
207
-
208
#####
xReal = DBLE( x ) ! taking the real part
209
#####
iLeft = 1
210
#####
iRight = NkTab
211
-
212
-
! Three cases: x left, in, or right of tabulated interval
213
-
214
#####
IF ( xReal < xTab( iLeft ) ) THEN
215
#####
f = fTab( iLeft )
216
#####
g = gTab( iLeft )
217
#####
iPower = iTab( iLeft )
218
#####
ELSE IF( xReal > xTab( iRight ) ) THEN
219
#####
f = fTab( iRight )
220
#####
g = gTab( iRight )
221
#####
iPower = iTab( iRight )
222
-
ELSE
223
-
224
-
! Search for bracketing abscissas:
225
-
! Log base 2 (NPts) stabs required for a bracket
226
-
227
#####
DO WHILE ( iLeft /= iRight-1 )
228
#####
iMid = ( iLeft + iRight )/2
229
#####
IF ( xTab( iMid ) > xReal ) THEN
230
#####
iRight = iMid
231
-
ELSE
232
#####
iLeft = iMid
233
-
END IF
234
-
END DO
235
-
236
-
! Extract the subset for interpolation and scale
237
-
238
#####
iLeft = MAX( iLeft - ( N - 2 ) / 2, 1 )
239
#####
iRight = MiN( iRight + ( N - 1 ) / 2, NkTab )
240
-
241
#####
NAct = iRight - iLeft + 1
242
#####
DO i = 1, NAct
243
#####
j = i + iLeft - 1
244
#####
iDel = iTab( j ) - iTab( iLeft )
245
#####
xT( i ) = xTab( j )
246
#####
fT( i ) = fTab( j ) * 10.0D0 ** iDel
247
#####
gT( i ) = gTab( j ) * 10.0D0 ** iDel
248
-
END DO
249
-
250
-
! Construct the interpolate
251
-
252
#####
f = Poly( x, xT, fT, NAct )
253
#####
g = Poly( x, xT, gT, NAct )
254
#####
iPower = iTab( iLeft )
255
-
END IF
256
-
257
#####
END SUBROUTINE InterpolateIRC
258
#####
END MODULE RefCoef

Legend:

Execution Count: Number of times each line was executed. ##### indicates never executed, - indicates non-executable.