implicit real*8 (a-h,o-z)                                         
      real*4 x1,x2,x3,x4,x5                                             
      common/b/cn(30),cd(30),mn,md,const                                
      call reset                                                        
      write(6,1000)                                                     
1000  format( ' '/                                                      
     &' f3=0-> low or highpass. f1=passband cutoff. f2=stopband cutoff'/
     &' f1<f2 -> lowpass.'/                                             
     &' f3>0 -> bandpass. f1,f2 are limits of passband. f3 is limit of'/
     &' either high or low stopband._ we require f1<f2.'/               
     &' ripple=passband ripple in db. atten=stopband attenuation in db.'
     &)                                                                 
      write(6,910)                                                      
910   format( ' the filter coefficients will be found in a file called
     $ fort.7'/
     $' the listing of the program output will be found in fort.8'/
     $ ' enter the following arguments in floating point form.')        
      write(6,900)                                                      
900   format( ' enter sampling rate.')                                  
      read(5,901) sr                                                    
901   format(g10.0)                                                     
      xnyq=sr/2.d0                                                      
      write(6,903)                                                      
903   format( ' enter f1')                                              
      read( 5,901) x1                                                   
      write(6,905)                                                      
905   format( ' enter f2')                                              
      read( 5,901) x2                                                   
      write(6,906)                                                      
906   format( ' enter f3')                                              
      read( 5,901) x3                                                   
      write(6,907)                                                      
907   format( ' enter ripple')                                          
      read( 5,901) x4
      write(6,908)                                                      
908   format( ' enter attenuation in db')                               
      read( 5,901) x5                                                   
c      write(8,98) x1,x2,x3,x4,x5                                       
      write(6,97) x1,x2,x3,x4,x5   
c98    format( ' f1=',f8.3,' f2=',f8.3,' f3=',f8.3,' ripple=',          
c     $f8.3,' atten=',f8.3)                                             
97    format('c     f1=',f8.3,' f2=',f8.3,' f3=',f8.3,' ripple=', 
     $f8.3,' atten=',f8.3)
      f1=x1                                                             
      f2=x2                                                             
      f3=x3                                                             
      ripple=x4                                                         
      atten=x5                                                          
      call ellips(f1,f2,f3,ripple,atten,sr)                             
      call fresp(200,sr,0.d0,xnyq,f1)                                   
      m2=mn/2                                                           
      write(7,1234) m2,x1,x2,x3,x4,x5                                   
1234  format('/*elliptical filter with',i3,' sections.'/                  
     &'f1,f2,f3=',3g10.4,' ripple=',g10.4,' db=',g10.4,'*/')                 
      write(7,899) (cn(i),cd(i),i=1,mn)                                 
899   format('ell(',4(e15.8,','))
      write(7,898) const
898   format(e15.8,')')
      stop                                                              
      end                                                               
      subroutine reset                                                  
      implicit real*8 (a-h,o-z)                                         
      common/b/cn(30),cd(30),mn,md,const                                
      mn=0                                                              
      md=0                                                              
      write(8,200)                                                      
  200 format('//')                                                      
      do 100 m=1,30                                                     
      cn(m)=0.                                                          
      cd(m)=0.                                                          
100   continue                                                          
      return                                                            
      end                                                               
      subroutine ellips(f1,f2,f3,ripple,atten,samr)                     
c   designs an elliptic filter. all parameters real*8 .                 
c   f3=0 -> lowpass or highpass. f1=passband cutoff. f2=stopband cutoff.
c   f1<f2 -> lowpass.                                                   
c   f3>0 -> bandpass. f1,f2 are limits of passband. f3 is limit of      
c   either high or low stopband. we require f1<f2.                      
c   ripple=passband ripple in db. atten=stopband attenuation in db.     
c   samr=sampling rate in hz.                                           
c    after gold+rader; written by bilofsky, revised by steiglitz        
c    pp.61-65 (elliptic filters), 72,76 (mappings                       
c    from s-plane to z-plane), 87 (approximation                        
c    for u0 and evaluation of elliptic functions).                      
      implicit real*8 (a-h,o-z)                                         
      real*8 k,k1,kay,kprime,k1prim ,nn,kk,kkp,kk1,kk1p                 
      common/ellipt/k,kprime,cosp0,w1,hpass                             
      prime(dummy)=dsqrt(1.d0-dummy**2)                                 
      bpt(w)=dabs((cosp0-dcos(w))/dsin(w))                              
      pi=3.14159265358979d0                                             
      w1=2.d0*pi*f1/samr                                                
      w2=2.d0*pi*f2/samr                                                
      w3=2.d0*pi*f3/samr                                                
      hpass=0.d0                                                        
      cosp0=0.d0                                                        
      if(f3.gt.0.d0)goto1                                               
      if(f1.lt.f2)goto2                                                 
c  modify frequencies for high pass.                                    
      w1=pi-w1                                                          
      w2=pi-w2                                                          
      hpass=1.d0                                                        
