/* 
   this is an elliptic filter design program converted from
   basic to c-- in a klunky way.  No attempt has been made to
   get rid of GOTOs, etc.

   from Cuthbert's adaption of Amstutz's Algorithm

   I received the basic version in response to a posting on
   usenet.  The original basic is included as a comment at the
   end of the c-program.

   Converted:  November 3, 1994 by Steven D. Swift, P.E.
                                   Seattle, Washington
*/

#include <stdio.h>
#include <math.h>

main() 

{
	double B[16],C[16],D[16],E[16],F[16];
	double DN;
	double PI=3.14159265359;
	double sqrt();
	double log();
	double exp();
	double pow();
	double sin();
	double cos();
	double fabs();
	double FS, FP; /* stopband edge and passband edge */
	int N;  /* number of attenuation peaks */
	double  R, P, Q, S, Y, W, Z, X ;
	double FC; /* geometric mean of FS and FP */
	int K, J, M, L;

	printf("\nSymmetrical Elliptical Filter Synthesis\n");

	DN = log(10.0)/10.0;

G2010:	printf("Passband edge (MHz): ");
		scanf("%lf", &FP);
		FP = FP*1000.0; /* convert to kHz */
		printf( "Stopband edge (MHz): ");
		scanf("%lf", &FS);
		FS = FS*1000.0; /* convert to kHz */
		if ( fabs(FS-FP) <= 0.0 ) goto G2010;
G5555:	printf("\nNumber of peaks (1-15)= ");
	scanf("%d", &N);
	if (N<=0) goto G2010;
M = 2*N+1;
FC = sqrt(FS*FP);
R=FC+FC;
for ( K=1;K<=2;K++) {
	S=FS+FP;
		for( J= 1; J<=6; J++) {
			P=sqrt(S*R);
			S=(S+R)/2.0;
			if (1.0e+08*(S-P) <S) goto G2170;
			R=P; }
G2170: if(K>=2) goto G2200;
	Q=( (double) M)/S;
	R=fabs(FS-FP); }

G2200: Q=Q*S;
	S=exp(-PI/Q);
	Y=S;
	printf("Critical Q= %lf\n", Q/(4.0*(1.0-S)*pow(S,(double) N)) );

	printf("Stopband rejection (dB) = ");
	scanf("%lf", &S);
	if (S<=0.0) goto G2010;
	S=exp(S*DN/2.0);
	R=exp(PI*Q);

	P=(log(1.0+(S*S-1.0)/( pow((R/4.0+1.0/R),2.0) )) )/DN;

	printf( "Passband Ripple(dB)= %lf\n",P);
	R=R/(2.0*(S+sqrt(S*S-1.0)));
	R=log(R+sqrt(R*R+1.0))/(2.0*Q);
	R=sin(R)/cos(R);
	W=R;
	printf( "3 dB (kHz) about= %lf\n", FP+(FS-FP)/(1+FC/(FP*R*R)) );
	printf( "Nominal ohms resistance=");
	scanf("%lf", &R);
	if (R<=0.0) goto G5555;
	Z=Y;
	E[N]=W;
	W=W*W;

for (J=1; J<= M-1; J++) { 
	F[J]=1; 
		}

K=1;
for (J=1; J <= 1024; J++) {
	F[K]=F[K]*(1-Z)/(1+Z);
	if (K<M-1) goto G2500;
	Z=Z*Y;
	X= pow( ((1-Z)/(1+Z)), 2.0);
	E[N]=E[N]*(W+X)/(1+W*X);
	K=0;
G2500: Z=Z*Y;
	if( Z<2.5E-19) goto G2530;
	K=K+1;
	} 

G2530: for (J=1; J <= N; J++) {
	F[J]=F[J]*F[M-J];
	F[M-J]=F[J];
	} 

for (J=1;J <= N; J++) {
D[J]=F[2*J]*(1.0-pow(F[J],4.0))/F[J];
B[J]=E[N]*F[J];
}

C[1]=1/B[N];

for (J= 1; J <= N-1; J++) {
C[J+1]=(C[J]-B[N-J])/(1+C[J]*B[N-J]);
E[N-J]= E[N+1-J] + E[N]*D[J]/(1+B[J]*B[J]); }

for (J= 1; J <= N; J++) {
B[J]=((1+C[J]*C[J])*E[J]/D[J]-C[J]/F[J])/2.0;
C[J]=C[J]*F[J];
D[J]=F[J]*F[J]; }

B[N+1]=B[N];
C[N+1]=C[N];
D[N+1]=D[N];

if (N==1) goto G6020;

L=1;
G5030: for(K=L+2; K <= N+1; K+=2) {
	 for(J=L; J <= K-2; J+=2) {
		 Y=C[J]-C[K];
		 Z=1.0/(Y/(B[J]*(D[K]-D[J]))-1.0);
		 B[K]=(B[K]-B[J])*Z*Z-B[J]*(1.0+Z+Z);
		 C[K]=Y*Z; }  }

if (L==2) goto G6010;

L=2;
goto G5030;

G6010:	S=B[N]/B[N+1]-1;

G6020:	Q=.0005/(PI*FC);

P=Q*R;
Q=Q/R;

if (FS<FP) goto G6150;
printf("\n\t** Low-Pass Filter **\n\n");

for (J=1; J<=N; J++) {
	C[J]=Q*C[J];
	D[J]=Q*B[J]*D[J];
	B[J]=P/B[J];
	F[J]=FC/F[J]; }

C[N+1]=Q*C[N+1];
goto G6230;

G6150: printf("\n\t** High-Pass Filter **\n\n");

for (J=1;J<=N;J++) {
	C[J]=Q/C[J];
	D[J]=Q/(B[J]*D[J]);
	B[J]=P*B[J];
	F[J]=FC*F[J]; }

C[N+1]=Q/C[N+1];

G6230: printf("kHz \t\tpicoFarads \t\tmicroHenrys\n\n");

for (J=1;J<=N;J+=2) {
printf("\t\t%0.1lf \n\n", C[J]/1.0e-12);
printf("%0.3lf   %d\t%0.1lf\t\t%0.3lf\n",F[J],J,D[J]/1.0e-12,B[J]/1.0e-6 );
}

printf("\t\t%0.1lf \n\n",  C[N+1]/1.0e-12);

if (N==1) exit() ;

L=(int) ((N+1)/2)*2;
K=M-1-L;

for (J=L+2; J<=M-1; J+=2) {
printf("%0.3lf   %d \t%0.1lf \t\t%0.3lf\n",F[K],K,D[K]/1.0e-12, B[K]/1.0e-6);
printf("\t\t%0.1lf \n\n",  C[K]/1.0e-12);
K=K-2; }

/* printf("\nprecision test: %lf\n",S); */

}

