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

/* structure which saves line parameters */

typedef struct {

	float duty_cycle; /* duty cycle of pulse as fraction of T */

	} DUTY_CYCLE;

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

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

			biphase()

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

Function computes the ABCD matrix of a transfer function corresponding
to a biphase pulse with adjustable duty cycle

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 biphase(current,f)
TWOPORT *current; 	/* twoport being modeled */
double f;		/* frequency */

{
	ABCD tmsect;
	POINTER calloc();
	DUTY_CYCLE *param;
	static COMPLEX tf;
	double alpha,temp1,temp2;


/* 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(DUTY_CYCLE))) == NULL)
			error_alloc();

		param = (DUTY_CYCLE *) current->param;

		if(sscanf(current->line,"biphase%f", &param->duty_cycle) != 1) {
			fprintf(stderr,
		"biphase(): cannot read parameters on line\n%s\n",
			current->line);
			exit(1);
			}

		fprintf(fdata,"Biphase pulse with duty cycle %f T\n", param->duty_cycle);
		}

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

	else 	{

		param = (DUTY_CYCLE *) current->param;

		/* Calculate transfer function of pulse */

		alpha = param->duty_cycle;
		temp1 = PI * f * T * alpha;
		if(f == 0) temp2 = T * alpha ;
		else	{
			temp2 = T * alpha * sin(temp1) / temp1;
			}

		/* correct for delay and and negative part */
		tf.real = temp2 * (cos(temp1) - cos(3.*temp1));
		tf.imag = - temp2 * (sin(temp1) - sin(3.*temp1));


	/* 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);
}