c  compute analog frequencies for low/high pass                         
    2 w1=dtan(.5d0*w1)                                                  
      w2=dtan(.5d0*w2)                                                  
      goto3                                                             
c  compute analog frequencies for band pass.                            
    1 cosp0=dcos((w1+w2)/2.d0)/dcos((w1-w2)/2.d0)                       
      w1=bpt(w1)                                                        
      de=w3-w2                                                          
      if (de.lt.0.d0) de=w1-w3                                          
      w2=dmin1(bpt(w1-de),bpt(w2+de))                                   
c  compute params for poles,zeros in lambda plane                       
3     k=w1/w2                                                           
      kprime=prime(k)                                                   
      eps=dsqrt(10.d0**(.1d0*ripple)-1.d0)                              
      a=10.d0**(.05d0*atten)                                            
      k1=eps/dsqrt(a*a-1.d0)                                            
      k1prim =prime(k1)                                                 
      kk=kay(k)                                                         
      kk1=kay(k1)                                                       
      kkp=kay(kprime)                                                   
      kk1p=kay(k1prim )                                                 
      n=idint(kk1p*kk/(kk1*kkp))+1                                      
      nn=n                                                              
    5 u0=-kkp*dlog((1.d0+dsqrt(1.d0+eps*eps))/eps)/kk1p                 
c  now compute poles,zeros in lambda plane,                             
c    transform one by one to z plane.                                   
      dd=kk/nn                                                          
      tt=kk-dd                                                          
      dd=dd+dd                                                          
      n2=(n+1)/2                                                        
      do 4 i=1,n2                                                       
      if (i*2.gt.n) tt=0.d0                                             
      call stuff1(-kkp,tt,'zero')                                       
      call stuff1(u0,tt,'pole')                                         
4     tt=tt-dd                                                          
      return                                                            
      end                                                               
      subroutine stuff1(q,r,whatsi )                                    
c    transforms poles and zeros to z-plane; stuffs coeff. array         
      implicit real*8 (a-h,o-z)                                         
      real*8 k,kprime                                                   
      common/b/cn(30),cd(30),mn,md,const                                
      character*4 whatsi                                                
      complex*16 dcmplx,cdsqrt,dconjg,z,s                               
      common/ellipt/k,kprime,cosp0,w1,hpass                             
      call djelf(snr,cnr,dnr,r,kprime*kprime)                           
      call djelf(snqp,cnqp,dnqp,q,k*k)                                  
      omega=1-snqp*snqp*dnr*dnr                                         
      if ( omega .eq. 0.d0 ) omega=1.d-30                               
      sigma=w1*snqp*cnqp*cnr*dnr/omega                                  
      omega=w1*snr*dnqp/omega                                           
      s=dcmplx(sigma,omega)                                             
      j=1                                                               
      if (cosp0.eq.0.d0) goto 1                                         
      j=-1                                                              
    4 z=(-cosp0+dfloat(j)*cdsqrt(cosp0*cosp0+s*s-1.d0))/(s-1.d0)        
      go to 3                                                           
    1 z=(1.d0+s)/(1.d0-s)                                               
      if(hpass.ne.0.d0)z=-z                                             
    3 if(dabs(dimag(z)).le.10.d-10) goto 2                              
      if(dimag(z).lt.0.d0) z=dconjg(z)                                  
      if(whatsi.eq.'pole')goto5                                         
      mn=mn+1                                                           
      cn(mn)=-2.d0*dreal(z)                                             
      mn=mn+1                                                           
      cn(mn)=dreal(z)**2+dimag(z)**2                                    
      goto6                                                             
    5 md=md+1                                                           
      cd(md)=-2.d0*dreal(z)                                             
      md=md+1                                                           
      cd(md)=dreal(z)**2+dimag(z)**2                                    
    6 write(8,202)whatsi,z                                              
  202 format(' complex ',a4,' pair at ',d17.9,' +-j',d17.9)             
      if(j.gt.0.or.r.eq.0.d0)return                                     
      j=1                                                               
      go to 4                                                           
    2 x=dreal(z)                                                        
      if(whatsi.eq.'pole')goto7                                         
      mn=mn+1                                                           
      cn(mn)=-x                                                         
      mn=mn+1                                                           
      cn(mn)=0.d0                                                       
      goto8                                                             
    7 md=md+1                                                           
      cd(md)=-x                                                         
      md=md+1                                                           
      cd(md)=0.d0                                                       
    8 write(8,201)whatsi,x                                              
  201 format(' real ',a4,' at ',d17.9)                                  
      if(j.gt.0) return                                                 
      j=1                                                               
      go to 4                                                           
      end                                                               
      subroutine fresp(k,samr,f1,f2,f3)                                 
