1C/MEMBR ADD NAME=CHEBY,SSI=0
2      subroutine cheby(ordr,demi,ieo,dp,x0,tam,win)
3c!
4c sous-programme: cheby
5c fenetre de dolph chebyshev
6c en double precision
7c acheve le 05/12/85
8c ecrit par philippe touron
9c
10c
11c                  parametres entrants
12c                 -------------------
13c ordr, l'ordre du filtre (entier)
14c demi, l'ordre de la demi fenetre (entier)
15c ieo, indicateur de parite vaut 0 si ordr pair et 1 sinon (entier)
16c dp, attenuation en absolu (reelle)
17c x0, constante de la fenetre de chebyshev fonction de
18c      la largeur du lobe principale (reelle)
19c tam, tableau de travail (reels) qui doit etre dans un programme
20c    appelant dimensionne a 3 fois ordr)
21c
22c                  parametres sortants
23c                 -------------------
24c win, tableau de valeurs de la fenetre (reelles) qui doit dans
25c    un programme appelant etre dimensionne a ordr
26c
27c                  variables internes
28c                 ------------------
29c c0,c1,c2, des buffers de calcul (reels)
30c compt,indi,xcompt,xindi, compteurs de boucles (entiers et reels)
31c teta,gama, coefficients de la fenetres de chebyshev (reels)
32c
33c sous programmes appele: coshin
34c
35c!
36      double precision win(*),tam(*)
37      double precision x0,dp,xordr,teta,gama,freq
38      double precision pi,twopi,twn,c0,c1,c2,somm,xcompt,xindi
39      double precision coshin
40      integer ordr,demi,compt,indi,icompt,jcompt
41c
42      xordr=dble(ordr)
43      teta=(x0+1.)/2.0d+0
44      gama=(x0-1.)/2.0d+0
45      pi=acos(-1.0d+0)
46      twopi=2.0d+0*pi
47      c0=(xordr-1.0d+0)/2.0d+0
48      do 40 compt=1,ordr
49      icompt=ordr+compt
50      jcompt=icompt+ordr
51      xcompt=dble(compt)-1.0d+0
52      c1=xcompt/xordr
53      freq=teta*cos(twopi*c1)+gama
54      CRES=abs(freq)-1.0d+0
55      if (CRES .le. 0) then
56         goto 10
57      else
58         goto 20
59      endif
6010    tam(icompt)=dp*cos(c0*acos(freq))
61      goto 30
6220    tam(icompt)=dp*cosh(c0*coshin(freq))
6330    tam(jcompt)=0.0d+0
64c
65c modification si filtre d'ordr pair
66c
67      if(ieo.eq.1)goto 40
68      tam(jcompt)=-tam(icompt)*sin(pi*c1)
69      tam(icompt)= tam(icompt)*cos(pi*c1)
70      if(compt.gt.(ordr/2+1))then
71      tam(icompt)=-tam(icompt)
72      tam(jcompt)=-tam(jcompt)
73      endif
7440    continue
75c
76c transformee de fourrier inverse
77c
78      twn=twopi/xordr
79      do 60 compt=1,demi
80      xcompt=dble(compt)-1.0d+0
81      somm=0.0d+0
82      do 50 indi=1,ordr
83      icompt=indi+ordr
84      jcompt=icompt+ordr
85      xindi=dble(indi)-1.0d+0
86      somm=somm+tam(icompt)*cos(twn*xindi*xcompt)
87     & +tam(jcompt)*sin(twn*xindi*xcompt)
8850    continue
89      win(compt)=somm
9060    continue
91      c2=win(1)
92      do 70 compt=1,demi
93      win(compt)=win(compt)/c2
9470    continue
95      return
96      end
97