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

/* This routine finds the self impedance of a dipole of arbitrary
length 'length' , radius 'r' at a wavelength 'wavelength' 

Both the real part R and the imaginary part X are returned. R is
found from Antenna Theory, Analysis and Design, by Balans (Published
by Harper and Row). 

The real part, Rin is found by computing the radiation resistance Rr 
(equ 4-70, pp 124) then computing Rin=Rr/(sin(k l/2)^2, as given by 
equ 7-30a, pp 294. 

The reactive part, Xin is found by computing Xm (equ 7-33, pp 294)
then using xin=Xm/(sin(k l/2)^2 (equ 70-30b pp 294). */

void self(double r, double length, double lambda, double *Rin, double *Xin)
{
	double beta, mu, epsilon,eta,sin_bl,cos_bl,current_scale_factor;
	double bl, ci_bl, si_bl, ci_2bl, si_2bl,Xm,Rr;

	beta=2*M_PI/lambda;
	mu=4*M_PI*1e-7;
	epsilon=8.854187818e-12;
	eta=sqrt(mu/epsilon);
	bl=beta*length;

	/* We often need si(bl) , si(2*bl), ci(bl), ci(2*bl). These are
	time consuming to calculate, so I calculate just once. The Numerical
	recipes routine cisi, calculates both si and ci at the same time,
	so I wont use my routines ci and si, which both simply call cisi,
	but was one of the results */
	cisi(bl, &ci_bl, &si_bl);
	cisi(2*bl, &ci_2bl, &si_2bl);
	sin_bl=sin(bl);
	cos_bl=cos(bl);
	/* imaginary part of zin given by Balans, pp 294 */
	Xm=(eta/(4*M_PI)) * 
	(2*si_bl+ cos_bl*( 2*si_bl - si_2bl ) 
	-sin_bl*( 2*ci_bl-ci_2bl-ci(2*beta*r*r/length)) ); 
	current_scale_factor=sin(bl/2)*sin(bl/2); /* sin(bl/2)^2 */
	*Xin=Xm/current_scale_factor; /* Xin=Xm/(sin(bl/2)*sin(bl/2)); */
	/* The radiation resistance, as described by Balans C. A. Antenna theory,
	pp124, Published by Harper and Row, (1982)  */

	/* Rr=(eta/(M_PI*2.0))*(EULER+log(bl)-ci_bl+0.5*sin_bl*
	(si_2bl-2.0*si_bl)+0.5*cos_bl
	*(EULER+log(bl/2)+ci_2bl-2*ci_bl)); */
	Rr=(eta/(2.0*M_PI))*(EULER+log(bl)-ci_bl+0.5*sin_bl*(si_2bl-2.0*si_bl)+0.5*cos_bl*(EULER+log(bl/2.0)+ci_2bl-2*ci_bl));
	*Rin=Rr/current_scale_factor;
	/* printf("Rr=%f Rin=%f csf=%f\n",Rr,*Rin,current_scale_factor); */
	/* I'm a little bothered that Rin falls to ~60 Ohms at 
	resonance. Lawson suggested using 73, which would imply a large error.
	A graph in RadCom by G3SEK would suggest the 60 Ohms is too low. A 
	graph in Krauss, would suggest ~70 Ohms, so what is going on????
	I thought I'd use an equation in Krauss, but since this is the 
	Brown and King formuala given in Balans, that too gives ~60 Ohms.
	The formula in Krass ignores the diameter, but gives much the same values
	as the Balan's method */
	/* Now compute by method in Krauss, pp 421. */
	/* cot_bl_over_2=1.0/(tan(bl/2.0)); */

	/* *Rin=30.0*(   (1.0-(cot_bl_over_2*cot_bl_over_2))*Cin(2.0*bl)+4*cot_bl_over_2*cot_bl_over_2*Cin(bl)+2.0*cot_bl_over_2*(si_2bl-2.0*si_bl)        ); 
	printf("Rin=%f\n",*Rin); */
}