lxsd / fvn

Blame view

fvn_fnlib/dbskes.f 1.42 KB
edit raw normal view history
38581db0c   daniau   git-svn-id: https... 
1
2
 
        subroutine dbskes (xnu, x, nin, bke)
  c july 1980 edition.   w. fullerton, c3, los alamos scientific lab.
00055ac08   kwagner   Define some const... 
3
 
        double precision xnu, x, bke(*), bknu1, v, vincr, vend, alnbig,
38581db0c   daniau   git-svn-id: https... 
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
 
       1  d1mach, dlog
        external d1mach
        data alnbig / 0.d0 /
  c
        if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2))
  c
        v = dabs(xnu)
        n = iabs(nin)
  c
        if (v.ge.1.d0) call seteru (29hdbskes  abs(xnu) must be lt 1, 29,
       1  1, 2)
        if (x.le.0.d0) call seteru (17hdbskes  x is le 0, 17, 2, 2)
        if (n.eq.0) call seteru (
       1  41hdbskes  n the number in the sequence is 0, 41, 3, 2)
  c
        call d9knus (v, x, bke(1), bknu1, iswtch)
        if (n.eq.1) return
  c
        vincr = sign (1.0, float(nin))
        direct = vincr
        if (xnu.ne.0.d0) direct = vincr*dsign(1.d0, xnu)
        if (iswtch.eq.1 .and. direct.gt.0.) call seteru (
       1  47hdbskes  x so small bessel k-sub-xnu+1 overflows, 47, 4, 2)
        bke(2) = bknu1
  c
        if (direct.lt.0.) call d9knus (dabs(xnu+vincr), x, bke(2), bknu1,
       1  iswtch)
        if (n.eq.2) return
  c
        vend = dabs (xnu+dble(float(nin))) - 1.0d0
        if (x + (vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)+.306d0)
       1  .gt.alnbig) call seteru (
       2               67hdbskes  x so small or abs(nu) so big that bessel
       3 k-sub-nu overflows, 67, 5, 2)
  c
        v = xnu
        do 10 i=3,n
          v = v + vincr
          bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2)
   10   continue
  c
        return
        end