#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define INCR 100
#define POINTX 500
#define POINTY (-500)

#define P 2
#define Y 1
#define X 0

double setup(double distsep,int numofant,double maxgrid,double point[2],int choice);
double circleout(double ant_arry[][3],double maxgrid,int numofant,int choice);
double get_point(double ant_arry[][3],double point[],int numofant);
void get_phase(double ant_arry[][3],double point[2],int numofant);
void int_ant(double ant_arry[][3],int numofant,double distsep);
void ant_info_out(double ant_arry[][3],int numofant,double point[2]);

int main(int argc,char *argv[])
{
  double f=0,g,gold,fold;
  int numofant=10;
  double distsep=M_PI;
  double maxgrid=2000;
  double min=0,max=1,inc=1;
  double point[2];
  int i,switchpoint=0;
  char filename[1024],term[1024];
  point[X]=POINTX;
  point[Y]=POINTY;
  if (argc>1&&atoi(argv[1])>0)
    {
      switchpoint=atoi(argv[1]);
      switch (switchpoint) 
	{
	case 1: //1 to 50 antennas
	  min=1;max=50;inc=1;
	  break;
	case 2: //distance between antennas from Pi/16 to 16Pi
	  min=M_PI/16;max=16*M_PI;inc=M_PI/16;
	  break;
	case 8:
	  sprintf(filename,"ant0.eps");
	  printf(
		 "set terminal postscript eps enhanced color solid\n"
		 "\n"
		 "plot '-' smooth csplines title '%d element phased array-directivity chart' with lines 3\n" ,filename,numofant);
	  min=M_PI/16;max=4*M_PI;inc=(double)1/16;
	  break;
	case 7: //distance between antennas from Pi/16 to 4Pi (measureing directivity)
	  sprintf(filename,"ant0.gif");
	  printf(
		 "set terminal gif size 500,500 xffffff x000000 x404040 xff0000 xffa500  x0000ff xdda0dd x9500d3\n"
		 "set output '%s'\n"
		 "plot '-' smooth csplines title '%d element phased array-directivity chart' with lines 3\n" ,filename,numofant);
	  min=M_PI/16;max=4*M_PI;inc=(double)1/16;
	  break;
	case 3: //view radius from 10^0 to 10^3.5
	  min=0;max=3.5;inc=3.5/250;
	  break;
	case 4: //Point goes round in a circle
	  min=0;max=2*M_PI;inc=(2*M_PI)/175;
	  break;
	case 5: //Point goes away from center 10^0 to 10^3.5
	  min=0;max=3.5;inc=3.5/250;
	  break;
	case 6: //Point and view go away from center 10^0 to 10^3.5
	  min=0;max=3.5;inc=3.5/250;
	  break;
	case 9:
	case 0: //no inc
	default:
	  break;	
	}
    }
  if (argc>2&&atoi(argv[2])>0)
    numofant=atoi(argv[2]);
  if (argc>3&&atof(argv[3])>0)
    distsep=atof(argv[3]);
  if (argc>4&&atof(argv[4])>0)
    maxgrid=atof(argv[4]);
  if (argc>5)
    point[X]=atof(argv[5]);
  if (argc>6)
    point[Y]=atof(argv[6]);
  fprintf(stderr,"%d,%f,%f (%f,%f)\n",numofant,distsep,maxgrid,point[X],point[Y]);
  
  
  for(f=min;f<max;f+=inc)
    {
      switch(switchpoint) 
	{
	case 1: //1 to 50 antennas
	  numofant=f;
	  break;
	case 8:
	case 7:
	case 2: //distance between antennas from Pi/16 to 16Pi
	  distsep=f;
	  break;
	case 3: //view radius from 10^0 to 10^3.5(radians)
	  fprintf(stderr,"%g,%g\n",f,maxgrid=pow(10,f));
	  break;
	case 4: //Point goes round in a circle at 500 radius
	  point[X]=POINTX*cos(f);
	  point[Y]=POINTX*sin(f);
	  break;
	case 5: //Point goes away from center 10^0 to 10^3.5 
	  point[X]=pow(10,f)*cos(-M_PI/4);
	  point[Y]=pow(10,f)*sin(-M_PI/4);
	  break;
	case 6: //Point and view go away from center 10^0 to 10^3.5 
	  point[X]=pow(10,f)*cos(-M_PI/4);
	  point[Y]=pow(10,f)*sin(-M_PI/4);
	  fprintf(stderr,"%g,%g\n",f,maxgrid=pow(10,f));
	  break;
	case 9:
	case 0: //no inc
	default:
	  break;	
	}
      if(switchpoint!=7&&switchpoint!=8)
	{ 
	  if(switchpoint==9)
	    {
	      sprintf(filename,"\n",f);
	      sprintf(term,"set terminal postscript eps enhanced color solid\n");
	    }
	  else
	    {
	      sprintf(filename,"set output 'ant%08.4f.gif'\n",f);
	      sprintf(term, "set terminal gif size 500,500 xffffff x000000 x404040 xff0000 xffa500  x0000ff xdda0dd x9500d3\n");
	    }
	  
	  printf(
		 "set polar\n"
		 "set xtics axis nomirror\n"
		 "set ytics axis nomirror\n"
		 "set grid polar\n"
		 "set size square\n"
		 "%s"
		 "%s"
		 "plot [-2*pi:2*pi] [-1:1] [-1:1] '-' title '(%d,%g,%g)(num,sep,view) at (%g,%g)(radians)' with lines 3\n" ,term,filename,numofant,distsep,maxgrid,point[X],point[Y]);
	  fprintf(stderr,"Directivity: %f\n",setup(distsep,numofant,maxgrid,point,1));  
	  printf( "e\n");
	} 
      else 
	{
	  printf( "%f %f\n", distsep,setup(distsep,numofant,maxgrid,point,0));
	}
    }

}
double setup(double distsep,int numofant,double maxgrid,double point[2],int choice)
{
  double ant_arry[numofant][3];
  double tphase;
  int i;
  double out;

  int_ant(ant_arry,numofant,distsep);
  get_phase(ant_arry,point,numofant);

  ant_info_out(ant_arry,numofant,point);
  out=circleout(ant_arry,maxgrid,numofant,choice);

  return out;
}
void ant_info_out(double ant_arry[][3],int numofant,double point[2])
{
  int i;
  double tphase;
  for(i=0;i<numofant;i++)
    { 
      fprintf(stderr,"ANTenna %d postion %f,%f,%f\n",i,ant_arry[i][X],ant_arry[i][Y],ant_arry[i][P]);
    }
  tphase=get_point(ant_arry,point,numofant);
  fprintf(stderr,"Point(%f,%f) value = %f\n",point[X],point[Y],tphase);
}
double circleout(double ant_arry[][3],double maxgrid,int numofant,int choice)
{
  double lpoint[2];
  int i;
  double kf,jf,tphase,rphase,sum=0,max=0;
  for(i=0,rphase=0;rphase<2*M_PI;rphase+=(2*M_PI)/(500),i++)
    {
      jf=maxgrid*cos(rphase);
      kf=maxgrid*sin(rphase);
      lpoint[X]=jf;
      lpoint[Y]=kf;
      sum+=(tphase=get_point(ant_arry,lpoint,numofant));
      if(choice==1)
	printf("%f %f\n",rphase,tphase);
      if(tphase>max) max=tphase;
    }

  return max/(sum/i);//directivty Pmax/Pavg
}

