Coverage Report: RefCoef.f90

Generated from GCOV analysis of Fortran source code

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

Legend:

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