1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4
5 /* According to the method given in the book:
6 Yagi Antenna Design, by Dr. James L. Lawson, W2PV, ARRL.
7 Essentially, the self impedance Zxx of an elemnt is found. Then the mutural
8 impedance between two elements Zxy is found. Then the currents and voltages
9 are computed since:
10
11 I1 Z11 + I2 Z12 + I3 Z13 = V1
12 I1 Z21 + I2 Z22 + I3 Z23 = V2
13 I1 Z31 + I2 Z32 + I3 Z33 = V3
14
15 For just element 1 being driven, V2 = V3 = 0.
16 The input impedance is Zin= V_driven / I_driven
17 gain and other parameters can be found, as required. See book for
18 formulae */
19
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include <errno.h>
25 #include "yagi.h"
26
27 extern int optind, opterr, errno;
28
main(int argc,char ** argv)29 int main(int argc, char **argv)
30 {
31 char *output_filename, *input_filename, *line;
32 double **driven_data, **parasitic_data, **z, frequency, *v, **A, *b, *x;
33 double d=0.0, real=0.0, imaginary=0.0, min_frequency=0.0, max_frequency=0.0;
34 double step_frequency=0.0;
35 double angular_step=0.0;
36 int helpflg=0, errflg=0, pflg=0, dflg=0;
37 int elements=-1, driven=-1, parasitic=-1, i=0, j=0, *indx;
38 int element_number=0, c=0, sflg=0;
39 FILE *ofp;
40 struct FCOMPLEX *current, *voltage;
41 /* necessary with gcc under DOS */
42 opterr=0;
43
44 while ((c = getoptions(argc, argv, "dhps")) != -1)
45 switch (c)
46 {
47 case 'd': /* display currents */
48 dflg=1;
49 break;
50 case 'h': /* help */
51 helpflg=1;
52 break;
53 case 'p':
54 pflg=1;
55 break;
56 case 's': /* suppress all messages*/
57 sflg=1;
58 break;
59 case '?':
60 errflg=1;
61 break;
62 }
63
64 if(errflg || helpflg || (argc-optind)!=1)
65 {
66 usage_yagi(argv[0]);
67 exit(1);
68 }
69 /* allocate memory, using routines similar to that in Numerical Recipes */
70 input_filename=string(0L,90L);
71 line=string(0L,MAX_LINE);
72 output_filename=string(0L,90L);
73 output_filename=get_data_filenames(optind, argv, input_filename);
74 if(errno)
75 fprintf(stderr, "Errno = %d in yagi.c\n", errno);
76 /* since the size of all the arrays depends on the number of elements,
77 lets first find out how many elements we have,including driven and
78 parasitic sepparately. Then try to allocate all the memory we need. If we
79 fail, then not much time is wasted. */
80 elements=get_number_of_elements(input_filename, &driven, ¶sitic);
81 indx=ivector(1,2L*elements);
82
83 /* allocate all memory */
84 driven_data=dmatrix(1L,(long)(driven),1L,6L); /* 6 bits of info on each driven ele */
85 if(parasitic > 0)
86 parasitic_data=dmatrix(1L,(long) (parasitic),1L,4L); /* 4 on each parasitic ele */
87 /* the impedance matrix for the impedance Z11 to Znn, needs an NxN
88 matrix. Since the data is complex, it needs 2*N by N, where N is
89 the total number of elelemnts */
90 z=dmatrix(1L,(long) elements, 1L, elements*2L);
91 A=dmatrix(1L, 2L*elements, 1L, 2L*elements); /* A.x=b for large matrices */
92 x=dvector(1L, 2L*elements);
93 b=dvector(1L, 2L*elements);
94 /* read text file containing yagi description into memory */
95 read_yagi_data(line, input_filename, &frequency, &min_frequency, &max_frequency,
96 &step_frequency, driven, driven_data, parasitic, parasitic_data, &angular_step);
97 /* write some data to a disk file. A header of 100 items will be written
98 first. I doubt if I will ever need 100 items, but its better to have too
99 many, rather than too few. */
100 if(errno)
101 fprintf(stderr, "Errno = %d in yagi.c 2\n", errno);
102 ofp=fopen(output_filename,"wb");
103 errno=0;
104 write_header_to_disk(ofp, elements, driven, parasitic,min_frequency, max_frequency, frequency, step_frequency, angular_step);
105 /* write the x and y coordinates of the centre of the element, and the
106 length in m */
107
108 write_coordinates_of_elements_to_disk(ofp, driven, parasitic, driven_data, parasitic_data);
109 /* Now we have all the data in memory about the elements, we can calculate
110 the impedance matrices. This will be a square matrix, of size NxN, where
111 N is the number of elements in the yagi. The components on the diaganol
112 Z11, Z22, Z33 etc, will be the self impedance of each element. The
113 off-diagonal components, Zij (Z12, Z13 etc are the mutual impedance between
114 elements. */
115 current=FCOMPLEXvector(1L,(long) elements);
116 voltage=FCOMPLEXvector(1L,(long) elements);
117
118
119 v=dvector(1L,2L*elements);
120
121 for(frequency=min_frequency; frequency <=max_frequency; frequency+=step_frequency)
122 {
123 fill_z_matrix(frequency,driven,parasitic,driven_data,parasitic_data, z);
124 /* Now we must fill the V vector, which gives the voltage at the centre of
125 each driven element. Since we know the magnitude of the voltage and phase,
126 we can calculate the real and imaginary components of voltage. NB */
127 fill_v_vector(driven, parasitic, driven_data, v);
128 /* We now have the voltage vector V, and the impedance matrix Z. All
129 we need to do now is solve a set of NxN equations, where N is the
130 number of elements, to get N values of current, which we can put in
131 the I matrix.
132 Unfortunately, this is not trivual, by any standards!!! It is
133 complicated by the fact the the N equations are all complex.
134 The method used us a brute-force approach mentioned by Press
135 in their 2nd Edition of the Numerical Recipes in C book. Here we
136 put the data in the form:
137
138 |Zr -Zc| |Ic| = |Vr| which we will call A.x=b
139 |Zc Zr| |Ic| |Vc|
140
141 so the Z data goes now in a 2Nx2N matrix. This is a bit wasteful of space
142 and time, but it will do here. */
143 /* Read impedance data from 'z' to matrix A */
144 for(i=1;i<=elements;++i)
145 {
146 for(j=1;j<=elements;++j)
147 {
148 A[i][j]=z[i][2*j-1]; /* real data, top left corner of Z */
149 A[i+elements][j+elements]=z[i][2*j-1]; /* Bot right, real data */
150 A[i][j+elements]=-z[i][2*j]; /* imaginary, top right */
151 A[i+elements][j]=z[i][2*j]; /* bottom left, imaginary */
152 }
153 }
154 /* The following function prints to stantard output the z matrix of the
155 antenna. It is really only used during debugging, so it can be commented
156 out normally. */
157 if(pflg)
158 print_z_matrix(frequency,elements,z);
159 /* read voltage data from v to b */
160 for(i=1;i<=elements;++i)
161 {
162 b[i]= v[2*i-1]; /* real data */
163 b[i+elements]=v[2*i]; /* complex data */
164 voltage[i].r=v[2*i-1];
165 voltage[i].i=v[2*i];
166 /*
167 if(isnand(voltage[i].r) || isnand (voltage[i].i) || isnand(b[i]) || isnand(b[i+elements]) ||
168 (voltage[i].r) < 0.001)
169 {
170 fprintf(stderr, "b[%d]=%f b[%d]=%f %f %f\n", i, b[i], i+elements, b[i+elements], voltage[i].r, voltage[i].i);
171 }
172 */
173 }
174 /* we now write the voltage data to disk, but only once, *not* at
175 every frequency, since its fixed */
176 if(frequency==min_frequency)
177 {
178 fwrite((char *) &voltage[1], sizeof(struct FCOMPLEX), elements, ofp);
179 }
180 /* Perform a LU decompositon of A */
181 #ifdef DEBUG2
182 for(test=1;test <= 2*elements; ++test)
183 {
184 for( =1; <= 2*elements; ++)
185 printf("A[%d][%d]=%f\n", test, , A[test][]);
186 }
187 #endif
188 ludcmp(A, elements*2, indx, &d);
189 /* We now have the voltages in b. After lubksb is run, we get the
190 currents in b */
191
192 lubksb(A, 2*elements, indx, b); /* current's in b after lubksb has run*/
193 /* Put currents an FCOMPLEX matrix current[element]
194 currents are stored in b as r,r,r,r ..... i,i,i,i */
195 for(element_number=1;element_number<=elements;element_number++)
196 {
197 current[element_number].r=b[element_number];
198 current[element_number].i=b[element_number+elements];
199 }
200 if(dflg)
201 display_antenna_currents(current,elements);
202 fwrite((char *) ¤t[1], sizeof(struct FCOMPLEX), elements, ofp);
203 real=b[1]/((b[1]*b[1])+(b[1+elements]*b[1+elements]));
204 imaginary=-b[1+elements]/((b[1]*b[1])+(b[1+elements]*b[1+elements]));
205 /* printf("f=%.2f MHz z = 1/(%f + i %f) = %f + i%lf \n", frequency/1e6, b[1], b[1+elements], real, imaginary); */
206
207 if(min_frequency < max_frequency && !sflg)
208 {
209 printf("%s completed %5.1f%% f=%f MHz\n",argv[0],100*(frequency-min_frequency)/(max_frequency-min_frequency), frequency);
210 }
211 }
212 if(errno)
213 fprintf(stderr, "Errno = %d in yagi.c 3\n", errno);
214 printf("\n");
215 fclose(ofp);
216
217 /* free strings */
218 free_string(input_filename,0L,90L);
219 free_string(line,0L,(long) MAX_LINE);
220 free_string(output_filename,0L,90L);
221
222 /* free vectors */
223 free_dvector(x,1L, 2L*elements);
224 free_dvector(b,1L, 2L*elements);
225 free_dvector(v,1L, 2L*elements);
226 free_ivector(indx,1,2L*elements);
227 free_FCOMPLEXvector(current,1L,(long) elements);
228 free_FCOMPLEXvector(voltage,1L,(long) elements);
229 /* free matrices */
230 free_dmatrix(driven_data,1L,(long)(driven),1L,6L);
231 // free_dmatrix(parasitic_data,1L,(long) parasitic,1L,4L);
232 free_dmatrix(z,1L,(long) elements, 1L, elements*2L);
233 free_dmatrix(A,1L, 2L*elements, 1L, 2L*elements);
234 if(errno)
235 {
236 fprintf(stderr, "Errno = %d in yagi.c 4\n", errno);
237 exit (errno);
238 }
239 else
240 exit(0);
241 }
242