1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4 /*
5 A few changes have been made, including:
6 1) Memory is freed, rather than lying on the OS
7 2) Memory is no longer repeatdly allocated, which meant excessive use
8 of system memory.
9 */
10 #include <stdio.h>
11 #ifdef HAVE_STDLIB_H
12
13 #endif
14 #include <math.h>
15 #include <string.h>
16 #include <sys/types.h>
17 #include "yagi.h"
18
19 double Zo=Z0; /* Z0 is defined in yagi.h */
20
21 extern int optind, opterr;
22 extern char *optarg;
23 double min_offset_from_peak, angular_stepsize_2; /* not needed */
24
25 double pin, design_f, normalised_f, f;
26 double gain_E_plane, gain_H_plane, peak_gain;
27 struct element_data *coordinates;
28 struct FCOMPLEX *current;
29 int elements;
30 double rrange_min=RRANGE_MIN;
31 double rrange_max=RRANGE_MAX;
32
33
main(int argc,char ** argv)34 int main(int argc, char **argv)
35 {
36 FILE *ifp, *ofp, *gain_fp;
37 FILE *gnuplot_filename_fp;
38 FILE *gnuplot_log_command_filename_fp;
39 FILE *gnuplot_lin_command_filename_fp;
40 char *input_filename, *output_filename, *temp;
41 char *gnuplot_log_command_filename;
42 char *gnuplot_lin_command_filename;
43 char *gain_filename;
44 char *gnuplot_filename;
45 char *original_filename;
46 double min_f, max_f, step_f, vswr, angular_step,theta;
47 double magnitude, phase,max_SL,E_max=E_MAX, H_max=H_MAX;
48 struct FCOMPLEX *voltage, input_impedance;
49 double fb_ratio,gain_E_plane_back;
50 double gain_H_plane_back;
51 struct flags flag;
52 struct pattern pattern_shape;
53 int c, step=0;
54 input_filename=string(0L, 100L);
55 gnuplot_filename=string(0L, 400L);
56 gnuplot_log_command_filename=string(0L, 400L);
57 gnuplot_lin_command_filename=string(0L, 400L);
58 output_filename=string(0L, 100L);
59 original_filename=string(0L, 100L);
60 gain_filename=string(0L, 100L);
61 memset((char *) &flag, 0, sizeof(flag));
62 while ((c = getoptions(argc,argv,"cehpsr:E:H:R:Z:")) != -1)
63 switch (c)
64 { case 'c':
65 flag.cflg=1;
66 break;
67 case 'p': /* Print data for gnuplot_log */
68 flag.pflg=1;
69 break;
70 case 'e': /* suppress E- plane 3dB BW */
71 flag.eflg=1;
72 break;
73 case 'h': /* suppress H- plane 3dB BW */
74 flag.hflg=1;
75 break;
76 case 'r': /* suppress H- plane 3dB BW */
77 rrange_min=atof(optarg);
78 break;
79 case 'R': /* suppress H- plane 3dB BW */
80 rrange_max=atof(optarg);
81 break;
82 case 's': /* suppress diagnostics */
83 flag.sflg=1;
84 break;
85 case 'E': /* Max angle in E plane to look for -3dB */
86 E_max=atof(optarg);
87 break;
88 case 'H':
89 H_max=atof(optarg); /* max angle in H plane to look for -3dB */
90 break;
91 case 'Z': /* Zo */
92 Zo=atof(optarg);
93 break;
94 case '?': /* default */
95 flag.errflg=1;
96 break;
97 }
98 if((argc-optind != 1) || flag.errflg)
99 {
100 usage_output(argv[0]);
101 exit(0);
102 }
103 /* The file read by 'output' is expected to have a .out extension.
104 'output' creates files with extensions:
105 .dat Containing the basic broperties of impedance, gain etc as a function
106 of frequency.
107 .gai Contains field patterns (gain) in the both H and E plane.
108
109 We here get all the file names, then open files. */
110 strcpy(input_filename, argv[optind]);
111 strcpy(original_filename, argv[optind]);
112 strcpy(output_filename, *(argv+optind));
113 strcpy(gain_filename, *(argv+optind));
114
115 strcat(input_filename,".out");
116 strcat(output_filename,".dat");
117 strcat(gain_filename,".gai");
118 ifp=fopen(input_filename,"rb");
119 ofp=fopen(output_filename,"wb");
120 gain_fp=fopen(gain_filename,"wt");
121 if(ifp == NULL)
122 {
123 fprintf(stderr,"Cant open %s\n", input_filename);
124 exit(10);
125 }
126
127 temp=string(0L,100L);
128 /* No longer need to refer to files by names, we just use the
129 file pointers, so free memory associated with names. */
130 free_string(output_filename,0L,100L);
131 free_string(gain_filename,0L,100L);
132 /* Read the header on the .out file, which is usually 100 Bytes,
133 but can be defined in yagi.h as HEADER_SIZE */
134 elements=read_header(ifp, ofp, &min_f, &max_f, &step_f, &design_f, &angular_step);
135 /* Creatre an array of structures, so contain the x and y locations of
136 each element, and the lengths of the elements. All dimensions are in m */
137 coordinates=element_data_vector(1,elements); /* x,y &l of centre of ele's */
138 /* read into memory the x,y, and lengths of all elements */
139 fread((char *) &coordinates[1], sizeof(struct element_data),elements,ifp);
140 /* Create arrays to hold real and imaginary components of voltage */
141 voltage=FCOMPLEXvector(1,elements);
142 current=FCOMPLEXvector(1,elements);
143 /* read voltages applied to elements into memory */
144 fread((char *) &voltage[1], sizeof(struct FCOMPLEX),elements, ifp);
145 /* Now print a one line header to put into .dat file */
146 fprintf(ofp," f(MHz) E(deg) H(deg) R jX VSWR Gain(dBi) FB(dB) SideLobes(dB)\n\n");
147 if(flag.pflg==1)
148 {
149 if(angular_step > 10)
150 {
151 fprintf(stderr,"The angular step site is set at less that 10 degress. Edit the file %s, change the ANGULAR_STEP to 10 degrees or less, then rerun 'yagi %s' and 'output -p %s'\n", original_filename, original_filename, original_filename);
152 exit(12);
153
154 }
155 sprintf(gnuplot_log_command_filename,"%s.glog",original_filename);
156 sprintf(gnuplot_lin_command_filename,"%s.glin",original_filename);
157 gnuplot_log_command_filename_fp=fopen(gnuplot_log_command_filename,"w");
158 gnuplot_lin_command_filename_fp=fopen(gnuplot_lin_command_filename,"w");
159 write_gnuplot_header(gnuplot_log_command_filename_fp,f,original_filename,step,LOG);
160 write_gnuplot_header(gnuplot_lin_command_filename_fp,f,original_filename,step,LIN);
161 }
162 /* compute data at every frequency of interest */
163 for(f=min_f; f<=max_f; f+=step_f)
164 {
165 normalised_f=f/design_f;
166 fread((char *) ¤t[1], sizeof(struct FCOMPLEX),elements, ifp);
167 z_input(voltage[1], current[1], &input_impedance);
168 /* Compute magnitude and phase of reflection coefficient */
169 reflection_coefficient(input_impedance, &magnitude, &phase);
170 /* Compute VSWR from magnitude of reflection coefficient */
171 vswr=calculate_vswr(magnitude);
172 /* Calculate power input to antenna */
173 pin=calculate_power_input(input_impedance.r,current[1]); /* in Watts */
174 if(pin < 0.0)
175 {
176 fprintf(stderr,"input power less than 0! Probably due to self or mutual impedance being inaccurate when current not sinusidal\n");
177 fprintf(stderr,"f=%f MHz Zin=%f +i%f\n",f/1e6, input_impedance.r, input_impedance.i);
178 }
179 /* compute gain at theta=90, phi=0 (forward direction) */
180 gain(90, 0, pin, normalised_f, coordinates, current, elements, &gain_E_plane, &gain_H_plane, f, design_f );
181 /* compute gain at theta -90 degrees (ie back - direction) */
182 gain(270, 0, pin, normalised_f, coordinates, current, elements, &gain_E_plane_back, &gain_H_plane_back, f, design_f);
183 fb_ratio=gain_E_plane-gain_E_plane_back; /* in dB, so subtract */
184 /* If the values of impedance are too large, the cant be printed in
185 the .dat file. Hence we truncate them. At this point, the numbers
186 are so large, that we are not likely to be able to use them anyway. */
187 if(input_impedance.r > 9999999)
188 input_impedance.r = 9999999;
189 if(input_impedance.i > 9999999)
190 input_impedance.i = 9999999;
191 if(input_impedance.r < -9999999)
192 input_impedance.r = -9999999;
193 if(input_impedance.i < -9999999)
194 input_impedance.i = -9999999;
195
196 /* Find the most importent bits of the pattern info */
197
198 /* Estimiate 3dB bandwidth using .See '-E' and '-H' options */
199 peak_gain=gain_E_plane;
200 if(flag.hflg || elements==1)
201 pattern_shape.three_dB_H=0.0;
202 else
203 pattern_shape.three_dB_H=(double) zbrent(error_3dB_H, 0.0,H_max,0.1);
204 if(flag.eflg)
205 pattern_shape.three_dB_E=0.0;
206 else
207 pattern_shape.three_dB_E=(double) zbrent(error_3dB_E, 90.0,E_max,0.1);
208 if(flag.cflg==1)
209 {
210 max_SL=find_max_sidelobe_slow(peak_gain, pin,coordinates,current,elements,f,design_f);
211 }
212 else
213 max_SL=0.0;
214 fprintf(ofp,"%9.3f %3.1f %3.1f %6.2f %6.2f %6.3f %10.3f %10.3f %10.3f\n",f/1e6,2*(pattern_shape.three_dB_E-90), 2*pattern_shape.three_dB_H,input_impedance.r, input_impedance.i, vswr, peak_gain, fb_ratio,max_SL);
215 write_gain_at_various_angles(gain_fp, angular_step, pin, normalised_f, f, coordinates, current, elements, design_f);
216 if((min_f < max_f) && !flag.sflg)
217 printf("%s completed %5.1f%%\n",argv[0],100*(f-min_f)/(max_f-min_f));
218 if(flag.pflg == 1) /* Plot results in files suitable for gnuplot */
219 {
220 step++;
221 sprintf(gnuplot_filename,"%s.%.4f.g",original_filename,f/1e6);
222 gnuplot_filename_fp=fopen(gnuplot_filename,"w");
223
224 for(theta=0;theta<360;theta+=angular_step)
225 {
226 gain(theta+90.0, 0, pin, normalised_f, coordinates, current, elements, &gain_E_plane, &gain_H_plane, f, design_f );
227 if(gain_E_plane < rrange_min)
228 gain_E_plane= rrange_min;
229 if(gain_H_plane < rrange_min)
230 gain_H_plane= rrange_min;
231 fprintf(gnuplot_filename_fp,"%11.3f %11.3f %11.3f %8.3f %11.3f %11.3f %11.3f %11.3f %11.3f\n",theta,gain_E_plane,gain_H_plane,vswr,pow(10.0,gain_E_plane/10.0),pow(10.0,gain_H_plane/10.0),0.0,0.0,0.0);
232 /*
233 if(theta<=90)
234 zz=-40.0;
235 else if (theta > 90 && theta <= 180)
236 zz=-20.0;
237 else if (theta > 180 && theta <= 270)
238 zz=0.0;
239 else
240 zz=10.0;
241 fprintf(gnuplot_log_filename_fp,"%11.3f %11.3f %11.3f %8.3f %11.3f %11.3f\n",theta,zz,gain_H_plane,vswr,0.0,0.0,0.0);
242 */
243 }
244 fprintf(gnuplot_log_command_filename_fp,"plot '%s' using 1:2\n",gnuplot_filename);
245 fprintf(gnuplot_lin_command_filename_fp,"plot '%s' using 1:5\n",gnuplot_filename);
246 fprintf(gnuplot_lin_command_filename_fp,"pause -1 \"Hit RETURN to see pattern at %.4f MHz on a linear scale\"\n",f);
247 fprintf(gnuplot_log_command_filename_fp,"pause -1 \"Hit RETURN to see pattern at %.4f MHz on a log scale\"\n",f);
248 }
249 }
250 if(flag.pflg==1)
251 {
252 fclose(gnuplot_filename_fp);
253 fclose(gnuplot_log_command_filename_fp);
254 fclose(gnuplot_lin_command_filename_fp);
255 }
256 output_filename=string(0L, 100L);
257 original_filename=string(0L, 100L);
258 gain_filename=string(0L, 100L);
259 /* free strings not already freed */
260 free_string(input_filename,0L,100L);
261 free_string(temp,0L,100L);
262 free_string(gnuplot_filename,0L, 400L);
263 free_string(gnuplot_log_command_filename,0L, 400L);
264 free_string(gnuplot_lin_command_filename,0L, 400L);
265 free_string(output_filename, 0L, 100L);
266 free_string(original_filename, 0L,100L);
267 free_string(gain_filename, 0L, 100L);
268 /* free vectors */
269 free_FCOMPLEXvector(voltage,1L,(long) elements);
270 free_FCOMPLEXvector(current,1L,(long) elements);
271
272 fclose(ifp);
273 fclose(gain_fp);
274 exit(0);
275 }
276
277
278
error_3dB_E(double x)279 double error_3dB_E(double x)
280 {
281 double ans;
282
283 gain(x, 0, pin, normalised_f, coordinates, current, elements, &gain_E_plane, &gain_H_plane, f, design_f );
284 ans=peak_gain-gain_E_plane-3.01;
285 return(ans);
286 }
287
error_3dB_H(double x)288 double error_3dB_H(double x)
289 {
290 double ans;
291
292 gain(0, x, pin, normalised_f, coordinates, current, elements, &gain_E_plane, &gain_H_plane, f, design_f );
293 ans=peak_gain-gain_H_plane-3.01;
294 return(ans);
295 }
296