← Back to branch summary

~daniele-bigoni/pyorthpol/pyorthpol-python3

« back to all changes in this revision

Viewing changes to PyORTHPOL/ORTHPOLxx/src/dmcheb.f

  • Committer: Daniele Bigoni
  • Date: 2014-12-16 10:43:03 UTC
  • Revision ID: dabi@dtu.dk-20141216104303-ba4utrt7ln5g39lq
started with launchpad

Show diffs side-by-side

added added

removed removed

Lines of Context:
PyORTHPOL/ORTHPOLxx/src/dmcheb.f
 
1
c
 
2
c
 
3
      subroutine dmcheb(n,ncapm,mcd,mp,dxp,dyp,dquad,deps,iq,
 
4
     *idelta,finld,finrd,dendl,dendr,dxfer,dwfer,da,db,dnu,dalpha,
 
5
     *dbeta,ncap,kount,ierrd,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2)
 
6
c
 
7
c This is a double-precision version of the routine  mccheb.
 
8
c
 
9
      double precision dxp,dyp,deps,dendl,dendr,dxfer,dwfer,da,db,
 
10
     *dnu,dalpha,dbeta,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2,dsum,dp1,
 
11
     *dp,dpm1
 
12
      dimension dxp(*),dyp(*),dendl(mcd),dendr(mcd),dxfer(ncapm),
 
13
     *dwfer(ncapm),da(*),db(*),dnu(*),dalpha(n),dbeta(n),dbe(n),
 
14
     *dx(ncapm),dw(ncapm),dxm(*),dwm(*),ds(n),ds0(*),ds1(*),ds2(*)
 
15
      logical finld,finrd
 
16
c
 
17
c The arrays  dxp,dyp  are assumed to have dimension  mp  if mp > 0,
 
18
c the arrays  da,db  dimension 2*n-1, the arrays  dnu,ds0,ds1,ds2
 
19
c dimension  2*n, and the arrays  dxm,dwm  dimension  mc*ncapm+mp.
 
20
c
 
21
      nd=2*n
 
22
      if(idelta.le.0) idelta=1
 
23
      if(n.lt.1) then
 
24
        ierrd=-1
 
25
        return
 
26
      end if
 
27
      incap=1
 
28
      kount=-1
 
29
      ierrd=0
 
30
      do 10 k=1,n
 
31
        dbeta(k)=0.d0
 
32
   10 continue
 
33
      ncap=(nd-1)/idelta
 
34
   20 do 30 k=1,n
 
35
        dbe(k)=dbeta(k)
 
36
   30 continue
 
37
      kount=kount+1
 
38
      if(kount.gt.1) incap=2**(kount/5)*n
 
39
      ncap=ncap+incap
 
40
      if(ncap.gt.ncapm) then
 
41
        ierrd=ncapm
 
42
        return
 
43
      end if
 
44
      mtncap=mcd*ncap
 
45
      do 50 i=1,mcd
 
46
        im1tn=(i-1)*ncap
 
47
        if(iq.eq.1) then
 
48
          call dquad(ncap,dx,dw,i,ierr)
 
49
        else
 
50
          call dqgp(ncap,dx,dw,i,ierr,mcd,finld,finrd,dendl,dendr,
 
51
     *      dxfer,dwfer)
 
52
        end if
 
53
        if(ierr.ne.0) then
 
54
          ierrd=i
 
55
          return
 
56
        end if
 
57
        do 40 k=1,ncap
 
58
          dxm(im1tn+k)=dx(k)
 
59
          dwm(im1tn+k)=dw(k)
 
60
   40   continue
 
61
   50 continue
 
62
      if(mp.ne.0) then
 
63
        do 60 k=1,mp
 
64
          dxm(mtncap+k)=dxp(k)
 
65
          dwm(mtncap+k)=dyp(k)
 
66
   60   continue
 
67
      end if
 
68
      mtnpmp=mtncap+mp
 
69
      do 90 k=1,nd
 
70
        km1=k-1
 
71
        dsum=0.d0
 
72
        do 80 i=1,mtnpmp
 
73
          dp1=0.d0
 
74
          dp=1.d0
 
75
          if(k.gt.1) then
 
76
            do 70 l=1,km1
 
77
              dpm1=dp1
 
78
              dp1=dp
 
79
              dp=(dxm(i)-da(l))*dp1-db(l)*dpm1
 
80
   70       continue
 
81
          end if
 
82
          dsum=dsum+dwm(i)*dp
 
83
   80   continue
 
84
        dnu(k)=dsum
 
85
   90 continue
 
86
      call dcheb(n,da,db,dnu,dalpha,dbeta,ds,ierr,ds0,ds1,ds2)
 
87
      do 100 k=1,n
 
88
        if(dabs(dbeta(k)-dbe(k)).gt.deps*dabs(dbeta(k))) goto 20
 
89
  100 continue
 
90
      return
 
91
      end
 
92