c    plots k pts. of freq. resp. from f1 to f2, norm. at f3             
      implicit real*8 (a-h,o-z)                                         
      complex*16 dcmplx,cdexp,tf,zm,zm2                                 
      common/b/cn(30),cd(30),mn,md,const                                
      pi=3.14159265358979d0                                             
      m2=mn/2                                                           
      write(8,200)m2,(cn(i),cd(i),i=1,mn)                               
  200 format('elliptic filter with ',i5,' sections'/4(d17.9))           
      w=pi*f3/(.5d0*samr)                                               
      zm=cdexp(dcmplx(0.d0,-1.d0*w))                                    
      zm2=zm*zm                                                         
      tf=(1.d0,0.d0)                                                    
      do 1 i=1,mn,2                                                     
    1 tf=tf*(1.d0+cn(i)*zm+cn(i+1)*zm2)/(1.d0+cd(i)*zm+cd(i+1)*zm2)     
      const=1.d0/cdabs(tf)                                              
      write(8,201)const                                                 
  201 format(' const=',d17.9)                                           
      write(8,205)                                                      
  205 format('/   freq     phase',10x,'    amp',10x,'    db.')          
      do 3 j=1,k                                                        
      freq=f1+(f2-f1)*dfloat(j-1)/dfloat(k-1)                           
      w=pi*freq/(.5d0*samr)                                             
      zm=cdexp(dcmplx(0.d0,-1.d0*w))                                    
      zm2=zm*zm                                                         
      tf=dcmplx(const,0.d0)                                             
      do 2 i=1,mn,2                                                     
    2 tf=tf*(1.d0+cn(i)*zm+cn(i+1)*zm2)/(1.d0+cd(i)*zm+cd(i+1)*zm2)     
      amp=cdabs(tf)                                                     
      if(amp.le.1.d-20)amp=1.d-20                                       
      x=dreal(tf)                                                       
      y=dimag(tf)                                                       
      phase=0.d0                                                        
      if(x.eq.0.d0 .and. y.eq.0.d0)goto4                                
      phase=(180.d0/pi)*datan2(y,x)                                     
    4 db=20.d0*dlog10(dmax1(amp,1.d-40))                                
    3 write(8,202)freq,phase,amp,db                                     
  202 format(' ',f10.2,2d17.9,f12.4)                                    
      return                                                            
      end                                                               
      double precision function kay(k)                                  
c    computes kay(k)=inverse sn(1)                                      
c    hastings, approx. for dig. comp., p. 172                           
      implicit real*8 (a-h,o-z)                                         
      double precision k,eta,peta,kk                                    
      dimension a(5),b(5)                                               
      data a/1.38629436112d0, .09666344259d0, .03590092383d0,           
     1    .03742563713d0, .01451196212d0/                               
      data b/.5d0, .12498593597d0, .06880248576d0, .03328355346d0,      
     1 .00441787012d0/                                                  
      kay=a(1)                                                          
      kk=b(1)                                                           
      eta=1.d0-k*k                                                      
      peta=eta                                                          
      do 1 i=2,5                                                        
      kay=kay+a(i)*peta                                                 
      kk=kk+b(i)*peta                                                   
1     peta=peta*eta                                                     
      kay=kay-kk*dlog(eta)                                              
      return                                                            
      end                                                               
      subroutine djelf(sn, cn, dn, x, sck)                              
c     ssp program: finds jacobian elliptic functions sn,cn,dn.          
      implicit real*8 (a-h,o-z)                                         
      dimension ari(12),geo(12)                                         
      double precision sn,cn,dn,x,sck,ari,geo,cm,y                      
     c,a,b,c,d                                                          
      cm=sck                                                            
      y=x                                                               
      if(sck)3,1,4                                                      
    1 d=dexp(x)                                                         
      a=1.d0/d                                                          
      b=a+d                                                             
      cn=2.d0/b                                                         
      dn=cn                                                             
      a=(d-a)/2.d0                                                      
      sn=a*cn                                                           
    2 return                                                            
    3 d=1.d0-sck                                                        
      cm=-sck/d                                                         
      d=dsqrt(d)                                                        
      y=d*x                                                             
    4 a=1.d0                                                            
      dn=1.d0                                                           
      do 6 i=1,12                                                       
      l=i                                                               
      ari(i)=a                                                          
      cm=dsqrt(cm)                                                      
      geo(i)=cm                                                         
      c=(a+cm)*.5d0                                                     
      if(dabs(a-cm)-1.d-9*a)7,7,5                                       
    5 cm=a*cm                                                           
    6 a=c                                                               
    7 y=c*y                                                             
      sn=dsin(y)                                                        
      cn=dcos(y)                                                        
      if(sn)8,13,8                                                      
    8 a=cn/sn                                                           
      c=a*c                                                             
      do 9 i=1,l                                                        
      k=l-i+1                                                           
      b=ari(k)                                                          
      a=c*a                                                             
      c=dn*c                                                            
      dn=(geo(k)+a)/(b+a)                                               
    9 a=c/b                                                             
      a=1.d0/dsqrt(c*c+1.d0)                                            
      if(sn)10,11,11                                                    
   10 sn=-a                                                             
      goto 12                                                           
   11 sn=a                                                              
   12 cn=c*sn                                                           
   13 if(sck)14,2,2                                                     
   14 a=dn                                                              
      dn=cn                                                             
      cn=a                                                              
      sn=sn/d                                                           
      return                                                            
      end