#ifdef HAVE_STDLIB_H
#include <stdlib.h>
#endif
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include "yagi.h"

extern int errno;

void mutual_impedance(int i, int j, double frequency, int driven, int parasitic, double **d, double **p, double **impedance)
{
	double xi, yi, xj, yj, li, lj, s, lamda, real, imag, mean_length;
	int col;
	/* first find location x,y of first element */
	if(i<=driven)				/* i is a  driven element*/
	{
		xi=d[i][X];
		yi=d[i][Y];
		li=d[i][LENGTH];
	}
	else if (i>driven)			/* i is a  parasitic  element */
	{
		xi=p[i-driven][X];
		yi=p[i-driven][Y];
		li=p[i-driven][LENGTH];
	}
	if(j<=driven)				/* j is a  driven element*/
	{
		xj=d[j][X];
		yj=d[j][Y];
		lj=d[j][LENGTH];
	}
	else if (j>driven)			/* i is a  parasitic  element */
	{
		col=j-driven;
		xj=p[col][X]; 
		yj=p[col][Y];
		lj=p[col][LENGTH];
	}
	/* compute distance from element i to element j */
	s=sqrt((xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)); /* in metres */
	mean_length=(li+lj)/2.0;
	lamda=3e8/frequency; 
	z21(lamda, s , mean_length , &real, &imag);
	/* fill in Zij and Zji at the same time, as matrix is symmetric. */
	impedance[i][(2*j)-1]=real; /* real; */
	impedance[i][2*j]=imag; /* imag; */
	impedance[j][2*i-1]=real; /* real; */
	impedance[j][2*i]=imag; /* imag; */

#ifdef DEBUG
	if(errno)
	{
		fprintf(stderr,"Errno =%d in  mutual.c\n", errno);
		exit(1);
	}
#endif
}

/* To find the mutual impedance between to arbitary length, thin elements
I used the equations on Krauss, Antennas, McGraw Hill, 1988, pp426 and
427. Original work from Brown and King, 'High Frequency Models in Antenna
Investigations', Proc IRE, vol 22, pp457-480, April 1934*/
void z21(double lamda, double d, double l, double *r21, double *x21)
{

	double  b, cos_bl, sin_bl, sin_bl_over_2, cos_bl_over_2, c, s ;
	double t1, t2, t3, t4;
	double si_t1, ci_t1, si_t4, ci_t4, ci_bd, si_bd;
	double ci_t2, si_t2, ci_t3, si_t3;

	b=M_PI*2/lamda;
	t1=b*(sqrt(d*d+l*l)+l);
	t2=0.5*b*(sqrt(4*d*d+l*l)-l);
	t3=0.5*b*(sqrt(4*d*d+l*l)+l);
	t4=b*(sqrt(d*d+l*l)-l);
	/* To save findinding the same slow trigometric and ever slower
	si and ci functions, I'll just look them up once */
	cos_bl=cos(b*l);
	sin_bl=sin(b*l);
	sin_bl_over_2=sin(b*l/2);
	cos_bl_over_2=cos(b*l/2);
	s=sin_bl_over_2*sin_bl_over_2;
	c=cos_bl_over_2*cos_bl_over_2;

	cisi(t1, &ci_t1, &si_t1);
	cisi(t2,&ci_t2, &si_t2);
	cisi(t3,&ci_t3, &si_t3);
	cisi(t4, &ci_t4, &si_t4);
	cisi(b*d,&ci_bd, &si_bd);
	/* Real part of mutual impedance, computed as equation on page
	426 of Kraus */
	*r21=(30/s)*(2*(2+cos_bl)*ci_bd
	-4*c*( ci_t2 + ci_t3 )
	+cos_bl*( ci_t4 + ci_t1 )
	+sin_bl* ( si_t1 - si_t4 -2*si_t3 +2*si_t2 ) );

	/* Imaginary part of mutual impedance, computed as equation on page
	427 of Kraus */
	*x21=(30/s)*(-2*(2+cos_bl)*si_bd
	+4*c*( si_t2 + si_t3 )
	-cos_bl*( si_t4 + si_t1 )
	+sin_bl* ( ci_t1 - ci_t4 -2*ci_t3 +2*ci_t2 ) );
#ifdef DEBUG
	if(errno)
	{
		fprintf(stderr,"Errno =%d in  z21() mutual.c\n", errno);
		exit(1);
	}
#endif
}