0010**************************************************************
0020 PRINT:"** ANTIMETRICAL ELLIPTIC FILTERS * TYPES A,B AND C **"
1000**************************************************************
1010 CHARACTER *1T
1020 DOUBLE PRECISION U,V,W,X,Y,Z,A0,E0,E8,FP,FS,SR,TQ,T0,
1030&B(16),C(16),D(16),E(30),F(16),R(15),S(15),DB(16),TB(16)
1040 DBN=.23025851;PI=3.1415927
1050 10 PRINT:;PRINT:;PRINT:
1060 PRINT:"REJECTION,RIPPLE(DB),HALF-DEGREE(2-15),TYPE(A,B OR C)"
1070 READ:AS,AP,M,T
1080 IF(AS.LE.AP)THEN STOP;N=2*M
2000*=============================================================
2010 EAS=EXP(DBN*AS)-1
2020 EAP=EXP(DBN*AP)-1
2030 V=SQRT(EAS/EAP)+SQRT(EAS/EAP-1)
2040 U=PI*PI/(2*ALOG(V+V))
2050 V=V/(SQRT(EAS)+SQRT(EAS+1))
2060 W=U*ALOG(V+SQRT(V*V+1))/PI
2070 W=SIN(W)/COS(W);A0=W;W=W*W
2080 Y=DEXP(-U);Z=Y;K=M-1
2090 DO 20 J=1,N
2100 20 E(J)=1
2110 DO 22 J=1,1024
2120 IF(K.NE.M)GO TO 21
2130 X=((1-Z)/(1+Z))**2
2140 A0=A0*(W+X)/(1+W*X)
2150 21 E(K)=E(K)*(1-Z)/(1+Z)
2160 Z=Z*Y;IF(Z.LT..25E-18)GO TO 23;K=K-1
2170 22 IF(K.EQ.0)K=N
2180 23 E(M)=0;E(N)=E(N)*E(N)
2190 PRINT:;PRINT:"U =",U,"  A0=",A0;PRINT:"EP=",E(N)
2200 DO 24 J=1,M-1
2210 E(J)=-E(J)*E(N-J)
2220 PRINT:"E =",-E(J)
2230 24 E(N-J)=-E(J)
2000*-------------------------------------------------------------
2250 X=DSQRT(A0*A0+1/(A0*A0)+E(N)*E(N)+1/(E(N)*E(N)))
2260 DO 25 J=1,M-1,2;K=(J+1)/2
2270 Y=A0*E(J);Y=Y+1/Y
2280 Z=E(N)*E(J)
2290 R(K)=E(M-J)*(1/Z-Z)/Y;S(K)=-X/Y
2300 PRINT:"RE=",R(K),"  SE=",S(K)
2310 R(M-K+1)=R(K)
2320 25 S(M-K+1)=-S(K)
2330 IF(K+K.EQ.M)GO TO 30
2340 R(K+1)=-A0;S(K+1)=0
2350 PRINT:"RE=",-A0
3000*=============================================================
3010 30 PRINT:;IT=2;IF(T.EQ."A")THEN IT=1
3020 E8=-E(1);IF(T.EQ."A")E8=E(N)
3030 EO=E(N);IF(T.EQ."C")EO=-E(1)
3040 FP=DSQRT((E(N)+EO)/(1+E(N)*E8))
3050 FS=DSQRT((1+E(N)*EO)/(E(N)+E8))
3060 D(1)=0;DO 31 J=IT,M
3070 D(J)=(E(2*J-1)+E8)/(1+E(2*J-1)*EO)
3080 31 F(J)=SQRT(1/D(J))
3090*-------------------------------------------------------------
3100 SR=0;TQ=0;T0=0;B(1)=0;I=1
3110 DO 32 J=1,M
3120 W=(A0**2+E(2*J-1)**2)/(1+(A0*E(2*J-1))**2)
3130 X=(1+EO*E8)*S(J)+EO+E8*W
3140 Y=EO**2+2*EO*S(J)+W
3150 Z=1+2*E8*S(J)+E8**2*W
3160 U=DSQRT(Y/Z);V=X/Z
3170 R(J)=DSQRT((U-V)/2);S(J)=DSQRT((U+V)/2)
3180 PRINT:"RF=",-R(J)/FP,"  SF=",S(J)/FP
3190*
3200 SR=SR+R(J)/U
3210 I=-I;W=I*R(J)/S(J)
3220 TQ=(TQ+W)/(1-TQ*W)
3230 IF(T.NE."A")GO TO 32
3240 U=(F(2)-S(J))/R(J);V=(F(2)+S(J))/R(J)
3250 W=I*(V-U)/(1+U*V)
3260 T0=(T0+W)/(1-T0*W)
3270 32 B(1)=B(1)+R(J)
4000*=============================================================
4010 IF(T.EQ."A")T0=T0/(1+DSQRT(1+T0*T0))
4020 DO 40 K=IT,M
4030 DB(K)=0;TB(K)=T0;I=1
4040 DO 40 J=1,M
4050 DB(K)=DB(K)+1/(R(J)+(F(K)-S(J))**2/R(J))
4060&           +1/(R(J)+(F(K)+S(J))**2/R(J))
4070 I=-I;W=(F(K)-I*S(J))/R(J)
4080 40 TB(K)=(TB(K)+W)/(1-TB(K)*W)
5000*=============================================================
5010 D(M+1)=D(M);F(M+1)=F(M);DB(M+1)=DB(M);TB(M+1)=TB(M);C(1)=0
5020 DO 50 J=1,M+1-IT,2
5030 50 TB(M+1-J)=-1/TB(M+1-J)
5040 DO 51 J=IT,M+1
5050 B(J)=(1+TB(J)**2)*DB(J)/(4*D(J))-TB(J)*F(J)/2
5060 51 C(J)=TB(J)/F(J)
6000*=============================================================
6010 DO 60 L=1,2
6020 DO 60 K=L+2,M+1,2
6030 DO 60 J=L,K-2,2
6040 U=C(J)-C(K)
6050 V=1/(U/(B(J)*(D(K)-D(J)))-1)
6060 B(K)=(B(K)-B(J))*V*V-B(J)*(V+V+1)
6070 60 C(K)=U*V
7000*=============================================================
7010 W=1;IF(T.NE."C")W=((1-TQ*T0)/(TQ+T0))**2
7020 DO 70 J=1,M+1,2;B(J)=B(J)*W
7030 70 C(J)=C(J)*W;PRINT:
7040 PRINT:"LOAD RESISTANCE=",SNGL(W),"   (",SNGL(1/W),"  )"
7050 PRINT:;PRINT:"                     L  ( C )        C  ( L )
7060&          PEAK"
7070 IF(T.NE."A")
7080&PRINT:"               1                ",SNGL(FP/B(1))
7090 V=0;DO 71 J=IT,M;V=V+C(J)
7100 71 PRINT:J,SNGL(FP*C(J)),SNGL(FP/B(J)),SNGL(F(J)/FP)
7110 PRINT:M+1;SNGL(FP*C(M+1)),"  STOPBAND EDGE=",SNGL(FS/FP)
7120 PRINT:
7130 PRINT:"TESTS",SNGL(B(M)/B(M+1)-1),SNGL((W+1)*SR-V-C(M+1))
7140 GO TO 10
