← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpbfout.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/interpolate/fitpack/fpbfout.f
 
1
      subroutine fpbfou(t,n,par,ress,resc)
 
2
c  subroutine fpbfou calculates the integrals
 
3
c                    /t(n-3)
 
4
c    ress(j) =      !        nj,4(x)*sin(par*x) dx    and
 
5
c              t(4)/
 
6
c                    /t(n-3)
 
7
c    resc(j) =      !        nj,4(x)*cos(par*x) dx ,  j=1,2,...n-4
 
8
c              t(4)/
 
9
c  where nj,4(x) denotes the cubic b-spline defined on the knots
 
10
c  t(j),t(j+1),...,t(j+4).
 
11
c
 
12
c  calling sequence:
 
13
c     call fpbfou(t,n,par,ress,resc)
 
14
c
 
15
c  input parameters:
 
16
c    t    : real array,length n, containing the knots.
 
17
c    n    : integer, containing the number of knots.
 
18
c    par  : real, containing the value of the parameter par.
 
19
c
 
20
c  output parameters:
 
21
c    ress  : real array,length n, containing the integrals ress(j).
 
22
c    resc  : real array,length n, containing the integrals resc(j).
 
23
c
 
24
c  restrictions:
 
25
c    n >= 10, t(4) < t(5) < ... < t(n-4) < t(n-3).
 
26
c  ..
 
27
c  ..scalar arguments..
 
28
      integer n
 
29
      real*8 par
 
30
c  ..array arguments..
 
31
      real*8 t(n),ress(n),resc(n)
 
32
c  ..local scalars..
 
33
      integer i,ic,ipj,is,j,jj,jp1,jp4,k,li,lj,ll,nmj,nm3,nm7
 
34
      real*8 ak,beta,con1,con2,c1,c2,delta,eps,fac,f1,f2,f3,one,quart,
 
35
     * sign,six,s1,s2,term
 
36
c  ..local arrays..
 
37
      real*8 co(5),si(5),hs(5),hc(5),rs(3),rc(3)
 
38
c  ..function references..
 
39
      real*8 cos,sin,abs
 
40
c  ..
 
41
c  initialization.
 
42
      one = 0.1e+01
 
43
      six = 0.6e+01
 
44
      eps = 0.1e-07
 
45
      quart = 0.25e0
 
46
      con1 = 0.5e-01
 
47
      con2 = 0.12e+03
 
48
      nm3 = n-3
 
49
      nm7 = n-7
 
50
      if(par.ne.0.) term = six/par
 
51
      beta = par*t(4)
 
52
      co(1) = cos(beta)
 
53
      si(1) = sin(beta)
 
54
c  calculate the integrals ress(j) and resc(j), j=1,2,3 by setting up
 
55
c  a divided difference table.
 
56
      do 30 j=1,3
 
57
        jp1 = j+1
 
58
        jp4 = j+4
 
59
        beta = par*t(jp4)
 
60
        co(jp1) = cos(beta)
 
61
        si(jp1) = sin(beta)
 
62
        call fpcsin(t(4),t(jp4),par,si(1),co(1),si(jp1),co(jp1),
 
63
     *  rs(j),rc(j))
 
64
        i = 5-j
 
65
        hs(i) = 0.
 
66
        hc(i) = 0.
 
67
        do 10 jj=1,j
 
68
          ipj = i+jj
 
69
          hs(ipj) = rs(jj)
 
70
          hc(ipj) = rc(jj)
 
71
  10    continue
 
72
        do 20 jj=1,3
 
73
          if(i.lt.jj) i = jj
 
74
          k = 5
 
75
          li = jp4
 
76
          do 20 ll=i,4
 
77
            lj = li-jj
 
78
            fac = t(li)-t(lj)
 
79
            hs(k) = (hs(k)-hs(k-1))/fac
 
80
            hc(k) = (hc(k)-hc(k-1))/fac
 
81
            k = k-1
 
82
            li = li-1
 
83
  20    continue
 
84
        ress(j) = hs(5)-hs(4)
 
85
        resc(j) = hc(5)-hc(4)
 
86
  30  continue
 
87
      if(nm7.lt.4) go to 160
 
88
c  calculate the integrals ress(j) and resc(j),j=4,5,...,n-7.
 
89
      do 150 j=4,nm7
 
90
        jp4 = j+4
 
91
        beta = par*t(jp4)
 
92
        co(5) = cos(beta)
 
