#ifdef HAVE_STDLIB_H
#include <stdlib.h>
#endif
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include "yagi.h"
#define TINY 1e-10

extern int errno;

/* This function finds the gain both in the E plane (xz plane) and the
H plane (xy plane) at angle (theta, phi). The method used is as described
on page 1-12 of 'Yagi Antenna Design' by Dr. Lawson , ARRL */

void gain(double theta, double phi, double pin, double F, struct element_data *coordinates, struct FCOMPLEX *current, int elements, double *gain_E_plane, double *gain_H_plane, double actual_frequency, double design_frequency) 
{
	int i;
	double *r_E, *r_H, *g_E, *g_H,integer_bit;
	double length, x, y, lamda_design, lamda, tmp;
	struct FCOMPLEX temp_E, temp_H, e_gain, h_gain;
	/* need to allocate space for FCOMPLEX types. Since there is no Numerical 
	Recipes routine, I'll use the standard method. Since I've always used
	elements positions 1 to N, I'll just wast the extra location */
	r_E=dvector(1L,(long) elements); 
	r_H=dvector(1L,(long) elements);
	g_E=dvector(1L,(long) elements);
	g_H=dvector(1L,(long) elements);
	e_gain.r=0;   /* set to zero. Necessary as we sum into these */
	e_gain.i=0;
	h_gain.r=0;
	h_gain.i=0;
	/* convert theta and thi to radians from degrees */
	theta=theta*M_PI/180;
	phi=phi*M_PI/180;
	lamda_design=3e8/design_frequency;
	lamda=3e8/actual_frequency;
	for(i=1;i<=elements;++i) 
	{
		length=coordinates[i].length/lamda;
		x=coordinates[i].x/lamda_design;
		y=coordinates[i].y/lamda_design;

		/* for E -plane */
		r_E[i]=sin(theta)*x;
		if(fabs(theta) < TINY)    /* avoid division by zero if theta=0 */
			g_E[i]=0;
		else
			g_E[i]=(cos(M_PI*length*cos(theta))-cos(M_PI*length))/sin(theta);

		/* for H -plane */
		r_H[i]=x*cos(phi)+y*sin(phi);
		g_H[i]=1-cos(M_PI*length);
		/* printf("g_H[%d]=%.16lf \n", i, g_H[i]);		 */
		temp_E.r=0;
		temp_E.i=2*M_PI*r_E[i]*F;
		temp_E=E_to_complex_power(temp_E); /* exp(j 2 pi r[i] F) */

		temp_H.r=0;
		temp_H.i=2*M_PI*r_H[i]*F;
		temp_H=E_to_complex_power(temp_H);
      /* printf("element %d temp_H.r=%f temp_H.i=%f\n",i, temp_H.r, temp_H.i); */
		/* get element currents */
		temp_E=Cmul(temp_E,current[i]);
		temp_H=Cmul(temp_H,current[i]);
		/* printf("element %d: current = %.10lf i%.10lf = %.10lf (mag) at %f degrees\n", i,current[i].r, current[i].i, sqrt(current[i].r*current[i].r+current[i].i*current[i].i), 180*atan2(current[i].i,current[i].r)/M_PI);         */
      /* printf("element %d I temp_H.r=%.10lf temp_H.i=%.10lf\n",i, temp_H.r, temp_H.i); */
		e_gain=Cadd(e_gain,RCmul(g_E[i],temp_E));
		h_gain=Cadd(h_gain,RCmul(g_H[i],temp_H));
		/* printf("h_gain=%.16lf + i %.16lf \n", h_gain.r, h_gain.i); */
	}
	/* there will be a divide  by zero in calculating
	equ 1.28, of Lawsons book, (see page 1-12)
	at theta=0,180,360 degree, infact anywhere where 
	sin(theta) is zero */

	if( modf(theta/M_PI,&integer_bit) < TINY || fabs(theta-M_PI)<TINY)
	{
		*gain_E_plane=-150.0;   /* a very very small number */
	}
	else
	{
		tmp=(Cabs(e_gain)*Cabs(e_gain)*60/pin);  
		*gain_E_plane=10*log10(tmp);  
#ifdef DEBUG
	if(errno)
	{
		fprintf(stderr,"Errno =%d in gain() of gain.c after E\n", errno);
		printf("tmp=%f pin=%f theta=%f gainE=%f\n", tmp, pin,theta,*gain_E_plane);
		exit(1);
	}
#endif
	}
	tmp=(Cabs(h_gain)*Cabs(h_gain)*60/pin);  
	*gain_H_plane=10*log10(tmp);  

#ifdef DEBUG
	if(errno)
	{
		fprintf(stderr,"Errno =%d in gain() of gain.c after H\n", errno);
		printf("tmp=%f pin=%f \n", tmp, pin);
		exit(1);
	}
#endif
	free_dvector(r_E,1L,(long) elements); 
	free_dvector(r_H,1L,(long) elements); 
	free_dvector(g_E,1L,(long) elements); 
	free_dvector(g_H,1L,(long) elements); 

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

struct FCOMPLEX E_to_complex_power(struct FCOMPLEX x)
{
	struct FCOMPLEX a, answer;
	double real;

	/* First real part */
	real=exp(x.r);
	/* now the imaginary part */
	/* Exp(i theta) = cos(theta) + i sin(theta) */
	a.r=cos(x.i);
	a.i=sin(x.i);
	answer=RCmul(real,a);
	return(answer);
}
#undef TINY