Compute parabolic cylinder function Vv(x) for small argument
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | Va |
Order |
||
| real(kind=wp), | intent(in) | :: | x |
Argument |
||
| real(kind=wp), | intent(out) | :: | Pv |
Vv(x) |
subroutine vvsa(Va,x,Pv) real(wp),intent(in) :: Va !! Order real(wp),intent(in) :: x !! Argument real(wp),intent(out) :: Pv !! Vv(x) real(wp) :: a0 , ep , fac , g1 , ga0 , gm , gw , & r , r1 , sv , sv0 , v1 , va0 , & vb0 , vm integer :: m real(wp),parameter :: eps = 1.0e-15_wp ep = exp(-0.25_wp*x*x) va0 = 1.0_wp + 0.5_wp*Va if ( x/=0.0_wp ) then a0 = 2.0_wp**(-0.5_wp*Va)*ep/(twopi) sv = sin(-(Va+0.5_wp)*pi) v1 = -0.5_wp*Va call gamma2(v1,g1) Pv = (sv+1.0_wp)*g1 r = 1.0_wp fac = 1.0_wp do m = 1 , 250 vm = 0.5_wp*(m-Va) call gamma2(vm,gm) r = r*sq2*x/m fac = -fac gw = fac*sv + 1.0_wp r1 = gw*r*gm Pv = Pv + r1 if ( abs(r1/Pv)<eps .and. gw/=0.0_wp ) exit enddo Pv = a0*Pv elseif ( va0<=0.0_wp .and. va0==int(va0) .or. Va==0.0_wp ) then Pv = 0.0_wp else vb0 = -0.5_wp*Va sv0 = sin(va0*pi) call gamma2(va0,ga0) Pv = 2.0_wp**vb0*sv0/ga0 endif end subroutine vvsa