← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/interpolate/fitpack/fourco.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/interpolate/fitpack/fourco.f
 
1
      subroutine fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)
 
2
c  subroutine fourco calculates the integrals
 
3
c                    /t(n-3)
 
4
c    ress(i) =      !        s(x)*sin(alfa(i)*x) dx    and
 
5
c              t(4)/
 
6
c                    /t(n-3)
 
7
c    resc(i) =      !        s(x)*cos(alfa(i)*x) dx, i=1,...,m,
 
8
c              t(4)/
 
9
c  where s(x) denotes a cubic spline which is given in its
 
10
c  b-spline representation.
 
11
c
 
12
c  calling sequence:
 
13
c     call fourco(t,n,c,alfa,m,ress,resc,wrk1,wrk2,ier)
 
14
c
 
15
c  input parameters:
 
16
c    t    : real array,length n, containing the knots of s(x).
 
17
c    n    : integer, containing the total number of knots. n>=10.
 
18
c    c    : real array,length n, containing the b-spline coefficients.
 
19
c    alfa : real array,length m, containing the parameters alfa(i).
 
20
c    m    : integer, specifying the number of integrals to be computed.
 
21
c    wrk1 : real array,length n. used as working space
 
22
c    wrk2 : real array,length n. used as working space
 
23
c
 
24
c  output parameters:
 
25
c    ress : real array,length m, containing the integrals ress(i).
 
26
c    resc : real array,length m, containing the integrals resc(i).
 
27
c    ier  : error flag:
 
28
c      ier=0 : normal return.
 
29
c      ier=10: invalid input data (see restrictions).
 
30
c
 
31
c  restrictions:
 
32
c    n >= 10
 
33
c    t(4) < t(5) < ... < t(n-4) < t(n-3).
 
34
c    t(1) <= t(2) <= t(3) <= t(4).
 
35
c    t(n-3) <= t(n-2) <= t(n-1) <= t(n).
 
36
c
 
37
c  other subroutines required: fpbfou,fpcsin
 
38
c
 
39
c  references :
 
40
c    dierckx p. : calculation of fouriercoefficients of discrete
 
41
c                 functions using cubic splines. j. computational
 
42
c                 and applied mathematics 3 (1977) 207-209.
 
43
c    dierckx p. : curve and surface fitting with splines, monographs on
 
44
c                 numerical analysis, oxford university press, 1993.
 
45
c
 
46
c  author :
 
47
c    p.dierckx
 
48
c    dept. computer science, k.u.leuven
 
49
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
 
50
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
 
51
c
 
52
c  latest update : march 1987
 
53
c
 
54
c  ..scalar arguments..
 
55
      integer n,m,ier
 
56
c  ..array arguments..
 
57
      real*8 t(n),c(n),wrk1(n),wrk2(n),alfa(m),ress(m),resc(m)
 
58
c  ..local scalars..
 
59
      integer i,j,n4
 
60
      real*8 rs,rc
 
61
c  ..
 
62
      n4 = n-4
 
63
c  before starting computations a data check is made. in the input data
 
64
c  are invalid, control is immediately repassed to the calling program.
 
65
      ier = 10
 
66
      if(n.lt.10) go to 50
 
67
      j = n
 
68
      do 10 i=1,3
 
69
        if(t(i).gt.t(i+1)) go to 50
 
70
        if(t(j).lt.t(j-1)) go to 50
 
71
        j = j-1
 
72
  10  continue
 
73
      do 20 i=4,n4
 
74
        if(t(i).ge.t(i+1)) go to 50
 
75
  20  continue
 
76
      ier = 0
 
77
c  main loop for the different alfa(i).
 
78
      do 40 i=1,m
 
79
c  calculate the integrals
 
80
c    wrk1(j) = integral(nj,4(x)*sin(alfa*x))    and
 
81
c    wrk2(j) = integral(nj,4(x)*cos(alfa*x)),  j=1,2,...,n-4,
 
82
c  where nj,4(x) denotes the normalised cubic b-spline defined on the
 
83
c  knots t(j),t(j+1),...,t(j+4).
 
84
         call fpbfou(t,n,alfa(i),wrk1,wrk2)
 
85
c  calculate the integrals ress(i) and resc(i).
 
86
         rs = 0.
 
87
         rc = 0.
 
88
         do 30 j=1,n4
 
89
            rs = rs+c(j)*wrk1(j)
 
90
            rc = rc+c(j)*wrk2(j)
 
91
  30     continue
 
92
         ress(i) = rs
 
93
         resc(i) = rc
 
94
  40  continue
 
95
  50  return
 
96
      end