#include <stdio.h>
#include <math.h>
#include "pi.h"
#include "type.h"

/* structure which saves line parameters */

typedef struct {

	int poles;	/* number of poles */
	float cutoff_f; /* cutoff frequency as a fraction of 1/T */

	} BUTTERWORTH;

extern double T;
extern FILE *fdata;
extern COMPLEX unity,zero;

/**********************************************************************

			butterworth()

************************************************************************

Function computes the ABCD matrix of a transfer function corresponding
to a butterworth lowpass filter with given number of poles and
cutoff frequency

The twoport does not affect the impedance, ie reproduces at its input
the impedance it sees at its output

Inputs: pointer to the data structure for the twoport and the frequency
Output: ABCD matrix for the twoport
*/

ABCD butterworth(current,f)
TWOPORT *current; 	/* twoport being modeled */
double f;		/* frequency */

{
	ABCD tmsect;
	POINTER calloc();
	BUTTERWORTH *param;
	static COMPLEX tf;
	float omeganorm;
	COMPLEX fac1,fac2,fac3;


/* the first time this function is called we simply read the parameters
   of the line from the input line */

	if(current->line != NULL) {

		/* allocate storage for the line parameters -- this
		   function may be called for many lines so it has to
		   be unique */

		if((current->param = calloc(1,sizeof(BUTTERWORTH))) == NULL)
			error_alloc();

		param = (BUTTERWORTH *) current->param;

		if(sscanf(current->line,"butterworth %d %f", &param->poles,
			&param->cutoff_f) != 2) {
			fprintf(stderr,
		"butterworth(): cannot read parameters on line\n%s\n",
			current->line);
			exit(1);
			}
		if(param->poles > 5 || param->poles < 1) {
			fprintf(stderr,"butterworth(): five poles maximum!");
			exit(2);
			}

		fprintf(fdata,"Butterworth filter with %d poles and\
			cutoff %f/T\n", param->poles,param->cutoff_f);
		}

/* every other time this function is called we calculate the ABCD matrix */

	else 	{

		param = (BUTTERWORTH *) current->param;

		/* Calculate transfer function of butterworth */

		/* first calculate the factors in the rational denominator */
		omeganorm = T * f / param->cutoff_f;
		switch(param->poles)  {
			case 1:
				fac1.real = 1.;
				fac1.imag = omeganorm;
				fac2 = fac3 = unity;
				break;
			case 2:
				fac1.real = 1. - omeganorm*omeganorm;
				fac1.imag = 1.4142 * omeganorm;
				fac2 = fac3 = unity;
				break;
			case 3:
				fac1.real = 1.;
				fac1.imag = omeganorm;
				fac2.real = 1. - omeganorm*omeganorm;
				fac2.imag = omeganorm;
				fac3 = unity;
				break;
			case 4:
				fac1.real = 1. - omeganorm*omeganorm;
				fac1.imag = .7654 * omeganorm;
				fac2.real = 1. - omeganorm*omeganorm;
				fac2.imag = 1.8478 * omeganorm;
				fac3 = unity;
				break;
			case 5:
				fac1.real = 1.;
				fac1.imag = omeganorm;
				fac2.real = 1. - omeganorm*omeganorm;
				fac2.imag = .6180 * omeganorm;
				fac3.real = 1. - omeganorm*omeganorm;
				fac3.imag = 1.6180 * omeganorm;
				break;
		}

		tf = cdiv(unity,fac1);
		tf = cdiv(tf,fac2);
		tf = cdiv(tf,fac3);


	/* ABCD matrix has 1/tf in A and D, zero otherwise --
	   this insures that the input impedance equals output impedance */

		/* if the transfer function is zero, 1/tf should be
		   make large to avoid float exception */
		if(tf.real*tf.real < 1.e-12 && tf.imag*tf.imag < 1.e-12)
			tmsect.a.real = tmsect.a.imag = 1.e6;
		else
			tmsect.a = cdiv(unity,tf);
		tmsect.d = tmsect.a;
		tmsect.b = tmsect.c = zero;
		}

	return(tmsect);
}
