      DOUBLE PRECISION B(16),C(16),D(16),E(15),F(30),W,X,Y,Z
C    *************************************************************
      PRINT *, "*:*:* SYMMETRICAL ELLIPTIC FILTERS *:*:*"
C    *************************************************************
      DBN=.23025851;PI=3.1415927
C    *------------------------------------------------------------
   20 PRINT *, "STOPBAND EDGE (KHZ) ?";READ *, FS
      PRINT *, "PASSBAND EDGE (KHZ) ?";READ *, FP
      IF(ABS(FS-FP).LE.0)GO TO 20
   21 PRINT *, "NUMBER OF PEAKS (1-15) ?";READ *, N
      IF(N.LE.0)GO TO 20
      M=2*N+1
C    *
      FC=SQRT(FS*FP)
      R=FC+FC
      DO 24 K=1,2
      S=FS+FP
      DO 22 J=1,6
      P=SQRT(S*R)
      S=(S+R)/2
      IF(1.E08*(S-P).LT.S)GO TO 23
   22 R=P
   23 IF(K.GE.2)GO TO 25
      Q=M/S
   24 R=ABS(FS-FP)
   25 Q=Q*S
      S=EXP(-PI/Q)
      Y=DBLE(S)
      PRINT *, "CRITICAL Q",Q/(4*(1-S)*S**N)
C    *
      PRINT *, "STOPBAND REJECTION (DB) ?";READ *, S
      IF(S.LE.0)GO TO 20
      S=EXP(S*DBN/2)
      R=EXP(PI*Q)
      P=(ALOG(1+(S*S-1)/(R/4+1/R)**2))/DBN
      PRINT *, "PASSBAND RIPPLE (DB)",P
      R=R/(2*(S+SQRT(S*S-1)))
      R=ALOG(R+SQRT(R*R+1))/(2*Q)
      R=SIN(R)/COS(R)
      W=DBLE(R)
      PRINT *, "3 DB ABOUT (KHZ)",FP+(FS-FP)/(1+FC/(FP*R*R))
      PRINT *, "NOMINAL RESISTANCE (OHM) ?";READ *, R
      IF(R.LE.0)GO TO 21
C    *
      Z=Y;E(N)=W;W=W*W
      DO 26 J=1,M-1
   26 F(J)=1
      K=1
      DO 28 J=1,1024
      F(K)=F(K)*(1-Z)/(1+Z)
      IF(K.LT.M-1)GO TO 27
      Z=Z*Y
      X=((1-Z)/(1+Z))**2
      E(N)=E(N)*(W+X)/(1+W*X)
      K=0
   27 Z=Z*Y
      IF(Z.LT..25D-18)GO TO 29
   28 K=K+1
   29 DO 30 J=1,N
      F(J)=F(J)*F(M-J)
   30 F(M-J)=F(J)
C    *------------------------------------------------------------
      DO 31 J=1,N
      D(J)=F(2*J)*(1-F(J)**4)/F(J)
   31 B(J)=E(N)*F(J)
      C(1)=1/B(N)
      DO 32 J=1,N-1
      C(J+1)=(C(J)-B(N-J))/(1+C(J)*B(N-J))
   32 E(N-J)=E(N+1-J)+E(N)*D(J)/(1+B(J)*B(J))
C    *------------------------------------------------------------
      DO 40 J=1,N
      B(J)=((1+C(J)*C(J))*E(J)/D(J)-C(J)/F(J))/2
      C(J)=C(J)*F(J)
   40 D(J)=F(J)*F(J)
      B(N+1)=B(N);C(N+1)=C(N);D(N+1)=D(N)
C    *------------------------------------------------------------
      IF(N.EQ.1)GO TO 60
      DO 50 L=1,2
      DO 50 K=L+2,N+1,2
      DO 50 J=L,K-2,2
      Y=C(J)-C(K)
      Z=1/(Y/(B(J)*(D(K)-D(J)))-1)
      B(K)=(B(K)-B(J))*Z*Z-B(J)*(1+Z+Z)
   50 C(K)=Y*Z
C    *------------------------------------------------------------
      S=SNGL(B(N)/B(N+1)-1)
   60 Q=.0005/(PI*FC)
      P=Q*R;Q=Q/R
      IF(FS.LT.FP)GO TO 62
C    *
      PRINT *, "               *:* LOW-PASS FILTER *:*"
      DO 61 J=1,N
      C(J)=Q*C(J)
      D(J)=Q*B(J)*D(J)
      B(J)=P/B(J)
   61 F(J)=FC/F(J)
      C(N+1)=Q*C(N+1)
      GO TO 64
C    *
   62 PRINT *, "               *:* HIGH-PASS FILTER *:*"
      DO 63 J=1,N
      C(J)=Q/C(J)
      D(J)=Q/(B(J)*D(J))
      B(J)=P*B(J)
   63 F(J)=FC*F(J)
      C(N+1)=Q/C(N+1)
C    *
   64 PRINT *, "   KILOHERTZ              FARAD            HENRY"
      DO 65 J=1,N,2
      PRINT 71,C(J)
   65 PRINT 72,F(J),J,D(J),B(J)
      PRINT 71,C(N+1)
      IF(N.EQ.1)GO TO 70
      L=((N+1)/2)*2;K=M-1-L
      DO 66 J=L+2,M-1,2
      PRINT 72,F(K),K,D(K),B(K)
      PRINT 71,C(K)
   66 K=K-2
      PRINT *, "PRECISION TEST:",S
C    *------------------------------------------------------------
   70 PRINT *, "END (1) OR REPEAT (0) ?";READ *, N
      IF(N.NE.1)GO TO 20;STOP
   71 FORMAT(20X,D17.8)
   72 FORMAT(1X,D15.8,I4,2D17.8)
      STOP
      END
