! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Copyright (c) 2015, Regents of the University of Colorado ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without modification, are ! permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of ! conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list ! of conditions and the following disclaimer in the documentation and/or other ! materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be ! used to endorse or promote products derived from this software without specific prior ! written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! History: ! July 2006: John Haynes - Initial version ! May 2015: Dustin Swales - Modified for COSPv2.0 ! ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% module math_lib USE COSP_KINDS, ONLY: wp use mod_cosp_error, ONLY: errorMessage implicit none contains ! ########################################################################## ! function PATH_INTEGRAL ! ########################################################################## function path_integral(f,s,i1,i2) use m_mrgrnk use array_lib implicit none ! ######################################################################## ! Purpose: ! evalues the integral (f ds) between f(index=i1) and f(index=i2) ! using the AVINT procedure ! ! Inputs: ! [f] functional values ! [s] abscissa values ! [i1] index of lower limit ! [i2] index of upper limit ! ! Returns: ! result of path integral ! ! Notes: ! [s] may be in forward or reverse numerical order ! ! Requires: ! mrgrnk package ! ! Created: ! 02/02/06 John Haynes (haynes@atmos.colostate.edu) ! ######################################################################## ! INPUTS real(wp),intent(in), dimension(:) :: & f, & ! Functional values s ! Abscissa values integer, intent(in) :: & i1, & ! Index of lower limit i2 ! Index of upper limit ! OUTPUTS real(wp) :: path_integral ! Internal variables real(wp) :: sumo, deltah, val integer :: nelm, j integer, dimension(i2-i1+1) :: idx real(wp), dimension(i2-i1+1) :: f_rev, s_rev nelm = i2-i1+1 if (nelm > 3) then call mrgrnk(s(i1:i2),idx) s_rev = s(idx) f_rev = f(idx) call avint(f_rev(i1:i2),s_rev(i1:i2),(i2-i1+1), & s_rev(i1),s_rev(i2), val) path_integral = val else sumo = 0._wp do j=i1,i2 deltah = abs(s(i1+1)-s(i1)) sumo = sumo + f(j)*deltah enddo path_integral = sumo endif return end function path_integral ! ########################################################################## ! subroutine AVINT ! ########################################################################## subroutine avint ( ftab, xtab, ntab, a_in, b_in, result ) implicit none ! ######################################################################## ! Purpose: ! estimate the integral of unevenly spaced data ! ! Inputs: ! [ftab] functional values ! [xtab] abscissa values ! [ntab] number of elements of [ftab] and [xtab] ! [a] lower limit of integration ! [b] upper limit of integration ! ! Outputs: ! [result] approximate value of integral ! ! Reference: ! From SLATEC libraries, in public domain ! !*********************************************************************** ! ! AVINT estimates the integral of unevenly spaced data. ! ! Discussion: ! ! The method uses overlapping parabolas and smoothing. ! ! Modified: ! ! 30 October 2000 ! 4 January 2008, A. Bodas-Salcedo. Error control for XTAB taken out of ! loop to allow vectorization. ! ! Reference: ! ! Philip Davis and Philip Rabinowitz, ! Methods of Numerical Integration, ! Blaisdell Publishing, 1967. ! ! P E Hennion, ! Algorithm 77, ! Interpolation, Differentiation and Integration, ! Communications of the Association for Computing Machinery, ! Volume 5, page 96, 1962. ! ! Parameters: ! ! Input, real ( kind = 8 ) FTAB(NTAB), the function values, ! FTAB(I) = F(XTAB(I)). ! ! Input, real ( kind = 8 ) XTAB(NTAB), the abscissas at which the ! function values are given. The XTAB's must be distinct ! and in ascending order. ! ! Input, integer NTAB, the number of entries in FTAB and ! XTAB. NTAB must be at least 3. ! ! Input, real ( kind = 8 ) A, the lower limit of integration. A should ! be, but need not be, near one endpoint of the interval ! (X(1), X(NTAB)). ! ! Input, real ( kind = 8 ) B, the upper limit of integration. B should ! be, but need not be, near one endpoint of the interval ! (X(1), X(NTAB)). ! ! Output, real ( kind = 8 ) RESULT, the approximate value of the integral. ! ########################################################################## ! INPUTS integer,intent(in) :: & ntab ! Number of elements of [ftab] and [xtab] real(wp),intent(in) :: & a_in, & ! Lower limit of integration b_in ! Upper limit of integration real(wp),intent(in),dimension(ntab) :: & ftab, & ! Functional values xtab ! Abscissa value ! OUTPUTS real(wp),intent(out) :: result ! Approximate value of integral ! Internal varaibles real(wp) :: a, atemp, b, btemp,ca,cb,cc,ctemp,sum1,syl,term1,term2,term3,x1,x2,x3 integer :: i,ihi,ilo,ind logical :: lerror lerror = .false. a = a_in b = b_in if ( ntab < 3 ) then call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Ntab is less than 3.') return end if do i = 2, ntab if ( xtab(i) <= xtab(i-1) ) then lerror = .true. exit end if end do if (lerror) then call errorMessage('FATAL ERROR(optics/math_lib.f90:AVINT): Xtab(i) is not greater than Xtab(i-1).') return end if !ds result = 0.0D+00 result = 0._wp if ( a == b ) then call errorMessage('WARNING(optics/math_lib.f90:AVINT): A=B => integral=0') return end if ! If B < A, temporarily switch A and B, and store sign. if ( b < a ) then syl = b b = a a = syl ind = -1 else syl = a ind = 1 end if ! Bracket A and B between XTAB(ILO) and XTAB(IHI). ilo = 1 ihi = ntab do i = 1, ntab if ( a <= xtab(i) ) then exit end if ilo = ilo + 1 end do ilo = max ( 2, ilo ) ilo = min ( ilo, ntab - 1 ) do i = 1, ntab if ( xtab(i) <= b ) then exit end if ihi = ihi - 1 end do ihi = min ( ihi, ntab - 1 ) ihi = max ( ilo, ihi - 1 ) ! Carry out approximate integration from XTAB(ILO) to XTAB(IHI). sum1 = 0._wp !ds sum1 = 0.0D+00 do i = ilo, ihi x1 = xtab(i-1) x2 = xtab(i) x3 = xtab(i+1) term1 = ftab(i-1) / ( ( x1 - x2 ) * ( x1 - x3 ) ) term2 = ftab(i) / ( ( x2 - x1 ) * ( x2 - x3 ) ) term3 = ftab(i+1) / ( ( x3 - x1 ) * ( x3 - x2 ) ) atemp = term1 + term2 + term3 btemp = - ( x2 + x3 ) * term1 & - ( x1 + x3 ) * term2 & - ( x1 + x2 ) * term3 ctemp = x2 * x3 * term1 + x1 * x3 * term2 + x1 * x2 * term3 if ( i <= ilo ) then ca = atemp cb = btemp cc = ctemp else ca = 0.5_wp * ( atemp + ca ) cb = 0.5_wp * ( btemp + cb ) cc = 0.5_wp * ( ctemp + cc ) !ds ca = 0.5D+00 * ( atemp + ca ) !ds cb = 0.5D+00 * ( btemp + cb ) !ds cc = 0.5D+00 * ( ctemp + cc ) end if sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3._wp & + cb * 0.5_wp * ( x2**2 - syl**2 ) + cc * ( x2 - syl ) !ds sum1 = sum1 + ca * ( x2**3 - syl**3 ) / 3.0D+00 & !ds + cb * 0.5D+00 * ( x2**2 - syl**2 ) + cc * ( x2 - syl ) ca = atemp cb = btemp cc = ctemp syl = x2 end do result = sum1 + ca * ( b**3 - syl**3 ) / 3._wp & + cb * 0.5_wp * ( b**2 - syl**2 ) + cc * ( b - syl ) !ds result = sum1 + ca * ( b**3 - syl**3 ) / 3.0D+00 & !ds + cb * 0.5D+00 * ( b**2 - syl**2 ) + cc * ( b - syl ) ! Restore original values of A and B, reverse sign of integral ! because of earlier switch. if ( ind /= 1 ) then ind = 1 syl = b b = a a = syl result = -result end if return end subroutine avint ! ###################################################################################### ! SUBROUTINE gamma ! Purpose: ! Returns the gamma function ! ! Input: ! [x] value to compute gamma function of ! ! Returns: ! gamma(x) ! ! Coded: ! 02/02/06 John Haynes (haynes@atmos.colostate.edu) ! (original code of unknown origin) ! ###################################################################################### function gamma(x) ! Inputs real(wp), intent(in) :: x ! Outputs real(wp) :: gamma ! Local variables real(wp) :: pi,ga,z,r,gr integer :: k,m1,m ! Parameters real(wp),dimension(26),parameter :: & g = (/1.0,0.5772156649015329, -0.6558780715202538, -0.420026350340952e-1, & 0.1665386113822915,-0.421977345555443e-1,-0.96219715278770e-2, & 0.72189432466630e-2,-0.11651675918591e-2, -0.2152416741149e-3, & 0.1280502823882e-3, -0.201348547807e-4, -0.12504934821e-5, 0.11330272320e-5, & -0.2056338417e-6, 0.61160950e-8,0.50020075e-8, -0.11812746e-8, 0.1043427e-9, & 0.77823e-11, -0.36968e-11, 0.51e-12, -0.206e-13, -0.54e-14, 0.14e-14, 0.1e-15/) !ds real(wp),dimension(26),parameter :: & !ds g = (/1.0d0,0.5772156649015329d0, -0.6558780715202538d0, -0.420026350340952d-1, & !ds 0.1665386113822915d0,-0.421977345555443d-1,-0.96219715278770d-2, & !ds 0.72189432466630d-2,-0.11651675918591d-2, -0.2152416741149d-3, & !ds 0.1280502823882d-3, -0.201348547807d-4, -0.12504934821d-5, 0.11330272320d-5, & !ds -0.2056338417d-6, 0.61160950d-8,0.50020075d-8, -0.11812746d-8, 0.1043427d-9, & !ds 0.77823d-11, -0.36968d-11, 0.51d-12, -0.206d-13, -0.54d-14, 0.14d-14, 0.1d-15/) pi = acos(-1._wp) if (x ==int(x)) then if (x > 0.0) then ga=1._wp m1=x-1 do k=2,m1 ga=ga*k enddo else ga=1._wp+300 endif else if (abs(x) > 1.0) then z=abs(x) m=int(z) r=1._wp do k=1,m r=r*(z-k) enddo z=z-m else z=x endif gr=g(26) do k=25,1,-1 gr=gr*z+g(k) enddo ga=1._wp/(gr*z) if (abs(x) > 1.0) then ga=ga*r if (x < 0.0) ga=-pi/(x*ga*sin(pi*x)) endif endif gamma = ga return end function gamma end module math_lib