Line data Source code
1 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 : ! Copyright (c) 2015, Regents of the University of Colorado
3 : ! All rights reserved.
4 : !
5 : ! Redistribution and use in source and binary forms, with or without modification, are
6 : ! permitted provided that the following conditions are met:
7 : !
8 : ! 1. Redistributions of source code must retain the above copyright notice, this list of
9 : ! conditions and the following disclaimer.
10 : !
11 : ! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 : ! of conditions and the following disclaimer in the documentation and/or other
13 : ! materials provided with the distribution.
14 : !
15 : ! 3. Neither the name of the copyright holder nor the names of its contributors may be
16 : ! used to endorse or promote products derived from this software without specific prior
17 : ! written permission.
18 : !
19 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 : ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 : ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24 : ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 : ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 : ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 : ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 : !
29 : ! History:
30 : ! July 2006: John Haynes - Initial version
31 : ! May 2015: Dustin Swales - Modified for COSPv2.0
32 : !
33 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 : module math_lib
35 : USE COSP_KINDS, ONLY: wp
36 : use mod_cosp_error, ONLY: errorMessage
37 : implicit none
38 :
39 : contains
40 : ! ##########################################################################
41 : ! function PATH_INTEGRAL
42 : ! ##########################################################################
43 0 : function path_integral(f,s,i1,i2)
44 : use m_mrgrnk
45 : use array_lib
46 : implicit none
47 : ! ########################################################################
48 : ! Purpose:
49 : ! evalues the integral (f ds) between f(index=i1) and f(index=i2)
50 : ! using the AVINT procedure
51 : !
52 : ! Inputs:
53 : ! [f] functional values
54 : ! [s] abscissa values
55 : ! [i1] index of lower limit
56 : ! [i2] index of upper limit
57 : !
58 : ! Returns:
59 : ! result of path integral
60 : !
61 : ! Notes:
62 : ! [s] may be in forward or reverse numerical order
63 : !
64 : ! Requires:
65 : ! mrgrnk package
66 : !
67 : ! Created:
68 : ! 02/02/06 John Haynes (haynes@atmos.colostate.edu)
69 : ! ########################################################################
70 :
71 : ! INPUTS
72 : real(wp),intent(in), dimension(:) :: &
73 : f, & ! Functional values
74 : s ! Abscissa values
75 : integer, intent(in) :: &
76 : i1, & ! Index of lower limit
77 : i2 ! Index of upper limit
78 :
79 : ! OUTPUTS
80 : real(wp) :: path_integral
81 :
82 : ! Internal variables
83 : real(wp) :: sumo, deltah, val
84 : integer :: nelm, j
85 0 : integer, dimension(i2-i1+1) :: idx
86 0 : real(wp), dimension(i2-i1+1) :: f_rev, s_rev
87 :
88 0 : nelm = i2-i1+1
89 0 : if (nelm > 3) then
90 0 : call mrgrnk(s(i1:i2),idx)
91 0 : s_rev = s(idx)
92 0 : f_rev = f(idx)
93 0 : call avint(f_rev(i1:i2),s_rev(i1:i2),(i2-i1+1), &
94 0 : s_rev(i1),s_rev(i2), val)
95 0 : path_integral = val
96 : else
97 : sumo = 0._wp
98 0 : do j=i1,i2
99 0 : deltah = abs(s(i1+1)-s(i1))
100 0 : sumo = sumo + f(j)*deltah
101 : enddo
102 : path_integral = sumo
103 : endif
104 :
105 : return
106 : end function path_integral
107 :
108 : ! ##########################################################################
109 : ! subroutine AVINT
110 : ! ##########################################################################
111 55900936 : subroutine avint ( ftab, xtab, ntab, a_in, b_in, result )
112 : implicit none
113 : ! ########################################################################
114 : ! Purpose:
115 : ! estimate the integral of unevenly spaced data
116 : !
117 : ! Inputs:
118 : ! [ftab] functional values
119 : ! [xtab] abscissa values
120 : ! [ntab] number of elements of [ftab] and [xtab]
121 : ! [a] lower limit of integration
122 : ! [b] upper limit of integration
123 : !
124 : ! Outputs:
125 : ! [result] approximate value of integral
126 : !
127 : ! Reference:
128 : ! From SLATEC libraries, in public domain
129 : !
130 : !***********************************************************************
131 : !
132 : ! AVINT estimates the integral of unevenly spaced data.
133 : !
134 : ! Discussion:
135 : !
136 : ! The method uses overlapping parabolas and smoothing.
137 : !
138 : ! Modified:
139 : !
140 : ! 30 October 2000
141 : ! 4 January 2008, A. Bodas-Salcedo. Error control for XTAB taken out of
142 : ! loop to allow vectorization.
143 : !
144 : ! Reference:
145 : !
146 : ! Philip Davis and Philip Rabinowitz,
147 : ! Methods of Numerical Integration,
148 : ! Blaisdell Publishing, 1967.
149 : !
150 : ! P E Hennion,
151 : ! Algorithm 77,
152 : ! Interpolation, Differentiation and Integration,
153 : ! Communications of the Association for Computing Machinery,
154 : ! Volume 5, page 96, 1962.
155 : !
156 : ! Parameters:
157 : !
158 : ! Input, real ( kind = 8 ) FTAB(NTAB), the function values,
159 : ! FTAB(I) = F(XTAB(I)).
160 : !
161 : ! Input, real ( kind = 8 ) XTAB(NTAB), the abscissas at which the
162 : ! function values are given. The XTAB's must be distinct
163 : ! and in ascending order.
164 : !
165 : ! Input, integer NTAB, the number of entries in FTAB and
166 : ! XTAB. NTAB must be at least 3.
167 : !
168 : ! Input, real ( kind = 8 ) A, the lower limit of integration. A should
169 : ! be, but need not be, near one endpoint of the interval
170 : ! (X(1), X(NTAB)).
171 : !
172 : ! Input, real ( kind = 8 ) B, the upper limit of integration. B should
173 : ! be, but need not be, near one endpoint of the interval
174 : ! (X(1), X(NTAB)).
175 : !
176 : ! Output, real ( kind = 8 ) RESULT, the approximate value of the integral.
177 : ! ##########################################################################
178 :
179 : ! INPUTS
180 : integer,intent(in) :: &
181 : ntab ! Number of elements of [ftab] and [xtab]
182 : real(wp),intent(in) :: &
183 : a_in, & ! Lower limit of integration
184 : b_in ! Upper limit of integration
185 : real(wp),intent(in),dimension(ntab) :: &
186 : ftab, & ! Functional values
187 : xtab ! Abscissa value
188 :
189 : ! OUTPUTS
190 : real(wp),intent(out) :: result ! Approximate value of integral
191 :
192 : ! Internal varaibles
193 : real(wp) :: a, atemp, b, btemp,ca,cb,cc,ctemp,sum1,syl,term1,term2,term3,x1,x2,x3
194 : integer :: i,ihi,ilo,ind
195 : logical :: lerror
196 :
197 55900936 : lerror = .false.
198 55900936 : a = a_in
199 55900936 : b = b_in
200 :
201 55900936 : if ( ntab < 3 ) then
202 0 : call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Ntab is less than 3.')
203 0 : return
204 : end if
205 :
206 4751579560 : do i = 2, ntab
207 4751579560 : if ( xtab(i) <= xtab(i-1) ) then
208 : lerror = .true.
209 : exit
210 : end if
211 : end do
212 :
213 55900936 : if (lerror) then
214 0 : call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Xtab(i) is not greater than Xtab(i-1).')
215 0 : return
216 : end if
217 :
218 : !ds result = 0.0D+00
219 55900936 : result = 0._wp
220 :
221 55900936 : if ( a == b ) then
222 0 : call errorMessage('WARNING(optics/math_lib.f90:AVINT): A=B => integral=0')
223 0 : return
224 : end if
225 :
226 : ! If B < A, temporarily switch A and B, and store sign.
227 55900936 : if ( b < a ) then
228 : syl = b
229 : b = a
230 : a = syl
231 : ind = -1
232 : else
233 55900936 : syl = a
234 55900936 : ind = 1
235 : end if
236 :
237 : ! Bracket A and B between XTAB(ILO) and XTAB(IHI).
238 55900936 : ilo = 1
239 55900936 : ihi = ntab
240 :
241 55900936 : do i = 1, ntab
242 55900936 : if ( a <= xtab(i) ) then
243 : exit
244 : end if
245 55900936 : ilo = ilo + 1
246 : end do
247 :
248 55900936 : ilo = max ( 2, ilo )
249 55900936 : ilo = min ( ilo, ntab - 1 )
250 :
251 55900936 : do i = 1, ntab
252 55900936 : if ( xtab(i) <= b ) then
253 : exit
254 : end if
255 55900936 : ihi = ihi - 1
256 : end do
257 :
258 55900936 : ihi = min ( ihi, ntab - 1 )
259 55900936 : ihi = max ( ilo, ihi - 1 )
260 :
261 : ! Carry out approximate integration from XTAB(ILO) to XTAB(IHI).
262 55900936 : sum1 = 0._wp
263 : !ds sum1 = 0.0D+00
264 :
265 4639777688 : do i = ilo, ihi
266 :
267 4583876752 : x1 = xtab(i-1)
268 4583876752 : x2 = xtab(i)
269 4583876752 : x3 = xtab(i+1)
270 :
271 4583876752 : term1 = ftab(i-1) / ( ( x1 - x2 ) * ( x1 - x3 ) )
272 4583876752 : term2 = ftab(i) / ( ( x2 - x1 ) * ( x2 - x3 ) )
273 4583876752 : term3 = ftab(i+1) / ( ( x3 - x1 ) * ( x3 - x2 ) )
274 :
275 4583876752 : atemp = term1 + term2 + term3
276 :
277 : btemp = - ( x2 + x3 ) * term1 &
278 : - ( x1 + x3 ) * term2 &
279 4583876752 : - ( x1 + x2 ) * term3
280 :
281 4583876752 : ctemp = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3
282 :
283 4583876752 : if ( i <= ilo ) then
284 : ca = atemp
285 : cb = btemp
286 : cc = ctemp
287 : else
288 4527975816 : ca = 0.5_wp * ( atemp + ca )
289 4527975816 : cb = 0.5_wp * ( btemp + cb )
290 4527975816 : cc = 0.5_wp * ( ctemp + cc )
291 : !ds ca = 0.5D+00 * ( atemp + ca )
292 : !ds cb = 0.5D+00 * ( btemp + cb )
293 : !ds cc = 0.5D+00 * ( ctemp + cc )
294 : end if
295 :
296 : sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3._wp &
297 4583876752 : + cb * 0.5_wp * ( x2**2 - syl**2 ) + cc * ( x2 - syl )
298 : !ds sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3.0D+00 &
299 : !ds + cb * 0.5D+00 * ( x2**2 - syl**2 ) + cc * ( x2 - syl )
300 :
301 4583876752 : ca = atemp
302 4583876752 : cb = btemp
303 4583876752 : cc = ctemp
304 :
305 4639777688 : syl = x2
306 :
307 : end do
308 :
309 : result = sum1 + ca * ( b**3 - syl**3 ) / 3._wp &
310 55900936 : + cb * 0.5_wp * ( b**2 - syl**2 ) + cc * ( b - syl )
311 : !ds result = sum1 + ca * ( b**3 - syl**3 ) / 3.0D+00 &
312 : !ds + cb * 0.5D+00 * ( b**2 - syl**2 ) + cc * ( b - syl )
313 :
314 : ! Restore original values of A and B, reverse sign of integral
315 : ! because of earlier switch.
316 55900936 : if ( ind /= 1 ) then
317 0 : ind = 1
318 0 : syl = b
319 0 : b = a
320 0 : a = syl
321 0 : result = -result
322 : end if
323 :
324 : return
325 : end subroutine avint
326 : ! ######################################################################################
327 : ! SUBROUTINE gamma
328 : ! Purpose:
329 : ! Returns the gamma function
330 : !
331 : ! Input:
332 : ! [x] value to compute gamma function of
333 : !
334 : ! Returns:
335 : ! gamma(x)
336 : !
337 : ! Coded:
338 : ! 02/02/06 John Haynes (haynes@atmos.colostate.edu)
339 : ! (original code of unknown origin)
340 : ! ######################################################################################
341 28140854 : function gamma(x)
342 : ! Inputs
343 : real(wp), intent(in) :: x
344 :
345 : ! Outputs
346 : real(wp) :: gamma
347 :
348 : ! Local variables
349 : real(wp) :: pi,ga,z,r,gr
350 : integer :: k,m1,m
351 :
352 : ! Parameters
353 : real(wp),dimension(26),parameter :: &
354 : g = (/1.0,0.5772156649015329, -0.6558780715202538, -0.420026350340952e-1, &
355 : 0.1665386113822915,-0.421977345555443e-1,-0.96219715278770e-2, &
356 : 0.72189432466630e-2,-0.11651675918591e-2, -0.2152416741149e-3, &
357 : 0.1280502823882e-3, -0.201348547807e-4, -0.12504934821e-5, 0.11330272320e-5, &
358 : -0.2056338417e-6, 0.61160950e-8,0.50020075e-8, -0.11812746e-8, 0.1043427e-9, &
359 : 0.77823e-11, -0.36968e-11, 0.51e-12, -0.206e-13, -0.54e-14, 0.14e-14, 0.1e-15/)
360 : !ds real(wp),dimension(26),parameter :: &
361 : !ds g = (/1.0d0,0.5772156649015329d0, -0.6558780715202538d0, -0.420026350340952d-1, &
362 : !ds 0.1665386113822915d0,-0.421977345555443d-1,-0.96219715278770d-2, &
363 : !ds 0.72189432466630d-2,-0.11651675918591d-2, -0.2152416741149d-3, &
364 : !ds 0.1280502823882d-3, -0.201348547807d-4, -0.12504934821d-5, 0.11330272320d-5, &
365 : !ds -0.2056338417d-6, 0.61160950d-8,0.50020075d-8, -0.11812746d-8, 0.1043427d-9, &
366 : !ds 0.77823d-11, -0.36968d-11, 0.51d-12, -0.206d-13, -0.54d-14, 0.14d-14, 0.1d-15/)
367 :
368 28140854 : pi = acos(-1._wp)
369 28140854 : if (x ==int(x)) then
370 27932981 : if (x > 0.0) then
371 27932981 : ga=1._wp
372 27932981 : m1=x-1
373 84006816 : do k=2,m1
374 84006816 : ga=ga*k
375 : enddo
376 : else
377 : ga=1._wp+300
378 : endif
379 : else
380 207873 : if (abs(x) > 1.0) then
381 207873 : z=abs(x)
382 207873 : m=int(z)
383 207873 : r=1._wp
384 1039365 : do k=1,m
385 1039365 : r=r*(z-k)
386 : enddo
387 207873 : z=z-m
388 : else
389 : z=x
390 : endif
391 207873 : gr=g(26)
392 5404698 : do k=25,1,-1
393 5404698 : gr=gr*z+g(k)
394 : enddo
395 207873 : ga=1._wp/(gr*z)
396 207873 : if (abs(x) > 1.0) then
397 207873 : ga=ga*r
398 207873 : if (x < 0.0) ga=-pi/(x*ga*sin(pi*x))
399 : endif
400 : endif
401 28140854 : gamma = ga
402 : return
403 : end function gamma
404 : end module math_lib
|
