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

		COMPUTATION FOR LINE MODEL

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

This file contains the computations of the line model,
including the ABCD matrix for the line and the FFT for
determining the impulse response
*/

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

/* declare external variables from linemod.c */

extern int mexp,nexp;
extern double T;
extern FILE *fimped,*ftime,*ffreq;

/* variables used by all the routines in this file */

ABCD *tm;	/* array of ABCD matricies, one for each frequency */

int n,m,nm;

COMPLEX unity = {1.0,0.0};
COMPLEX zero = {0.0,0.0};

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

			abcd()

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

Function calculates the ABCD matrix for the line model at
each frequency

The frequencies are chosen to be appropriate for the subsequent
FFT calculation of the impulse response
*/

abcd(first)
TWOPORT *first;	/* points to first twoport in line to be characterized
		   by ABCD parameters */
{
	POINTER calloc();
	double f,f0;	/* frequency */
	int i;

/* calculate 2**mexp and 2**nexp  */

	m=n=1;
	for(i=0; i<mexp; ++i) m *= 2;
	for(i=0; i<nexp; ++i) n *= 2;

	nm = n*m;	/* no. points in FFT */

/* allocate memory for array of ABCD matricies, one at each freq. */

	tm = (ABCD *) calloc(nm,sizeof(ABCD));

/* calculate the ABCD matrix at each frequency */

	f0 = 1.0/(m*T);
	for(i=0; i<nm; ++i) {

		f = i*f0;    /* frequency corresponding to matrix */

		/* routine abcdf() calculates the ABCD matrix at
		   a single frequency */
		abcdf(first,f,tm+i);
		}
	return;
}

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

			abcdf()

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

Function computes the ABCD matrix at a single frequency
*/

abcdf(first,f,tf)
TWOPORT *first;	/* first twoport in the transmission line */
double f;	/* frequency */
ABCD *tf;	/* points to ABCD matrix that is to be calculated */

{
	ABCD tmsect;  /* ABCD matrix for twoport pointed to by first */
	ABCD tmbridge;  /* ABCD matrix for a bridge tap */

/* calculate the ABCD matrix for the first twoport in the line */

	/* first case is where first twoport is not a bridge tap */
	if(first->tap == NULL)
		tmsect = (*first->function)(first,f);

	/* otherwise it is a bridge tap */
	else	{

		/* determine ABCD matrix for the bridge tap by
		   calling abcdf() recursively */
		abcdf(first->tap,f,&tmbridge); /* result stored in tmbridge */

		/* since the bridge tap looks like a shunt impedance,
		   the bridge tap admittance goes in C, unity in
		   A and D, and zero in B */
		tmsect.c = cdiv(tmbridge.c,tmbridge.a); /* admittance */
		tmsect.a = tmsect.d = unity;
		tmsect.b = zero;
		}

/* if twoport is the last in the line, replace *tf by tmsect */

	if(first->right == NULL)
		*tf = tmsect;

/* if twoport is not last in the line, call abcdf() recursively
   for the line on the right, and then multiply *tf by tmsect
   and replace *tf by the result */

	else {
		abcdf(first->right,f,tf);

		/* *tf now has the ABCD matrix for the line on
		   the right, now multiply by tmsect and put
		   result temporarily in tmbridge */
		tmbridge.a = cadd(cmpy(tmsect.a,tf->a),cmpy(tmsect.b,tf->c));
		tmbridge.b = cadd(cmpy(tmsect.a,tf->b),cmpy(tmsect.b,tf->d));
		tmbridge.c = cadd(cmpy(tmsect.c,tf->a),cmpy(tmsect.d,tf->c));
		tmbridge.d = cadd(cmpy(tmsect.c,tf->b),cmpy(tmsect.d,tf->d));

		/* now store the result back in *tf */
		*tf = tmbridge;
		}

	return;
}


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

			imped()

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

Function calculates the impedance of the line at each frequency

Note: this is the impedance including any source impedance;
	the impedance of the line itself can be obtained by
	subtracting the source impedance
*/

imped()
{
	int i;
	COMPLEX imped;
        COMPLEX num, den, R;

	R.real=110; R.imag=0;
	for(i=0; i<nm; ++i) {

	        num=cadd(tm[i].a, cdiv(tm[i].b, R));
		den=cadd(tm[i].c, cdiv(tm[i].d, R));
		imped = cdiv(num,den);

		fprintf(fimped,"%f %f %f\n",i/(m*T),imped.real,imped.imag);
		}
	return;
		
}

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

			ir()

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

Function calculates the voltage transfer function of the line from
the ABCD matrix, and prints out the log magnitude of that transfer
function to a file

The transfer function is multiplied by the complex-valued transfer
function filt_tf() provided by the user, which can be used to
include the effects of an input signal other than an inpulse
or an output equalizer

The impulse response, as modifed by filt_tf(), is then calculated
using the FFT and printed out to a file
*/

ir()
{
	POINTER calloc();
	COMPLEX filt_tf();
	double mag;
	double area;
	int i;
	COMPLEX *vt;	/* voltage transfer function of cable */
	COMPLEX temp;

/* allocate memory for transfer function array */

	vt = (COMPLEX *) calloc(nm,sizeof(COMPLEX));

/* determine the transfer function from the ABCD matrix --
   simultaneously print freq and trans fn in file ffreq */

	for(i=0; i<nm; ++i) {

		/* transfer function of line alone */
		vt[i] = cdiv(unity,tm[i].a);

		/* multiply by filt_tf()  */
		vt[i] = cmpy(vt[i],filt_tf(i/(m*T)));

		/* determine magnitude squared of transfer function */
		temp = cmpy(vt[i],conj(vt[i]));

		/* calculate magnitude of transfer function in dB --
		   but dont attempt to calculate attenuations larger
		   than 100 dB  */
		if(temp.real > 1.0e-20)
			mag = 10.*log10(temp.real);
		else
			mag = -200.;

		/* print out the voltage transfer function and
		   magnitude transfer function */
		fprintf(ffreq,"%f %e %e %f\n",i/(m*T),vt[i].real,vt[i].imag,mag);
		}

/* note the area of the pulse, which is vt[0] */

	area = vt[0].real;

/* calculate the impulse response */

	fft(vt,mexp+nexp);

/* print out time and impulse response at each frequency */

	for(i=0; i<nm; ++i)

		fprintf(ftime,"%-15g %-15g\n",i*T/n,
			(2.*vt[i].real - area)/(m*T));

	return;
}
