~ubuntu-branches/ubuntu/precise/yagiuda/precise

#ifdef HAVE_STDLIB_H

#include <stdlib.h>

#endif

#include <stdio.h>

#include <math.h>

#include <malloc.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