← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/integrate/quadpack/dqc25c.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/integrate/quadpack/dqc25c.f
 
1
      subroutine dqc25c(f,a,b,c,result,abserr,krul,neval)
 
2
c***begin prologue  dqc25c
 
3
c***date written   810101   (yymmdd)
 
4
c***revision date  830518   (yymmdd)
 
5
c***category no.  h2a2a2,j4
 
6
c***keywords  25-point clenshaw-curtis integration
 
7
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
 
8
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
 
9
c***purpose  to compute i = integral of f*w over (a,b) with
 
10
c            error estimate, where w(x) = 1/(x-c)
 
11
c***description
 
12
c
 
13
c        integration rules for the computation of cauchy
 
14
c        principal value integrals
 
15
c        standard fortran subroutine
 
16
c        double precision version
 
17
c
 
18
c        parameters
 
19
c           f      - double precision
 
20
c                    function subprogram defining the integrand function
 
21
c                    f(x). the actual name for f needs to be declared
 
22
c                    e x t e r n a l  in the driver program.
 
23
c
 
24
c           a      - double precision
 
25
c                    left end point of the integration interval
 
26
c
 
27
c           b      - double precision
 
28
c                    right end point of the integration interval, b.gt.a
 
29
c
 
30
c           c      - double precision
 
31
c                    parameter in the weight function
 
32
c
 
33
c           result - double precision
 
34
c                    approximation to the integral
 
35
c                    result is computed by using a generalized
 
36
c                    clenshaw-curtis method if c lies within ten percent
 
37
c                    of the integration interval. in the other case the
 
38
c                    15-point kronrod rule obtained by optimal addition
 
39
c                    of abscissae to the 7-point gauss rule, is applied.
 
40
c
 
41
c           abserr - double precision
 
42
c                    estimate of the modulus of the absolute error,
 
43
c                    which should equal or exceed abs(i-result)
 
44
c
 
45
c           krul   - integer
 
46
c                    key which is decreased by 1 if the 15-point
 
47
c                    gauss-kronrod scheme has been used
 
48
c
 
49
c           neval  - integer
 
50
c                    number of integrand evaluations
 
51
c
 
52
c.......................................................................
 
53
c***references  (none)
 
54
c***routines called  dqcheb,dqk15w,dqwgtc
 
55
c***end prologue  dqc25c
 
56
c
 
57
      double precision a,abserr,ak22,amom0,amom1,amom2,b,c,cc,centr,
 
58
     *  cheb12,cheb24,dabs,dlog,dqwgtc,f,fval,hlgth,p2,p3,p4,resabs,
 
59
     *  resasc,result,res12,res24,u,x
 
60
      integer i,isym,k,kp,krul,neval
 
61
c
 
62
      dimension x(11),fval(25),cheb12(13),cheb24(25)
 
63
c
 
64
      external f,dqwgtc
 
65
c
 
66
c           the vector x contains the values cos(k*pi/24),
 
67
c           k = 1, ..., 11, to be used for the chebyshev series
 
68
c           expansion of f
 
69
c
 
70
      data x(1) / 0.9914448613 7381041114 4557526928 563d0 /
 
71
      data x(2) / 0.9659258262 8906828674 9743199728 897d0 /
 
72
      data x(3) / 0.9238795325 1128675612 8183189396 788d0 /
 
73
      data x(4) / 0.8660254037 8443864676 3723170752 936d0 /
 
74
      data x(5) / 0.7933533402 9123516457 9776961501 299d0 /
 
75
      data x(6) / 0.7071067811 8654752440 0844362104 849d0 /
 
76
      data x(7) / 0.6087614290 0872063941 6097542898 164d0 /
 
77
      data x(8) / 0.5000000000 0000000000 0000000000 000d0 /
 
78
      data x(9) / 0.3826834323 6508977172 8459984030 399d0 /
 
79
      data x(10) / 0.2588190451 0252076234 8898837624 048d0 /
 