93
        si(5) = sin(beta)
 
94
        delta = t(jp4)-t(j)
 
95
c  the way of computing ress(j) and resc(j) depends on the value of
 
96
c  beta = par*(t(j+4)-t(j)).
 
97
        beta = delta*par
 
98
        if(abs(beta).le.one) go to 60
 
99
c  if !beta! > 1 the integrals are calculated by setting up a divided
 
100
c  difference table.
 
101
        do 40 k=1,5
 
102
          hs(k) = si(k)
 
103
          hc(k) = co(k)
 
104
  40    continue
 
105
        do 50 jj=1,3
 
106
          k = 5
 
107
          li = jp4
 
108
          do 50 ll=jj,4
 
109
            lj = li-jj
 
110
            fac = par*(t(li)-t(lj))
 
111
            hs(k) = (hs(k)-hs(k-1))/fac
 
112
            hc(k) = (hc(k)-hc(k-1))/fac
 
113
            k = k-1
 
114
            li = li-1
 
115
  50    continue
 
116
        s2 = (hs(5)-hs(4))*term
 
117
        c2 = (hc(5)-hc(4))*term
 
118
        go to 130
 
119
c  if !beta! <= 1 the integrals are calculated by evaluating a series
 
120
c  expansion.
 
121
  60    f3 = 0.
 
122
        do 70 i=1,4
 
123
          ipj = i+j
 
124
          hs(i) = par*(t(ipj)-t(j))
 
125
          hc(i) = hs(i)
 
126
          f3 = f3+hs(i)
 
127
  70    continue
 
128
        f3 = f3*con1
 
129
        c1 = quart
 
130
        s1 = f3
 
131
        if(abs(f3).le.eps) go to 120
 
132
        sign = one
 
133
        fac = con2
 
134
        k = 5
 
135
        is = 0
 
136
        do 110 ic=1,20
 
137
          k = k+1
 
138
          ak = k
 
139
          fac = fac*ak
 
140
          f1 = 0.
 
141
          f3 = 0.
 
142
          do 80 i=1,4
 
143
            f1 = f1+hc(i)
 
144
            f2 = f1*hs(i)
 
145
            hc(i) = f2
 
146
            f3 = f3+f2
 
147
  80      continue
 
148
          f3 = f3*six/fac
 
149
          if(is.eq.0) go to 90
 
150
          is = 0
 
151
          s1 = s1+f3*sign
 
152
          go to 100
 
153
  90      sign = -sign
 
154
          is = 1
 
155
          c1 = c1+f3*sign
 
156
 100      if(abs(f3).le.eps) go to 120
 
157
 110    continue
 
158
 120    s2 = delta*(co(1)*s1+si(1)*c1)
 
159
        c2 = delta*(co(1)*c1-si(1)*s1)
 
160
 130    ress(j) = s2
 
161
        resc(j) = c2
 
162
        do 140 i=1,4
 
163
          co(i) = co(i+1)
 
164
          si(i) = si(i+1)
 
165
 140    continue
 
166
 150  continue
 
167
c  calculate the integrals ress(j) and resc(j),j=n-6,n-5,n-4 by setting
 
168
c  up a divided difference table.
 
169
 160  do 190 j=1,3
 
170
        nmj = nm3-j
 
171
        i = 5-j
 
172
        call fpcsin(t(nm3),t(nmj),par,si(4),co(4),si(i-1),co(i-1),
 
173
     *  rs(j),rc(j))
 
174
        hs(i) = 0.
 
175
        hc(i) = 0.
 
176
        do 170 jj=1,j
 
177
          ipj = i+jj
 
178
          hc(ipj) = rc(jj)
 
179
          hs(ipj) = rs(jj)
 
180
 170    continue
 
181
        do 180 jj=1,3
 
182
          if(i.lt.jj) i = jj
 
183
          k = 5
 
184
          li = nmj
 
185
          do 180 ll=i,4
 
186
            lj = li+jj
 
187
            fac = t(lj)-t(li)
 
188
            hs(k) = (hs(k-1)-hs(k))/fac
 
189
            hc(k) = (hc(k-1)-hc(k))/fac
 
190
            k = k-1
 
191
            li = li+1
 
192
 180    continue
 
193
        ress(nmj) = hs(4)-hs(5)
 
194
        resc(nmj) = hc(4)-hc(5)
 
195
 190  continue
 
196
      return
 
197
      end