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

extern int errno;

int main(int argc, char **argv)
{
	double r,x,d;
	for(d=0.0000001;d<=1;d+=0.02)
	{
		z21(1.0,d,0.5, &r, &x);
		printf("%f %f %f %f \n", d,r,x, sqrt(r*r+x*x));
	}
	exit(0);
}

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

double ci(double x)			/* cosine integral */
{
	double null, result;

	cisi(x, &result, &null);
	return(result);
}