double get_point(double ant_arry[][3],double point[],int numofant)
{
  int i,j,k;
  double lf;
  double out=0;
  double tout=0;
  double max=pow(numofant,2)/2;
  double tconst[numofant];
  for(i=0;i<numofant;i++)
    tconst[i]=hypot(point[X]-(ant_arry[i][X]),point[Y]-(ant_arry[i][Y]))+ant_arry[i][P];
  for(k=0,i=numofant-1;i>0;i--)
    {
      for(j=i;j>0;j--,k++)
	out+=cos(tconst[i]-tconst[j-1]);
    }
  out+=numofant/2;

  return (out/max);
}

void get_phase(double ant_arry[][3],double point[2],int numofant)
{
  int i;
  double phase0;
  for(i=0;i<numofant;i++)
    {
      ant_arry[i][P]=hypot(ant_arry[i][X]-point[X],ant_arry[i][Y]-point[Y]);
      if (i==0) phase0=ant_arry[0][P];
      ant_arry[i][P]=phase0-ant_arry[i][P];
    }
}
void int_ant(double ant_arry[][3],int numofant,double distsep)
{
  int i;
  double center[2];
  double angle,radius;
  center[X]=0;
  center[Y]=0;
  for(i=0;i<numofant;i++)
    {
      ant_arry[i][P]=0;
    }
  angle=(double)(2*M_PI)/numofant;
  radius=(double)1/(sin(angle/2)/(distsep/2));
   for(i=0;i<numofant;i++)
    {
      ant_arry[i][X]=cos(i*angle)*radius+center[X];
      ant_arry[i][Y]=sin(i*angle)*radius+center[Y];
    }
   if (numofant==1) 
     {
       ant_arry[0][X]=0+center[X];
       ant_arry[0][Y]=0+center[Y]; 
     }

}