Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

math_lib.f90 @ 5268

Last change on this file since 5268 was 5268, checked in by abarral, 8 days ago
.f90 <-> .F90 depending on cpp key use
File size: 12.5 KB

Line
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	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	integer, dimension(i2-i1+1) :: idx
86	real(wp), dimension(i2-i1+1) :: f_rev, s_rev
87
88	nelm = i2-i1+1
89	if (nelm > 3) then
90	call mrgrnk(s(i1:i2),idx)
91	s_rev = s(idx)
92	f_rev = f(idx)
93	call avint(f_rev(i1:i2),s_rev(i1:i2),(i2-i1+1), &
94	s_rev(i1),s_rev(i2), val)
95	path_integral = val
96	else
97	sumo = 0._wp
98	do j=i1,i2
99	deltah = abs(s(i1+1)-s(i1))
100	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	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	lerror = .false.
198	a = a_in
199	b = b_in
200
201	if ( ntab < 3 ) then
202	call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Ntab is less than 3.')
203	return
204	end if
205
206	do i = 2, ntab
207	if ( xtab(i) <= xtab(i-1) ) then
208	lerror = .true.
209	exit
210	end if
211	end do
212
213	if (lerror) then
214	call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Xtab(i) is not greater than Xtab(i-1).')
215	return
216	end if
217
218	!ds result = 0.0D+00
219	result = 0._wp
220
221	if ( a == b ) then
222	call errorMessage('WARNING(optics/math_lib.f90:AVINT): A=B => integral=0')
223	return
224	end if
225
226	! If B < A, temporarily switch A and B, and store sign.
227	if ( b < a ) then
228	syl = b
229	b = a
230	a = syl
231	ind = -1
232	else
233	syl = a
234	ind = 1
235	end if
236
237	! Bracket A and B between XTAB(ILO) and XTAB(IHI).
238	ilo = 1
239	ihi = ntab
240
241	do i = 1, ntab
242	if ( a <= xtab(i) ) then
243	exit
244	end if
245	ilo = ilo + 1
246	end do
247
248	ilo = max ( 2, ilo )
249	ilo = min ( ilo, ntab - 1 )
250
251	do i = 1, ntab
252	if ( xtab(i) <= b ) then
253	exit
254	end if
255	ihi = ihi - 1
256	end do
257
258	ihi = min ( ihi, ntab - 1 )
259	ihi = max ( ilo, ihi - 1 )
260
261	! Carry out approximate integration from XTAB(ILO) to XTAB(IHI).
262	sum1 = 0._wp
263	!ds sum1 = 0.0D+00
264
265	do i = ilo, ihi
266
267	x1 = xtab(i-1)
268	x2 = xtab(i)
269	x3 = xtab(i+1)
270
271	term1 = ftab(i-1) / ( ( x1 - x2 ) * ( x1 - x3 ) )
272	term2 = ftab(i) / ( ( x2 - x1 ) * ( x2 - x3 ) )
273	term3 = ftab(i+1) / ( ( x3 - x1 ) * ( x3 - x2 ) )
274
275	atemp = term1 + term2 + term3
276
277	btemp = - ( x2 + x3 ) * term1 &
278	- ( x1 + x3 ) * term2 &
279	- ( x1 + x2 ) * term3
280
281	ctemp = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3
282
283	if ( i <= ilo ) then
284	ca = atemp
285	cb = btemp
286	cc = ctemp
287	else
288	ca = 0.5_wp * ( atemp + ca )
289	cb = 0.5_wp * ( btemp + cb )
290	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 * ( x23 - syl3 ) / 3._wp &
297	+ cb * 0.5_wp * ( x22 - syl2 ) + cc * ( x2 - syl )
298	!ds sum1 = sum1 + ca * ( x23 - syl3 ) / 3.0D+00 &
299	!ds + cb * 0.5D+00 * ( x22 - syl2 ) + cc * ( x2 - syl )
300
301	ca = atemp
302	cb = btemp
303	cc = ctemp
304
305	syl = x2
306
307	end do
308
309	result = sum1 + ca * ( b3 - syl3 ) / 3._wp &
310	+ cb * 0.5_wp * ( b2 - syl2 ) + cc * ( b - syl )
311	!ds result = sum1 + ca * ( b3 - syl3 ) / 3.0D+00 &
312	!ds + cb * 0.5D+00 * ( b2 - syl2 ) + cc * ( b - syl )
313
314	! Restore original values of A and B, reverse sign of integral
315	! because of earlier switch.
316	if ( ind /= 1 ) then
317	ind = 1
318	syl = b
319	b = a
320	a = syl
321	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	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	pi = acos(-1._wp)
369	if (x ==int(x)) then
370	if (x > 0.0) then
371	ga=1._wp
372	m1=x-1
373	do k=2,m1
374	ga=ga*k
375	enddo
376	else
377	ga=1._wp+300
378	endif
379	else
380	if (abs(x) > 1.0) then
381	z=abs(x)
382	m=int(z)
383	r=1._wp
384	do k=1,m
385	r=r*(z-k)
386	enddo
387	z=z-m
388	else
389	z=x
390	endif
391	gr=g(26)
392	do k=25,1,-1
393	gr=gr*z+g(k)
394	enddo
395	ga=1._wp/(gr*z)
396	if (abs(x) > 1.0) then
397	ga=ga*r
398	if (x < 0.0) ga=-pi/(xgasin(pi*x))
399	endif
400	endif
401	gamma = ga
402	return
403	end function gamma
404	end module math_lib

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format