80
      data x(11) / 0.1305261922 2005159154 8406227895 489d0 /
 
81
c
 
82
c           list of major variables
 
83
c           ----------------------
 
84
c           fval   - value of the function f at the points
 
85
c                    cos(k*pi/24),  k = 0, ..., 24
 
86
c           cheb12 - chebyshev series expansion coefficients,
 
87
c                    for the function f, of degree 12
 
88
c           cheb24 - chebyshev series expansion coefficients,
 
89
c                    for the function f, of degree 24
 
90
c           res12  - approximation to the integral corresponding
 
91
c                    to the use of cheb12
 
92
c           res24  - approximation to the integral corresponding
 
93
c                    to the use of cheb24
 
94
c           dqwgtc - external function subprogram defining
 
95
c                    the weight function
 
96
c           hlgth  - half-length of the interval
 
97
c           centr  - mid point of the interval
 
98
c
 
99
c
 
100
c           check the position of c.
 
101
c
 
102
c***first executable statement  dqc25c
 
103
      cc = (0.2d+01*c-b-a)/(b-a)
 
104
      if(dabs(cc).lt.0.11d+01) go to 10
 
105
c
 
106
c           apply the 15-point gauss-kronrod scheme.
 
107
c
 
108
      krul = krul-1
 
109
      call dqk15w(f,dqwgtc,c,p2,p3,p4,kp,a,b,result,abserr,
 
110
     *  resabs,resasc)
 
111
      neval = 15
 
112
      if (resasc.eq.abserr) krul = krul+1
 
113
      go to 50
 
114
c
 
115
c           use the generalized clenshaw-curtis method.
 
116
c
 
117
   10 hlgth = 0.5d+00*(b-a)
 
118
      centr = 0.5d+00*(b+a)
 
119
      neval = 25
 
120
      fval(1) = 0.5d+00*f(hlgth+centr)
 
121
      fval(13) = f(centr)
 
122
      fval(25) = 0.5d+00*f(centr-hlgth)
 
123
      do 20 i=2,12
 
124
        u = hlgth*x(i-1)
 
125
        isym = 26-i
 
126
        fval(i) = f(u+centr)
 
127
        fval(isym) = f(centr-u)
 
128
   20 continue
 
129
c
 
130
c           compute the chebyshev series expansion.
 
131
c
 
132
      call dqcheb(x,fval,cheb12,cheb24)
 
133
c
 
134
c           the modified chebyshev moments are computed by forward
 
135
c           recursion, using amom0 and amom1 as starting values.
 
136
c
 
137
      amom0 = dlog(dabs((0.1d+01-cc)/(0.1d+01+cc)))
 
138
      amom1 = 0.2d+01+cc*amom0
 
139
      res12 = cheb12(1)*amom0+cheb12(2)*amom1
 
140
      res24 = cheb24(1)*amom0+cheb24(2)*amom1
 
141
      do 30 k=3,13
 
142
        amom2 = 0.2d+01*cc*amom1-amom0
 
143
        ak22 = (k-2)*(k-2)
 
144
        if((k/2)*2.eq.k) amom2 = amom2-0.4d+01/(ak22-0.1d+01)
 
145
        res12 = res12+cheb12(k)*amom2
 
146
        res24 = res24+cheb24(k)*amom2
 
147
        amom0 = amom1
 
148
        amom1 = amom2
 
149
   30 continue
 
150
      do 40 k=14,25
 
151
        amom2 = 0.2d+01*cc*amom1-amom0
 
152
        ak22 = (k-2)*(k-2)
 
153
        if((k/2)*2.eq.k) amom2 = amom2-0.4d+01/(ak22-0.1d+01)
 
154
        res24 = res24+cheb24(k)*amom2
 
155
        amom0 = amom1
 
156
        amom1 = amom2
 
157
   40 continue
 
158
      result = res24
 
159
      abserr = dabs(res24-res12)
 
160
   50 return
 
161
      end