#ifdef HAVE_STDLIB_H
#include <stdlib.h>
#endif
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include <memory.h>
#include <string.h>
#include <errno.h>
#include <ctype.h>
#include "nrutil.h"
#include "yagi.h"
#include "globals.h"

int main(int argc, char **argv)
{

	FILE *fp, *update_fp, *fp_out;
	char *input_filename, notes[1000];
	char *output_filename, *update_filename, *line, *ofile;
	int elements, driven, parasitic;
	int  ii;
	int  *indx, i, better=FALSE;
	double **A, *x,*v,*b,**z,pin,old_perf=-10000.0, new_perf=-1000.0;
	double **driven_data, **parasitic_data, angular_step, scale_factor=1.0;
	double **driven_data_tmp, **parasitic_data_tmp;
	double **driven_data_global_best, **parasitic_data_global_best;
	double design_frequency, min_frequency, max_frequency, step_frequency;
	double E_fwd, E_back, H_fwd, H_back, frequency, old_SD_of_currents=1e100;
	struct flags flag;
	struct FCOMPLEX *voltage, *current, input_impedance, z_centre;
	struct element_data *coordinates;
	struct performance_data mean, best_performance,start,changes,worst;
	/* Zero all the flags in the structure flag. Some compilers allow
         me to zero them when I define it in the header file, but not all.
         Hope this method is less compiler specific. */
	memset(  (char *) &flag , 0, sizeof(flag));
	/* Set the structs to have sensible values to start with. */
   set_performance_structures(&weight,&max,&best_performance,&worst);
	opterr=0;
	ofile=string(0L,100L);
	/* Since there are so many command line options to get, as of version 1.10
	of the software, all the options are got from a sepparate file, using lots
	of global variables! */
	flag.eflg=31;    /* set elements to move */
	flag.oflg=32805; /* Set things to optimise for */
	test_for_stop_file(); /* Exit if the file "stop" exists */
	get_command_line_options(argc, argv,&flag);
	if(!flag.Tflg && !flag.tflg)
		optimising_for(flag);
	/* a function check_flags() does all the checking to make sure flags 
	set (by user adding options) are reasonable. It also acts on many of them. */
	check_flags(flag, argc, optind, argv[0]);
	iterations=atoi(argv[optind+1]);
	if(flag.Aflg)
		iterations=2;
	/* allocate memory for strings */
	input_filename=string(0L,100L);
	line=string(0L,MAX_LINE);
	output_filename=string(0L,100L);
	update_filename=string(0L,100L);
	seedRNG(); 
	/* Create a file showing best gain, fb-ratio etc achieved to date */
	/* This might typically have the name 'antenna.up' */ 
	strcpy(update_filename,argv[optind]); 
	strcat(update_filename,".up");
	update_fp=fopen(update_filename,"wt"); /* Remove filname.up */
	fclose(update_fp);
	strcpy(input_filename,argv[optind]);
	fp=fopen(input_filename,"rt");
	if(fp == NULL)
	{
		fprintf(stderr,"Cant open %s\n", input_filename);
		exit(10);
	}
	else
		fclose(fp);
	strcpy(output_filename,input_filename);
	strcat(output_filename,".bes");
	fp_out=fopen(output_filename,"wb");  /* Remove filname.bes */
	fclose(fp_out);
	/* since the size of all the arrays depends on the number of elements,
	lets first find out how many elements we have, including driven and 
	parasitic sepparately.  Then try to allocate all the memory we need. If we 
	fail, then not much time is wasted. */
	elements=get_number_of_elements(input_filename, &driven, &parasitic);
	/* allocate all memory */
	driven_data=dmatrix(1L,(long)(driven),1L,6L); 
	parasitic_data=dmatrix(1L,(long) (parasitic),1L,4L); 
	driven_data_global_best=dmatrix(1L,(long)(driven),1L,6L); 
	parasitic_data_global_best=dmatrix(1L,(long) (parasitic),1L,4L); 
	driven_data_tmp=dmatrix(1L,(long)(driven),1L,6L); 
	parasitic_data_tmp=dmatrix(1L,(long) (parasitic),1L,4L); 
	z=dmatrix(1L,(long) elements, 1L, elements*2L);
	A=dmatrix(1L, 2L*elements, 1L, 2L*elements);   /* A.x=b for large matrices */
	x=dvector(1L, 2L*elements);
	b=dvector(1L, 2L*elements);

	current=FCOMPLEXvector(1L,(long) elements); 
	voltage=FCOMPLEXvector(1L,(long) elements);

	coordinates=element_data_vector(1,elements); /* x,y &l of centre of ele's */
	v=dvector(1L,2L*elements); 
	indx=ivector(1,2L*elements);

	read_yagi_data(line, input_filename, &design_frequency, &min_frequency,                &max_frequency, &step_frequency, driven, driven_data, parasitic,                parasitic_data,&angular_step);

	/* make a copy of the element data, as we will randomise it 
	Copy from 'parasitic_data' to 'parasitic_data_tmp' etc */
	errno=0;
	/* cp _data to _data_global_best */
	copy_matrix(driven,6,driven_data_global_best, driven_data);
	copy_matrix(parasitic,4,parasitic_data_global_best,parasitic_data); 
	/* cp _data to _data_tmp */
	copy_matrix(driven,6,driven_data_tmp, driven_data);
	copy_matrix(parasitic,4,parasitic_data_tmp,parasitic_data); 
	/* Find the maximum likely gain, given the boom length */
	max.gain=determine_maximum_gain2(elements);
	if(flag.Kflg==0)
		K_times_max=1;
	/* If a genetic algorithm is used, call the genetic algorithm
	and exit */
	if(flag.gflg) /* Do this if using the genetic algorithm */
	{
		genetic_algorithm(output_filename, update_filename, flag,\
		 design_frequency, min_frequency, max_frequency, step_frequency, \
		 angular_step, driven, parasitic, driven_data, parasitic_data, v, z, \
		 &pin, voltage, current, &input_impedance, coordinates, A, b, indx, \
		 &mean); 
		printf("The best design is in a file \"%s\". You should check it thoroughly\n",output_filename);
		printf("and if its better than %s, copy %s to %s\n",input_filename,\
		 output_filename, input_filename);
		exit (0);
	}
	if(flag.wflg) /* Do for wide band use */
		sprintf(notes,"This has been run through optimise and optimised for \
		 wide-band use over the frequency range %.3f MHz to %.3fMHz. \n", \
		 min_frequency/1e6, max_frequency/1e6);
	for(i=1; i<=iterations; ++i) /* main loop */
	{
		percent=change_max_percentage_changes(i,iterations,original_percent);
		/* By default flag.wflg is set to zero, hence the following loop
		will be exectude only once. If the flag was set to +2, the loop
		is executed 3 times and we get the average performance */
		for(ii=-1; ii<flag.wflg;++ii)
		{
			if(ii==-1)
			{
				frequency=design_frequency;
				/* Zero mean structure */
				memset((char *) &mean,0,sizeof(mean));
			}
			if(ii==0)
				frequency=min_frequency;
			if(ii==1)
				frequency=max_frequency;
			solve_equations(frequency, driven, parasitic, driven_data, parasitic_data, v, z, &pin, voltage, current, &input_impedance, coordinates, A, b, indx);
			if(ii==-1)
				z_centre=input_impedance;
		   /* compute gain at theta=90, phi=0 (forward direction) */
			gain(90,0,pin,frequency/design_frequency,coordinates,current,elements,&E_fwd,&H_fwd,frequency,design_frequency);
	   	/* now compute gain in the reverse direction */
			gain(270,0,pin,frequency/design_frequency,coordinates,current, \
			elements,&E_back,&H_back,frequency,design_frequency);
			/* put the performance in a structure performance_data */
			set_mean_structure(input_impedance,E_fwd,E_back,flag, \
			pin, coordinates, current, elements,frequency,design_frequency,&mean);
			if(i==1)
				start=mean;
		}
		if(flag.wflg)
		{
			mean.gain/=3.0; mean.fb/=3.0; mean.swr/=3.0; 
		}
		if(flag.Wflg || flag.gflg)
		{
			old_perf=performance(flag, best_performance, max, weight,start);
			new_perf=performance(flag, mean, max, weight,start);
			if(new_perf > old_perf)
				better=TRUE;
			else
				better=FALSE;
			if((better==TRUE) && flag.Cflg)
				better=linear_current_optimisation_test(current, &old_SD_of_currents, elements, parasitic, flag);
		}
		else if(flag.Aflg)
			better=TRUE;
		else
		{
			better=is_it_better(flag.oflg,mean,best_performance);
			if((better==TRUE) && flag.Cflg)
				better=linear_current_optimisation_test(current, &old_SD_of_currents, elements, parasitic, flag);
		}
		if(better==TRUE && !flag.Tflg)
		{
			/* The function do_since_better() does most things that need to
			be done if we have a better design. These include writing the 
			new antenna description to disk, print the new performance to disk 
			and to stdout, etc */  

			/* We write the input impedance at the centre of the band - not
			averaging, but the swr is averaged */

			do_since_better(i,output_filename, update_filename, z_centre, \
			 mean, flag, notes, design_frequency, min_frequency, \
			 max_frequency, step_frequency, elements, driven, parasitic, \
			 angular_step,driven_data,parasitic_data,scale_factor, new_perf); 
			/* now set the best performance to the performance of this iteration*/

			best_performance=mean;
			/* Copy current new design (_data) to the place we keep the best
			design (_data_tmp) */
			/* copy _data to _data_global_best */
			copy_matrix(driven,6,driven_data_global_best,driven_data);
			copy_matrix(parasitic,4,parasitic_data_global_best,parasitic_data);
			K_times=0;
		}
		else if (better == FALSE || flag.Tflg) /* no improvement */
		{
			/* Restored original data, without any randomisation */
			K_times++;
			if(K_times<K_times_max) /* usual case */
			{
				/* Retriev data from tmp into data */
				copy_matrix(driven,6,driven_data,driven_data_tmp);
				copy_matrix(parasitic,4,parasitic_data,parasitic_data_tmp);
			}
			else
			{
				/* Retrieve data from _global_best to data, then put
				data in tmp */
				copy_matrix(driven,6,driven_data,driven_data_global_best);
				copy_matrix(parasitic,4,parasitic_data,parasitic_data_global_best);
				copy_matrix(driven,6,driven_data_tmp,driven_data);
				copy_matrix(parasitic,4,parasitic_data_tmp,parasitic_data);
				K_times=0;
			}
			if(flag.Tflg) /* Want to find the sensitivity to constuction errors */
			{
				if(mean.gain<worst.gain)
					worst.gain=mean.gain;
				if(mean.fb<worst.fb)
					worst.fb=mean.fb;
				if(mean.swr > worst.swr)
					worst.swr=mean.swr;
				if(mean.sidelobe < worst.sidelobe)
					worst.sidelobe=mean.sidelobe;
			}
		}
      /* The following line looks for a file 'stop' in the current directory
		and stops next time if it finds it. For  DOS use, calling the function
		kbhit() - keyboard hit, might be better, so the user kits the keyboard
		to stop the program. We pause after 100 iterations. With a much faster/
		slower machine, you might want to change this */
		end_if_stop_exists(&i,iterations,100);
		dynamic_changing_of_weights(i,1,&weight);
		if(!flag.Tflg && !flag.Aflg)/*user doesnt want to test sensitivity */  
			randomise(flag.eflg, design_frequency/1e6, percent, driven_data, parasitic_data, driven, parasitic); 
		else if(flag.Aflg)
			automatic_enhancement(flag,design_frequency,driven_data,parasitic_data, driven,parasitic,voltage, current, &input_impedance, v,z,A,b,indx,coordinates);
		else if(flag.Tflg) /* Check sensitivity to construction errors */
			sensitivity(boom_sd, length_sd, driven_data, parasitic_data, driven, parasitic); 

	} /* End of main loop, now finished trying to better the design */
	if(!flag.Tflg)
	{
		printf("The best design is in a file \"%s\". You should check it thoroughly\n",output_filename);
		printf("and if its better than %s, copy %s to %s\n",input_filename,output_filename, input_filename);
		printf("For your information, the original data on the antenna was:\n");
		print_relavent_performance_data(stdout,"Start data:",0,flag,start,0.0,0,0);
		print_relavent_performance_data(stdout,"Final data:",0,flag,best_performance,0.0,0,0); 
		changes=subtract_structures(best_performance,start);
		print_relavent_performance_data(stdout,"Changes:   ",0,flag,changes,0.0,0,0); 
	}
	else if (flag.Tflg) /* find worst design */
	{
		printf("For your inforation, the original data on the antenna was:\n");
		print_relavent_performance_data(stdout,"Start data:",0,flag,start,0.0,0,0);
		print_relavent_performance_data(stdout,"Worst data:",0,flag,worst,0.0,0,0); 
		changes=subtract_structures(worst,start);
		print_relavent_performance_data(stdout,"Changes:   ",0,flag,changes,0.0,0,0); 
		printf("\nTolerance parameters were: Length_SD(t)=%.3fmm Boom_SD(T) = %.3fmm\n",length_sd,boom_sd);
	} 

	/* free strings */
	free_string(input_filename, 0L,100L);
	free_string(line, 0L, MAX_LINE);
	free_string(update_filename, 0L, 100L);
	free_string(output_filename,0L, 100L);
	free_string(ofile,0L,100L);
	/* free vectors */
	free_dvector(x,1L, 2L*elements);
	free_dvector(b,1L, 2L*elements);
	free_dvector(v,1L,2L*elements); 
	free_ivector(indx, 1L, 2L*elements);
	/* free matrices */
	free_dmatrix(z,1L,(long) elements, 1L, elements*2L);
	free_dmatrix(A,1L, 2L*elements, 1L, 2L*elements);
	free_dmatrix(driven_data,1L,(long)(driven),1L,6L); 
	free_dmatrix(parasitic_data,1L,(long) (parasitic),1L,4L); 
	free_dmatrix(parasitic_data_global_best,1L,(long) (parasitic),1L,4L); 
	free_dmatrix(driven_data_global_best,1L,(long) (parasitic),1L,6L); 
	free_dmatrix(driven_data_tmp,1L,(long)(driven),1L,6L); 
	free_dmatrix(parasitic_data_tmp,1L,(long) (parasitic),1L,4L); 
	/* free FCOMPLEX vectors */
	free_FCOMPLEXvector(current, 1L,(long) elements); 
	free_FCOMPLEXvector(voltage, 1L,(long) elements);
	free_element_data_vector(coordinates,1,elements);
	exit(0);
}