/* --------------------------------------------------------------- */

/*  Original basic program is below

11 ' from Cuthbert's adaption of Amstutz's Algorithm
1010 DIM B(16),C(16),D(16),E(16),F(16)
1015 CLS
1018 PRINT"Symmetrical Elliptical Filter Synthesis"
1020 DN=LOG(10)/10:PI=3.1415926#
2010 PRINT "Stopband edge (KHz): ":INPUT FS
2020 PRINT "Passband edge (KHz): ":INPUT FP
2030 IF ABS(FS-FP)<=0 GOTO 2010
2040 PRINT"Number of peaks (1-15)= ";:INPUT N
2050 IF N<=0 GOTO 2010
2060 M=2*N+1
2080 FC=SQR(FS*FP)
2090 R=FC+FC
2100 FOR K= 1 TO 2
2110 S=FS+FP
2120 FOR J= 1 TO 6
2130 P=SQR(S*R)
2140 S=(S+R)/2
2150 IF 1E+08 * (S-P) <S GOTO 2170
2160 R=P:NEXT J
2170 IF K>=2 GOTO 2200
2180 Q=M/S
2190 R=ABS(FS-FP): NEXT K
2200 Q=Q*S
2210 S=EXP(-PI/Q)
2220 Y=S
2230 PRINT"Critical Q=";Q/(4*(1-S)*S^N)
2250 PRINT"Stopband rejection (dB) = ";:INPUT S
2260 IF S<=0 GOTO 2010
2270 S=EXP(S*DN/2)
2280 R=EXP(PI*Q)
2290 P=(LOG(1+(S*S-1)/(R/4+1/R)^2))/DN
2300 PRINT "Passband Ripple(dB)=";P
2310 R=R/(2*(S+SQR(S*S-1)))
2320 R=LOG(R+SQR(R*R+1))/(2*Q)
2330 R=SIN(R)/COS(R)
2340 W=R
2350 PRINT "3 dB (KHz) about=";FP+(FS-FP)/(1+FC/(FP*R*R))
2360 PRINT "Nominal ohms resistance=";:INPUT R
2370 IF R<=0 GOTO 2040
2390 Z=Y:E(N)=W:W=W*W
2400 FOR J= 1 TO M-1
2410 F(J)=1: NEXT J
2420 K=1
2430 FOR J = 1 TO 1024
2440 F(K)=F(K)*(1-Z)/(1+Z)
2450 IF K<M-1 GOTO 2500
2460 Z=Z*Y
2470 X=((1-Z)/(1+Z))^2
2480 E(N)=E(N)*(W+X)/(1+W*X)
2490 K=0
2500 Z=Z*Y
2510 IF Z<2.5E-19 GOTO 2530
2520 K=K+1: NEXT J
2530 FOR J= 1 TO N
2540 F(J)=F(J)*F(M-J)
2550 F(M-J)=F(J):NEXT J
3010 FOR J= 1 TO N
3020 D(J)=F(2*J)*(1-F(J)^4)/F(J)
3030 B(J)=E(N)*F(J):NEXT J
3040 C(1)=1/B(N)
3050 FOR J= 1 TO N-1
3060 C(J+1)=(C(J)-B(N-J))/(1+C(J)*B(N-J))
3070 E(N-J)=E(N+1-J)+E(N)*D(J)/(1+B(J)*B(J)):NEXT J
4010 FOR J= 1 TO N
4020 B(J)=((1+C(J)*C(J))*E(J)/D(J)-C(J)/F(J))/2
4030 C(J)=C(J)*F(J)
4040 D(J)=F(J)*F(J):NEXT J
4050 B(N+1)=B(N):C(N+1)=C(N):D(N+1)=D(N)
5010 IF N=1 GOTO 6020
5020 L=1
5030 FOR K= L+2 TO N+1 STEP 2
5040 FOR J= L TO K-2 STEP 2
5050 Y=C(J)-C(K)
5060 Z=1/(Y/(B(J)*(D(K)-D(J)))-1)
5070 B(K)=(B(K)-B(J))*Z*Z-B(J)*(1+Z+Z)
5080 C(K)=Y*Z:NEXT J
5081 NEXT K
5082 IF L= 2 GOTO 6010
5083 L=2:GOTO 5030
6010 S=B(N)/B(N+1)-1
6020 Q=.0005/(PI*FC)
6030 P=Q*R:Q=Q/R
6040 IF FS<FP GOTO 6150
6060 PRINT"      ** Low-Pass Filter **"
6070 FOR J= 1 TO N
6080 C(J)=Q*C(J)
6090 D(J)=Q*B(J)*D(J)
6100 B(J)=P/B(J)
6110 F(J)=FC/F(J):NEXT J
6120 C(N+1)=Q*C(N+1)
6130 GOTO 6230
6150 PRINT"      ** High-Pass Filter **"
6160 FOR J= 1 TO N
6170 C(J)=Q/C(J)
6180 D(J)=Q/(B(J)*D(J))
6190 B(J)=P*B(J)
6200 F(J)=FC*F(J): NEXT J
6210 C(N+1)=Q/C(N+1)
6230 PRINT "      KHz                          microFarads         microHenrys"
6240 FOR J= 1 TO N STEP 2
6250 PRINT TAB(30);:PRINT USING "#########.#####";C(J)/.000001
6260 PRINT USING"#########.#####"; F(J);J;D(J)/.000001;B(J)/.000001:NEXT J
6270 PRINT TAB(30);:PRINT USING "#########.#####";C(N+1)/.000001
6280 IF N=1 THEN STOP
6290 L=(INT((N+1)/2))*2:K=M-1-L
6300 FOR J= L+2 TO M-1 STEP 2
6310 PRINT USING "#########.#####";F(K);K;D(K)/.000001;B(K)/.000001
6320 PRINT TAB(30);:PRINT USING "#########.#####";C(K)/.000001
6330 K=K-2:NEXT J
6340 PRINT:PRINT:PRINT" precision test:";S
7020 INPUT Q$
7030 IF Q$="q" OR Q$="Q" THEN RUN "menu" ELSE GOTO 1015

*/