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

« back to all changes in this revision

Viewing changes to scipy/integrate/quadpack/dqmomo.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
* 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:
 
1
      subroutine dqmomo(alfa,beta,ri,rj,rg,rh,integr)
 
2
c***begin prologue  dqmomo
 
3
c***date written   820101   (yymmdd)
 
4
c***revision date  830518   (yymmdd)
 
5
c***category no.  h2a2a1,c3a2
 
6
c***keywords  modified chebyshev moments
 
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  this routine computes modified chebsyshev moments. the k-th
 
10
c            modified chebyshev moment is defined as the integral over
 
11
c            (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev
 
12
c            polynomial of degree k.
 
13
c***description
 
14
c
 
15
c        modified chebyshev moments
 
16
c        standard fortran subroutine
 
17
c        double precision version
 
18
c
 
19
c        parameters
 
20
c           alfa   - double precision
 
21
c                    parameter in the weight function w(x), alfa.gt.(-1)
 
22
c
 
23
c           beta   - double precision
 
24
c                    parameter in the weight function w(x), beta.gt.(-1)
 
25
c
 
26
c           ri     - double precision
 
27
c                    vector of dimension 25
 
28
c                    ri(k) is the integral over (-1,1) of
 
29
c                    (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
 
30
c
 
31
c           rj     - double precision
 
32
c                    vector of dimension 25
 
33
c                    rj(k) is the integral over (-1,1) of
 
34
c                    (1-x)**beta*t(k-1,x), k = 1, ..., 25.
 
35
c
 
36
c           rg     - double precision
 
37
c                    vector of dimension 25
 
38
c                    rg(k) is the integral over (-1,1) of
 
39
c                    (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
 
40
c
 
41
c           rh     - double precision
 
42
c                    vector of dimension 25
 
43
c                    rh(k) is the integral over (-1,1) of
 
44
c                    (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
 
45
c
 
46
c           integr - integer
 
47
c                    input parameter indicating the modified
 
48
c                    moments to be computed
 
49
c                    integr = 1 compute ri, rj
 
50
c                           = 2 compute ri, rj, rg
 
51
c                           = 3 compute ri, rj, rh
 
52
c                           = 4 compute ri, rj, rg, rh
 
53
c
 
54
c***references  (none)
 
55
c***routines called  (none)
 
56
c***end prologue  dqmomo
 
57
c
 
58
      double precision alfa,alfp1,alfp2,an,anm1,beta,betp1,betp2,ralf,
 
59
     *  rbet,rg,rh,ri,rj
 
60
      integer i,im1,integr
 
61
c
 
62
      dimension rg(25),rh(25),ri(25),rj(25)
 
63
c
 
64
c
 
65
c***first executable statement  dqmomo
 
66
      alfp1 = alfa+0.1d+01
 
67
      betp1 = beta+0.1d+01
 
68
      alfp2 = alfa+0.2d+01
 
69
      betp2 = beta+0.2d+01
 
70
      ralf = 0.2d+01**alfp1
 
71
      rbet = 0.2d+01**betp1
 
72
c
 
73
c           compute ri, rj using a forward recurrence relation.
 
74
c
 
75
      ri(1) = ralf/alfp1
 
76
      rj(1) = rbet/betp1
 
77
      ri(2) = ri(1)*alfa/alfp2
 
78
      rj(2) = rj(1)*beta/betp2
 
79
      an = 0.2d+01
 
80
      anm1 = 0.1d+01
 
81
      do 20 i=3,25
 
82
        ri(i) = -(ralf+an*(an-alfp2)*ri(i-1))/(anm1*(an+alfp1))
 
83
        rj(i) = -(rbet+an*(an-betp2)*rj(i-1))/(anm1*(an+betp1))
 
84
        anm1 = an
 
85
        an = an+0.1d+01
 
86
   20 continue
 
87
      if(integr.eq.1) go to 70
 
88
      if(integr.eq.3) go to 40
 
89
c
 
90
c           compute rg using a forward recurrence relation.
 
91
c
 
92
      rg(1) = -ri(1)/alfp1
 
93
      rg(2) = -(ralf+ralf)/(alfp2*alfp2)-rg(1)
 
94
      an = 0.2d+01
 
95
      anm1 = 0.1d+01
 
96
      im1 = 2
 
97
      do 30 i=3,25
 
98
        rg(i) = -(an*(an-alfp2)*rg(im1)-an*ri(im1)+anm1*ri(i))/
 
99
     *  (anm1*(an+alfp1))
 
100
        anm1 = an
 
101
        an = an+0.1d+01
 
102
        im1 = i
 
103
   30 continue
 
104
      if(integr.eq.2) go to 70
 
105
c
 
106
c           compute rh using a forward recurrence relation.
 
107
c
 
108
   40 rh(1) = -rj(1)/betp1
 
109
      rh(2) = -(rbet+rbet)/(betp2*betp2)-rh(1)
 
110
      an = 0.2d+01
 
111
      anm1 = 0.1d+01
 
112
      im1 = 2
 
113
      do 50 i=3,25
 
114
        rh(i) = -(an*(an-betp2)*rh(im1)-an*rj(im1)+
 
115
     *  anm1*rj(i))/(anm1*(an+betp1))
 
116
        anm1 = an
 
117
        an = an+0.1d+01
 
118
        im1 = i
 
119
   50 continue
 
120
      do 60 i=2,25,2
 
121
        rh(i) = -rh(i)
 
122
   60 continue
 
123
   70 do 80 i=2,25,2
 
124
        rj(i) = -rj(i)
 
125
   80 continue
 
126
   90 return
 
127
      end