1 REM MININEC(3) ***************** NOSC CODE 822 (JCL CHANGE 12) 1-20-88
'
' INCLUDES THE ACUTE ANGLE KERNEL
'
' BASED ON MININEC(3) - CHANGE 11 OF 7-31-87
'
' NOTES: 1. MUST USE /MBF OPTION WITH QB87
'
2 DEFINT I-K, N
3 DIM K!(6, 2), Q(14)
REM ----- MAXIMUM NUMBER OF PULSES = MP
REM $DYNAMIC
MP = 50
4 DIM C%(MP, 2), CI(MP), CR(MP), P(MP), W%(MP)
5 DIM ZR(MP, MP), ZI(MP, MP)
REM ----- MAXIMUM NUMBER OF WIRES = MW
MW = 50
DIM A(MW), CA(MW), CB(MW), CG(MW), J1(MW), J2(MW, 2), N(MW, 2), S(MW)
REM ----- MAXIMUM NUMBER OF SEGMENTS (PULSES + 2 * WIRES) = MS
MS = MP + 2 * MW
6 DIM X(MS), Y(MS), Z(MS)
10 REM ----- MAXIMUM NUMBER OF LOADS = ML
11 ML = 12
12 REM ----- MAXIMUM ORDER OF S-DOMAIN LOADS = MA
13 MA = ML
14 DIM LA(2, ML, ML), LP(ML), LS(ML)
15 REM ----- MAXIMUM NUMBER OF MEDIA = 6
16 MM = 6
17 REM ----- H MUST BE DIMENSIONED AT LEAST = MM
18 DIM H(MM), T(MM), U(MM), V(MM), Z1(MM), Z2(MM)
23 REM ---- ARRAYS E,L & M DIMENSIONED TO MW+MP
24 DIM E(MW + MP), L(MW + MP), M(MW + MP)
25 COLOR 2, 0
26 GOTO 1497
27 REM ********** KERNEL EVALUATION OF INTEGRALS I2 & I3 **********
28 IF K < 0 THEN 33
29 X3 = X2 + T * (V1 - X2)
30 Y3 = Y2 + T * (V2 - Y2)
31 Z3 = Z2 + T * (V3 - Z2)
32 GOTO 36
33 X3 = V1 + T * (X2 - V1)
34 Y3 = V2 + T * (Y2 - V2)
35 Z3 = V3 + T * (Z2 - V3)
36 D3 = X3 * X3 + Y3 * Y3 + Z3 * Z3
37 REM ----- MOD FOR SMALL RADIUS TO WAVELENGTH RATIO
38 IF A(P4) <= SRM THEN D = SQR(D3): GOTO 49
39 D = D3 + A2
40 IF D > 0 THEN D = SQR(D)
41 REM ----- CRITERIA FOR USING REDUCED KERNEL
42 IF I6! = 0 THEN 49
REM ----- ADDITION FOR ACUTE ANGLE BETWEEN WIRES -- 1-88
IF IFLG = 0 THEN 44
REM ----- EXACT KERNEL FOR DIFFERENT WIRES
LL = X3 * CA(P4) + Y3 * CB(P4) + Z3 * CG(P4)
LL2 = LL * LL
RHO = SQR(D3 - LL2 + AM2)
D = SQR(LL2 + RHO * RHO)
RN = RHO - A(P4)
RP = RHO + A(P4)
RN2 = RN * RN
RP2 = RP * RP
B = (LL2 + RN2) / (LL2 + RP2)
TT = LOG((LL2 + RN2) / RP2 / 16) / RP
GOTO 45
43 REM ----- EXACT KERNEL CALCULATION WITH ELLIPTIC INTEGRAL
44 B = D3 / (D3 + 4 * A2)
REM ----- ADDITION FOR ACUTE ANGLE BETWEEN WIRES -- 1-88
TT = LOG(D3 / 64 / A2) / 2 / A(P4)
45 W0 = C0 + B * (C1 + B * (C2 + B * (C3 + B * C4)))
46 W1 = C5 + B * (C6 + B * (C7 + B * (C8 + B * C9)))
47 V0 = (W0 - W1 * LOG(B)) * SQR(1 - B)
REM ----- ADDITION FOR ACUTE ANGLE BETWEEN WIRES -- 1-88
IF IFLG = 1 THEN V0 = V0 / SQR(RHO * A(P4)) ELSE V0 = V0 / A(P4)
T3 = T3 + (V0 + TT) / P - 1 / D
49 B1 = D * W
50 REM ----- EXP(-J*K*R)/R
51 T3 = T3 + COS(B1) / D
52 T4 = T4 - SIN(B1) / D
53 RETURN
54 REM ***** PSI(P1,P2,P3) = T1 + J * T2 **********
55 REM ----- ENTRIES REQUIRED FOR NEAR FIELD CALCULATION
56 X1 = X0 + P1 * T5 / 2
57 Y1 = Y0 + P1 * T6 / 2
58 Z1 = Z0 + P1 * T7 / 2
59 X2 = X1 - X(P2)
60 Y2 = Y1 - Y(P2)
61 Z2 = Z1 - K * Z(P2)
62 V1 = X1 - X(P3)
63 V2 = Y1 - Y(P3)
64 V3 = Z1 - K * Z(P3)
65 GOTO 135
66 I4 = INT(P2)
67 I5 = I4 + 1
68 X2 = X0 - (X(I4) + X(I5)) / 2
69 Y2 = Y0 - (Y(I4) + Y(I5)) / 2
70 Z2 = Z0 - K * (Z(I4) + Z(I5)) / 2
71 V1 = X0 - X(P3)
72 V2 = Y0 - Y(P3)
73 V3 = Z0 - K * Z(P3)
74 GOTO 135
75 X2 = X0 - X(P2)
76 Y2 = Y0 - Y(P2)
77 Z2 = Z0 - K * Z(P2)
78 I4 = INT(P3)
79 I5 = I4 + 1
80 V1 = X0 - (X(I4) + X(I5)) / 2
81 V2 = Y0 - (Y(I4) + Y(I5)) / 2
82 V3 = Z0 - K * (Z(I4) + Z(I5)) / 2
83 GOTO 135
84 REM ----- ENTRIES REQUIRED FOR IMPEDANCE MATRIX CALCULATION
85 REM ----- S(M) GOES IN (X1,Y1,Z1) FOR SCALAR POTENTIAL
86 REM ----- MOD FOR SMALL RADIUS TO WAVE LENGTH RATIO
87 FVS = 1
88 IF K < 1 THEN 94
89 IF A(P4) > SRM THEN 94
90 IF (P3 = P2 + 1 AND P1 = (P2 + P3) / 2) THEN 91 ELSE 94
91 T1 = 2 * LOG(S(P4) / A(P4))
92 T2 = -W * S(P4)
93 RETURN
94 I4 = INT(P1)
95 I5 = I4 + 1
96 X1 = (X(I4) + X(I5)) / 2
97 Y1 = (Y(I4) + Y(I5)) / 2
98 Z1 = (Z(I4) + Z(I5)) / 2
99 GOTO 113
100 REM ----- S(M) GOES IN (X1,Y1,Z1) FOR VECTOR POTENTIAL
101 REM ----- MOD FOR SMALL RADIUS TO WAVE LENGTH RATIO
102 FVS = 0
103 IF K < 1 THEN 109
104 IF A(P4) >= SRM THEN 109
105 IF (I = J AND P3 = P2 + .5) THEN 106 ELSE 109
106 T1 = LOG(S(P4) / A(P4))
107 T2 = -W * S(P4) / 2
108 RETURN
109 X1 = X(P1)
110 Y1 = Y(P1)
111 Z1 = Z(P1)
112 REM ----- S(U)-S(M) GOES IN (X2,Y2,Z2)
113 I4 = INT(P2)
114 IF I4 = P2 THEN 120
115 I5 = I4 + 1
116 X2 = (X(I4) + X(I5)) / 2 - X1
117 Y2 = (Y(I4) + Y(I5)) / 2 - Y1
118 Z2 = K * (Z(I4) + Z(I5)) / 2 - Z1
119 GOTO 124
120 X2 = X(P2) - X1
121 Y2 = Y(P2) - Y1
122 Z2 = K * Z(P2) - Z1
123 REM ----- S(V)-S(M) GOES IN (V1,V2,V3)
124 I4 = INT(P3)
125 IF I4 = P3 THEN 131
126 I5 = I4 + 1
127 V1 = (X(I4) + X(I5)) / 2 - X1
128 V2 = (Y(I4) + Y(I5)) / 2 - Y1
129 V3 = K * (Z(I4) + Z(I5)) / 2 - Z1
130 GOTO 135
131 V1 = X(P3) - X1
132 V2 = Y(P3) - Y1
133 V3 = K * Z(P3) - Z1
134 REM ----- MAGNITUDE OF S(U) - S(M)
135 D0 = X2 * X2 + Y2 * Y2 + Z2 * Z2
136 REM ----- MAGNITUDE OF S(V) - S(M)
137 IF D0 > 0 THEN D0 = SQR(D0)
138 D3 = V1 * V1 + V2 * V2 + V3 * V3
139 IF D3 > 0 THEN D3 = SQR(D3)
140 REM ----- SQUARE OF WIRE RADIUS
141 A2 = A(P4) * A(P4)
142 REM ----- MAGNITUDE OF S(V) - S(U)
143 S4 = (P3 - P2) * S(P4)
144 REM ----- ORDER OF INTEGRATION
145 REM ----- LTH ORDER GAUSSIAN QUADRATURE
146 T1 = 0
147 T2 = 0
148 I6! = 0
149 F2 = 1
150 L = 7
151 T = (D0 + D3) / S(P4)
152 REM ----- CRITERIA FOR EXACT KERNEL
153 IF T > 1.1 THEN 165
REM ----- ADDITION FOR ACUTE ANGLE BETWEEN WIRES -- 1-88
IFLG = 0
F2 = 2 * (P3 - P2)
FVS1 = 1
IF FVS = 0 THEN FVS1 = 2
IF (D0 + D3) > ((S(P4) + A(P4)) / FVS1) THEN 154 ELSE 161
REM ----- I1 FOR DIFFERENT WIRES
154 LV = -(V1 * CA(P4) + V2 * CB(P4) + V3 * CG(P4))
LU = -(X2 * CA(P4) + Y2 * CB(P4) + Z2 * CG(P4))
RHOV = SQR(D3 * D3 - LV * LV + AM2)
RHOU = SQR(D0 * D0 - LU * LU + AM2)
RHO = RHOU
LL = LU
160 RN = RHO - A(P4)
RP = RHO + A(P4)
TF = 0
IF LL <> 0 THEN TF = LL * (LOG(LL * LL + RN * RN) / 2 - LOG(4 * RP) - 1)
X = 0
IF RN <> 0 THEN X = LL / RN
IF X < 100 THEN TT = RN * ATN(X) ELSE TT = RN * (P * SGN(LL) / 2 - RN / LL)
I6! = I6! + 2 * (TF + TT) / P / RP / S4
IF IFLG = 1 THEN 167
RHO = RHOV
LL = LV
I6! = -I6!
IFLG = 1
GOTO 160
REM ----- SELF TERM
161 IF A(P4) > SRM THEN 163
IF FVS = 1 THEN 91 ELSE 106
163 I6! = (1 - LOG(S4 / F2 / 8 / A(P4))) / P / A(P4)
164 GOTO 167
165 IF T > 6 THEN L = 3
166 IF T > 10 THEN L = 1
167 I5 = L + L
168 T3 = 0
169 T4 = 0
170 T = (Q(L) + .5) / F2
171 GOSUB 28
172 T = (.5 - Q(L)) / F2
173 GOSUB 28
174 L = L + 1
175 T1 = T1 + Q(L) * T3
176 T2 = T2 + Q(L) * T4
177 L = L + 1
178 IF L < I5 THEN 168
179 T1 = S4 * (T1 + I6!)
180 T2 = S4 * T2
181 RETURN
182 REM ********** COMPLEX SQUARE ROOT **********
183 REM ----- W6+I*W7=SQR(Z6+I*Z7)
184 T6 = SQR((ABS(Z6) + SQR(Z6 * Z6 + Z7 * Z7)) / 2)
185 T7 = ABS(Z7) / 2 / T6
186 IF Z6 < 0 THEN 191
187 W6 = T6
188 W7 = T7
189 IF Z7 < 0 THEN W7 = -T7
190 RETURN
191 W6 = T7
192 W7 = T6
193 IF Z7 < 0 THEN W7 = -T6
194 RETURN
195 REM ********** IMPEDANCE MATRIX CALCULATION **********
196 IF FLG = 1 THEN 428
197 IF FLG = 2 THEN 477
198 REM ----- BEGIN MATRIX FILL TIME CALCULATION
199 OT$ = TIME$
200 Q$ = "MATRIX FILL "
201 PRINT
202 PRINT "BEGIN "; Q$
203 REM ----- ZERO IMPEDANCE MATRIX
204 FOR I = 1 TO N
205 FOR J = 1 TO N
206 ZR(I, J) = 0
207 ZI(I, J) = 0
208 NEXT J
209 NEXT I
210 REM ----- COMPUTE ROW I OF MATRIX (OBSERVATION LOOP)
211 FOR I = 1 TO N
212 I1 = ABS(C%(I, 1))
213 I2 = ABS(C%(I, 2))
214 F4 = SGN(C%(I, 1)) * S(I1)
215 F5 = SGN(C%(I, 2)) * S(I2)
216 REM ----- R(M + 1/2) - R(M - 1/2) HAS COMPONENTS (T5,T6,T7)
217 T5 = F4 * CA(I1) + F5 * CA(I2)
218 T6 = F4 * CB(I1) + F5 * CB(I2)
219 T7 = F4 * CG(I1) + F5 * CG(I2)
220 IF C%(I, 1) = -C%(I, 2) THEN T7 = S(I1) * (CG(I1) + CG(I2))
REM ----- ADDITION FOR ACUTE ANGLE BETWEEN WIRES -- 1-88
REM ----- SQUARE OF LARGEST MATCH POINT RADIUS
AM2 = A(I1)
IF AM2 < A(I2) THEN AM2 = A(I2)
AM2 = AM2 * AM2
221 REM ----- COMPUTE COLUMN J OF ROW I (SOURCE LOOP)
222 FOR J = 1 TO N
223 J1 = ABS(C%(J, 1))
224 J2 = ABS(C%(J, 2))
225 F4 = SGN(C%(J, 1))
226 F5 = SGN(C%(J, 2))
227 F6 = 1
228 F7 = 1
229 REM ----- IMAGE LOOP
230 FOR K = 1 TO G STEP -2
231 IF C%(J, 1) <> -C%(J, 2) THEN 235
232 IF K < 0 THEN 332
233 F6 = F4
234 F7 = F5
235 F8 = 0
236 IF K < 0 THEN 248
237 REM ----- SET FLAG TO AVOID REDUNANT CALCULATIONS
238 IF I1 <> I2 THEN 246
239 IF (CA(I1) + CB(I1)) = 0 THEN 241
240 IF C%(I, 1) <> C%(I, 2) THEN 246
241 IF J1 <> J2 THEN 246
242 IF (CA(J1) + CB(J1)) = 0 THEN 244
243 IF C%(J, 1) <> C%(J, 2) THEN 246
244 IF I1 = J1 THEN F8 = 1
245 IF I = J THEN F8 = 2
246 IF ZR(I, J) <> 0 THEN 317
247 REM ----- COMPUTE PSI(M,N,N+1/2)
248 P1 = 2 * W%(I) + I - 1
249 P2 = 2 * W%(J) + J - 1
250 P3 = P2 + .5
251 P4 = J2
252 GOSUB 102
253 U1 = F5 * T1
254 U2 = F5 * T2
255 REM ----- COMPUTE PSI(M,N-1/2,N)
256 P3 = P2
257 P2 = P2 - .5
258 P4 = J1
259 IF F8 < 2 THEN GOSUB 102
260 V1 = F4 * T1
261 V2 = F4 * T2
262 REM ----- S(N+1/2)*PSI(M,N,N+1/2) + S(N-1/2)*PSI(M,N-1/2,N)
263 X3 = U1 * CA(J2) + V1 * CA(J1)
264 Y3 = U1 * CB(J2) + V1 * CB(J1)
265 Z3 = (F7 * U1 * CG(J2) + F6 * V1 * CG(J1)) * K
266 REM ----- REAL PART OF VECTOR POTENTIAL CONTRIBUTION
267 D1 = W2 * (X3 * T5 + Y3 * T6 + Z3 * T7)
268 X3 = U2 * CA(J2) + V2 * CA(J1)
269 Y3 = U2 * CB(J2) + V2 * CB(J1)
270 Z3 = (F7 * U2 * CG(J2) + F6 * V2 * CG(J1)) * K
271 REM ----- IMAGINARY PART OF VECTOR POTENTIAL CONTRIBUTION
272 D2 = W2 * (X3 * T5 + Y3 * T6 + Z3 * T7)
273 REM ----- COMPUTE PSI(M+1/2,N,N+1)
274 P1 = P1 + .5
275 IF F8 = 2 THEN P1 = P1 - 1
276 P2 = P3
277 P3 = P3 + 1
278 P4 = J2
279 IF F8 <> 1 THEN 283
280 U5 = F5 * U1 + T1
281 U6 = F5 * U2 + T2
282 GOTO 291
283 GOSUB 87
284 IF F8 < 2 THEN 288
285 U1 = (2 * T1 - 4 * U1 * F5) / S(J1)
286 U2 = (2 * T2 - 4 * U2 * F5) / S(J1)
287 GOTO 314
288 U5 = T1
289 U6 = T2
290 REM ----- COMPUTE PSI(M-1/2,N,N+1)
291 P1 = P1 - 1
292 GOSUB 87
293 U1 = (T1 - U5) / S(J2)
294 U2 = (T2 - U6) / S(J2)
295 REM ----- COMPUTE PSI(M+1/2,N-1,N)
296 P1 = P1 + 1
297 P3 = P2
298 P2 = P2 - 1
299 P4 = J1
300 GOSUB 87
301 U3 = T1
302 U4 = T2
303 REM ----- COMPUTE PSI(M-1/2,N-1,N)
304 IF F8 < 1 THEN 308
305 T1 = U5
306 T2 = U6
307 GOTO 311
308 P1 = P1 - 1
309 GOSUB 87
310 REM ----- GRADIENT OF SCALAR POTENTIAL CONTRIBUTION
311 U1 = U1 + (U3 - T1) / S(J1)
312 U2 = U2 + (U4 - T2) / S(J1)
313 REM ----- SUM INTO IMPEDANCE MATRIX
314 ZR(I, J) = ZR(I, J) + K * (D1 + U1)
315 ZI(I, J) = ZI(I, J) + K * (D2 + U2)
316 REM ----- AVOID REDUNANT CALCULATIONS
317 IF J < I THEN 332
318 IF F8 = 0 THEN 332
319 ZR(J, I) = ZR(I, J)
320 ZI(J, I) = ZI(I, J)
321 REM ----- SEGMENTS ON SAME WIRE SAME DISTANCE APART HAVE SAME Z
322 P1 = J + 1
323 IF P1 > N THEN 332
324 IF C%(P1, 1) <> C%(P1, 2) THEN 332
325 IF C%(P1, 2) = C%(J, 2) THEN 328
326 IF C%(P1, 2) <> -C%(J, 2) THEN 332
327 IF (CA(J2) + CB(J2)) <> 0 THEN 332
328 P2 = I + 1
329 IF P2 > N THEN 332
330 ZR(P2, P1) = ZR(I, J)
331 ZI(P2, P1) = ZI(I, J)
332 NEXT K
333 NEXT J
334 PCT = I / N
335 GOSUB 1599
336 NEXT I
337 REM ----- END MATRIX FILL TIME CALCULATION
338 T$ = TIME$
339 GOSUB 1589
340 PRINT #3, " "
341 PRINT #3, "FILL MATRIX : "; T$
342 REM ********** ADDITION OF LOADS **********
343 IF NL = 0 THEN 377
344 F5 = 2 * P * F * 1000000!
345 FOR I = 1 TO NL
346 IF LS(I) < 1 THEN 366
347 REM ----- S-DOMAIN LOADS
348 U1 = 0
349 U2 = 0
350 D1 = 0
351 D2 = 0
352 S = -1
353 FOR J = 0 TO LS(I) STEP 2
354 S = -S
355 U1 = U1 + LA(1, I, J) * S * F5 ^ J
356 D1 = D1 + LA(2, I, J) * S * F5 ^ J
357 L = J + 1
358 U2 = U2 + LA(1, I, L) * S * F5 ^ L
359 D2 = D2 + LA(2, I, L) * S * F5 ^ L
360 NEXT J
361 J = LP(I)
362 D = D1 * D1 + D2 * D2
363 LI = (U2 * D1 - D2 * U1) / D
364 LR = (U1 * D1 + U2 * D2) / D
365 GOTO 369
366 LR = LA(1, I, 1)
367 LI = LA(2, I, 1)
368 J = LP(I)
369 F2 = 1 / M
370 IF C%(J, 1) <> -C%(J, 2) THEN 372
371 IF K < 0 THEN F2 = 2 / M
372 ZR(J, J) = ZR(J, J) + F2 * LI
373 ZI(J, J) = ZI(J, J) - F2 * LR
374 NEXT I
375 REM ********** IMPEDANCE MATRIX FACTORIZATION **********
376 REM ----- BEGIN MATRIX FACTOR TIME CALCULATION
377 OT$ = TIME$
378 Q$ = "FACTOR MATRIX"
379 PRINT
380 PRINT "BEGIN "; Q$;
381 X = N
382 PCTN = X * (X - 1) * (X + X - 1)
383 FOR K = 1 TO N - 1
384 REM ----- SEARCH FOR PIVOT
385 T = ZR(K, K) * ZR(K, K) + ZI(K, K) * ZI(K, K)
386 I1 = K
387 FOR I = K + 1 TO N
388 T1 = ZR(I, K) * ZR(I, K) + ZI(I, K) * ZI(I, K)
389 IF T1 < T THEN 392
390 I1 = I
391 T = T1
392 NEXT I
393 REM ----- EXCHANGE ROWS K AND I1
394 IF I1 = K THEN 403
395 FOR J = 1 TO N
396 T1 = ZR(K, J)
397 T2 = ZI(K, J)
398 ZR(K, J) = ZR(I1, J)
399 ZI(K, J) = ZI(I1, J)
400 ZR(I1, J) = T1
401 ZI(I1, J) = T2
402 NEXT J
403 P(K) = I1
404 REM ----- SUBTRACT ROW K FROM ROWS K+1 TO N
405 FOR I = K + 1 TO N
406 REM ----- COMPUTE MULTIPLIER L(I,K)
407 T1 = (ZR(I, K) * ZR(K, K) + ZI(I, K) * ZI(K, K)) / T
408 T2 = (ZI(I, K) * ZR(K, K) - ZR(I, K) * ZI(K, K)) / T
409 ZR(I, K) = T1
410 ZI(I, K) = T2
411 REM ----- SUBTRACT ROW K FROM ROW I
412 FOR J = K + 1 TO N
413 ZR(I, J) = ZR(I, J) - (ZR(K, J) * T1 - ZI(K, J) * T2)
414 ZI(I, J) = ZI(I, J) - (ZR(K, J) * T2 + ZI(K, J) * T1)
415 NEXT J
416 NEXT I
417 X = N - K
418 PCT = 1 - X * (X - 1) * (X + X - 1) / PCTN
419 GOSUB 1599
420 NEXT K
421 REM ----- END MATRIX FACTOR TIME CALCULATION
422 T$ = TIME$
423 GOSUB 1589
424 PRINT
425 PRINT #3, "FACTOR MATRIX: "; T$
426 REM ********** SOLVE **********
427 REM ----- COMPUTE RIGHT HAND SIDE
428 FOR I = 1 TO N
429 CR(I) = 0
430 CI(I) = 0
431 NEXT I
432 FOR J = 1 TO NS
433 F2 = 1 / M
434 IF C%(E(J), 1) = -C%(E(J), 2) THEN F2 = 2 / M
435 CR(E(J)) = F2 * M(J)
436 CI(E(J)) = -F2 * L(J)
437 NEXT J
438 REM ----- PERMUTE EXCITATION
439 FOR K = 1 TO N - 1
440 I1 = P(K)
441 IF I1 = K THEN 448
442 T1 = CR(K)
443 T2 = CI(K)
444 CR(K) = CR(I1)
445 CI(K) = CI(I1)
446 CR(I1) = T1
447 CI(I1) = T2
448 NEXT K
449 REM ----- FORWARD ELIMINATION
450 FOR I = 2 TO N
451 T1 = 0
452 T2 = 0
453 FOR J = 1 TO I - 1
454 T1 = T1 + ZR(I, J) * CR(J) - ZI(I, J) * CI(J)
455 T2 = T2 + ZR(I, J) * CI(J) + ZI(I, J) * CR(J)
456 NEXT J
457 CR(I) = CR(I) - T1
458 CI(I) = CI(I) - T2
459 NEXT I
460 REM ----- BACK SUBSTITUTION
461 FOR I = N TO 1 STEP -1
462 T1 = 0
463 T2 = 0
464 IF I = N THEN 469
465 FOR J = I + 1 TO N
466 T1 = T1 + ZR(I, J) * CR(J) - ZI(I, J) * CI(J)
467 T2 = T2 + ZR(I, J) * CI(J) + ZI(I, J) * CR(J)
468 NEXT J
469 T = ZR(I, I) * ZR(I, I) + ZI(I, I) * ZI(I, I)
470 T1 = CR(I) - T1
471 T2 = CI(I) - T2
472 CR(I) = (T1 * ZR(I, I) + T2 * ZI(I, I)) / T
473 CI(I) = (T2 * ZR(I, I) - T1 * ZI(I, I)) / T
474 NEXT I
475 FLG = 2
476 REM ********** SOURCE DATA **********
477 PRINT #3, " "
478 PRINT #3, B$; " SOURCE DATA "; B$
479 PWR = 0
480 FOR I = 1 TO NS
481 CR = CR(E(I))
482 CI = CI(E(I))
483 T = CR * CR + CI * CI
484 T1 = (L(I) * CR + M(I) * CI) / T
485 T2 = (M(I) * CR - L(I) * CI) / T
486 O2 = (L(I) * CR + M(I) * CI) / 2
487 PWR = PWR + O2
488 PRINT #3, "PULSE "; E(I), "VOLTAGE = ("; L(I); ","; M(I); "J)"
489 PRINT #3, " ", "CURRENT = ("; CR; ","; CI; "J)"
490 PRINT #3, " ", "IMPEDANCE = ("; T1; ","; T2; "J)"
491 PRINT #3, " ", "POWER = "; O2; " WATTS"
492 NEXT I
493 IF NS > 1 THEN PRINT #3, " "
494 IF NS > 1 THEN PRINT #3, "TOTAL POWER = "; PWR; "WATTS"
495 RETURN
496 REM ********** PRINT CURRENTS **********
497 GOSUB 196
498 S$ = "N"
499 PRINT #3, " "
500 PRINT #3, B$; " CURRENT DATA "; B$
501 FOR K = 1 TO NW
502 IF S$ = "Y" THEN 507
503 PRINT #3, " "
504 PRINT #3, "WIRE NO. "; K; ":"
505 PRINT #3, "PULSE", "REAL", "IMAGINARY", "MAGNITUDE", "PHASE"
506 PRINT #3, " NO.", "(AMPS)", "(AMPS)", "(AMPS)", "(DEGREES)"
507 N1 = N(K, 1)
508 N2 = N(K, 2)
509 I = N1
510 C = C%(I, 1)
511 IF (N1 = 0 AND N2 = 0) THEN C = K
512 IF G = 1 THEN 515
513 IF (J1(K) = -1 AND N1 > N2) THEN N2 = N1
514 IF J1(K) = -1 THEN 525
515 E% = 1
516 GOSUB 572
517 I2! = I1!
518 J2! = J1!
519 GOSUB 607
520 IF S$ = "N" THEN PRINT #3, I$, I1!; TAB(29); J1!; TAB(43); S1; TAB(57); S2
521 IF S$ = "Y" THEN PRINT #1, I1!; ","; J1!; ","; S1; ","; S2
522 IF N1 = 0 THEN 532
523 IF C = K THEN 525
524 IF I$ = "J" THEN N1 = N1 + 1
525 FOR I = N1 TO N2 - 1
526 I2! = CR(I)
527 J2! = CI(I)
528 GOSUB 607
529 IF S$ = "N" THEN PRINT #3, I, CR(I); TAB(29); CI(I); TAB(43); S1; TAB(57); S2
530 IF S$ = "Y" THEN PRINT #1, CR(I); ","; CI(I); ","; S1; ","; S2
531 NEXT I
532 I = N2
533 C = C%(I, 2)
534 IF (N1 = 0 AND N2 = 0) THEN C = K
535 IF G = 1 THEN 537
536 IF J1(K) = 1 THEN 543
537 E% = 2
538 GOSUB 572
539 IF (N1 = 0 AND N2 = 0) THEN 549
540 IF N1 > N2 THEN 549
541 IF C = K THEN 543
542 IF I$ = "J" THEN 549
543 I2! = CR(N2)
544 J2! = CI(N2)
545 GOSUB 607
546 IF S$ = "N" THEN PRINT #3, N2, CR(N2); TAB(29); CI(N2); TAB(43); S1; TAB(57); S2
547 IF S$ = "Y" THEN PRINT #1, CR(N2); ","; CI(N2); ","; S1; ","; S2
548 IF J1(K) = 1 THEN 554
549 I2! = I1!
550 J2! = J1!
551 GOSUB 607
552 IF S$ = "N" THEN PRINT #3, I$, I1!; TAB(29); J1!; TAB(43); S1; TAB(57); S2
553 IF S$ = "Y" THEN PRINT #1, I1!; ","; J1!; ","; S1; ","; S2
554 IF S$ = "Y" THEN PRINT #1, " 1 , 1 , 1 , 1"
555 NEXT K
556 IF S$ = "Y" THEN 569
557 PRINT
558 INPUT "SAVE CURRENTS TO A FILE (Y/N) "; S$
559 IF S$ = "N" THEN 569
560 IF S$ <> "Y" THEN 557
561 PRINT #3, " "
562 INPUT "FILENAME (NAME.OUT) "; F$
563 IF LEFT$(RIGHT$(F$, 4), 1) = "." THEN 564 ELSE F$ = F$ + ".OUT"
564 IF O$ > "C" THEN PRINT #3, "FILENAME (NAME.OUT): "; F$
565 OPEN F$ FOR OUTPUT AS #1
566 PRINT #3, " "
567 PRINT #1, NW; ","; PWR; ",C"
568 GOTO 501
569 CLOSE #1
570 RETURN
571 REM ----- SORT JUNCTION CURRENTS
572 I$ = "E"
573 I1! = 0!
574 J1! = 0!
575 IF (C = K OR C = 0) THEN 580
576 I$ = "J"
577 I1! = CR(I)
578 J1! = CI(I)
579 REM ----- CHECK FOR OTHER OVERLAPPING WIRES
580 FOR J = 1 TO NW
581 IF J = K GOTO 604
582 L1 = N(J, 1)
583 L2 = N(J, 2)
584 IF E% = 2 THEN 590
585 CO = C%(L1, 1)
586 CT = C%(L2, 2)
587 L3 = L1
588 L4 = L2
589 GOTO 594
590 CO = C%(L2, 2)
591 CT = C%(L1, 1)
592 L3 = L2
593 L4 = L1
594 IF CO = -K THEN 596
595 GOTO 599
596 I1! = I1! - CR(L3)
597 J1! = J1! - CI(L3)
598 I$ = "J"
599 IF CT = K THEN 601
600 GOTO 604
601 I1! = I1! + CR(L4)
602 J1! = J1! + CI(L4)
603 I$ = "J"
604 NEXT J
605 RETURN
606 REM ----- CALCULATE S1 AND S2
607 I3! = I2! * I2!
608 J3! = J2! * J2!
609 IF (I3! > 0 OR J3! > 0) THEN 612
610 S1 = 0!
611 GOTO 613
612 S1 = SQR(I3! + J3!)
613 IF I2! <> 0 THEN 616
614 S2 = 0!
615 RETURN
616 S2 = ATN(J2! / I2!) / P0
617 IF I2! > 0 THEN RETURN
618 S2 = S2 + SGN(J2!) * 180
619 RETURN
620 REM ********** FAR FIELD CALCULATION **********
621 IF FLG < 2 THEN GOSUB 196
622 O2 = PWR
623 REM ----- TABULATE IMPEDANCE
624 IF NM = 0 THEN 634
625 FOR I = 1 TO NM
626 Z6 = T(I)
627 Z7 = -V(I) / (2 * P * F * 8.85E-06)
628 REM ----- FORM IMPEDANCE=1/SQR(DIELECTRIC CONSTANT)
629 GOSUB 184
630 D = W6 * W6 + W7 * W7
631 Z1(I) = W6 / D
632 Z2(I) = -W7 / D
633 NEXT I
634 PRINT #3, " "
635 PRINT #3, B$; " FAR FIELD "; B$
636 PRINT #3, " "
637 REM ----- INPUT VARIABLES FOR FAR FIELD CALCULATION
638 INPUT "CALCULATE PATTERN IN DBI OR VOLTS/METER (D/V)"; P$
639 IF P$ = "D" THEN 655
640 IF P$ <> "V" THEN 638
641 F1 = 1
642 PRINT
643 PRINT "PRESENT POWER LEVEL = "; PWR; " WATTS"
644 INPUT "CHANGE POWER LEVEL (Y/N) "; A$
645 IF A$ = "N" THEN 650
646 IF A$ <> "Y" THEN 644
647 INPUT "NEW POWER LEVEL (WATTS) "; O2
648 IF O$ > "C" THEN PRINT #3, "NEW POWER LEVEL = "; O2
649 GOTO 644
650 IF (O2 < 0 OR O2 = 0) THEN O2 = PWR
651 F1 = SQR(O2 / PWR)
652 PRINT
653 INPUT "RADIAL DISTANCE (METERS) "; RD
654 IF RD < 0 THEN RD = 0
655 A$ = "ZENITH ANGLE : INITIAL,INCREMENT,NUMBER"
656 PRINT A$;
657 INPUT ZA, ZC, NZ
658 IF NZ = 0 THEN NZ = 1
659 IF O$ > "C" THEN PRINT #3, A$; ": "; ZA; ","; ZC; ","; NZ
660 A$ = "AZIMUTH ANGLE: INITIAL,INCREMENT,NUMBER"
661 PRINT A$;
662 INPUT AA, AC, NA
663 IF NA = 0 THEN NA = 1
664 IF O$ > "C" THEN PRINT #3, A$; ": "; AA; ","; AC; ","; NA
665 PRINT #3, " "
666 REM ********** FILE FAR FIELD DATA **********
667 INPUT "FILE PATTERN (Y/N)"; S$
668 IF S$ = "N" THEN 676
669 IF S$ <> "Y" THEN 667
670 PRINT #3, " "
671 INPUT "FILENAME (NAME.OUT)"; F$
672 IF LEFT$(RIGHT$(F$, 4), 1) = "." THEN 673 ELSE F$ = F$ + ".OUT"
673 IF O$ > "C" THEN PRINT #3, "FILENAME (NAME.OUT): "; F$
674 OPEN F$ FOR OUTPUT AS #1
675 PRINT #1, NA * NZ; ","; O2; ","; P$
676 PRINT #3, " "
677 K9! = .016678 / PWR
678 REM ----- PATTERN HEADER
679 PRINT #3, B$; " PATTERN DATA "; B$
680 IF P$ = "V" GOTO 685
681 PRINT #3, "ZENITH", "AZIMUTH", "VERTICAL", "HORIZONTAL", "TOTAL"
682 A$ = "PATTERN (DB)"
683 PRINT #3, " ANGLE", " ANGLE", A$, A$, A$
684 GOTO 692
685 IF RD > 0 THEN PRINT #3, TAB(15); "RADIAL DISTANCE = "; RD; " METERS"
686 PRINT #3, TAB(15); "POWER LEVEL = "; PWR * F1 * F1; " WATTS"
687 PRINT #3, "ZENITH AZIMUTH", " E(THETA) ", " E(PHI)"
688 A$ = " MAG(V/M) PHASE(DEG)"
689 PRINT #3, " ANGLE ANGLE", A$, A$
690 IF S$ = "Y" THEN PRINT #1, RD
691 REM ----- LOOP OVER AZIMUTH ANGLE
692 Q1 = AA
693 FOR I1 = 1 TO NA
694 U3 = Q1 * P0
695 V1 = -SIN(U3)
696 V2 = COS(U3)
697 REM ----- LOOP OVER ZENITH ANGLE
698 Q2 = ZA
699 FOR I2 = 1 TO NZ
700 U4 = Q2 * P0
701 R3 = COS(U4)
702 T3 = -SIN(U4)
703 T1 = R3 * V2
704 T2 = -R3 * V1
705 R1 = -T3 * V2
706 R2 = T3 * V1
707 X1 = 0
708 Y1 = 0
709 Z1 = 0
710 X2 = 0
711 Y2 = 0
712 Z2 = 0
713 REM ----- IMAGE LOOP
714 FOR K = 1 TO G STEP -2
715 FOR I = 1 TO N
716 IF K > 0 THEN 718
717 IF C%(I, 1) = -C%(I, 2) THEN 813
718 J = 2 * W%(I) - 1 + I
719 REM ----- FOR EACH END OF PULSE COMPUTE A CONTRIBUTION TO E-FIELD
720 FOR F5 = 1 TO 2
721 L = ABS(C%(I, F5))
722 F3 = SGN(C%(I, F5)) * W * S(L) / 2
723 IF C%(I, 1) <> -C%(I, 2) THEN 725
724 IF F3 < 0 THEN 812
725 IF K = 1 THEN 728
726 IF NM <> 0 THEN 747
727 REM ----- STANDARD CASE
728 S2 = W * (X(J) * R1 + Y(J) * R2 + Z(J) * K * R3)
729 S1 = COS(S2)
730 S2 = SIN(S2)
731 B1 = F3 * (S1 * CR(I) - S2 * CI(I))
732 B2 = F3 * (S1 * CI(I) + S2 * CR(I))
733 IF C%(I, 1) = -C%(I, 2) THEN 742
734 X1 = X1 + K * B1 * CA(L)
735 X2 = X2 + K * B2 * CA(L)
736 Y1 = Y1 + K * B1 * CB(L)
737 Y2 = Y2 + K * B2 * CB(L)
738 Z1 = Z1 + B1 * CG(L)
739 Z2 = Z2 + B2 * CG(L)
740 GOTO 812
741 REM ----- GROUNDED ENDS
742 Z1 = Z1 + 2 * B1 * CG(L)
743 Z2 = Z2 + 2 * B2 * CG(L)
744 GOTO 812
745 REM ----- REAL GROUND CASE
746 REM ----- BEGIN BY FINDING SPECULAR DISTANCE
747 T4 = 100000!
748 IF R3 = 0 THEN 750
749 T4 = -Z(J) * T3 / R3
750 B9 = T4 * V2 + X(J)
751 IF TB = 1 THEN 755
752 B9 = B9 * B9 + (Y(J) - T4 * V1) ^ 2
753 IF B9 > 0 THEN B9 = SQR(B9) ELSE 755
754 REM ----- SEARCH FOR THE CORRESPONDING MEDIUM
755 J2 = NM
756 FOR J1 = NM TO 1 STEP -1
757 IF B9 > U(J1) THEN 759
758 J2 = J1
759 NEXT J1
760 REM ----- OBTAIN IMPEDANCE AT SPECULAR POINT
761 Z4 = Z1(J2)
762 Z5 = Z2(J2)
763 REM ----- IF PRESENT INCLUDE GROUND SCREEN IMPEDANCE IN PARALLEL
764 IF NR = 0 THEN 776
765 IF B9 > U(1) THEN 776
766 R = B9 + NR * RR
767 Z8 = W * R * LOG(R / (NR * RR)) / NR
768 S8 = -Z5 * Z8
769 S9 = Z4 * Z8
770 T8 = Z4
771 T9 = Z5 + Z8
772 D = T8 * T8 + T9 * T9
773 Z4 = (S8 * T8 + S9 * T9) / D
774 Z5 = (S9 * T8 - S8 * T9) / D
775 REM ----- FORM SQR(1-Z^2*SIN^2)
776 Z6 = 1 - (Z4 * Z4 - Z5 * Z5) * T3 * T3
777 Z7 = -(2 * Z4 * Z5) * T3 * T3
778 GOSUB 184
779 REM ----- VERTICAL REFLECTION COEFFICIENT
780 S8 = R3 - (W6 * Z4 - W7 * Z5)
781 S9 = -(W6 * Z5 + W7 * Z4)
782 T8 = R3 + (W6 * Z4 - W7 * Z5)
783 T9 = W6 * Z5 + W7 * Z4
784 D = T8 * T8 + T9 * T9
785 V8 = (S8 * T8 + S9 * T9) / D
786 V9 = (S9 * T8 - S8 * T9) / D
787 REM ----- HORIZONTAL REFLECTION COEFFICIENT
788 S8 = W6 - R3 * Z4
789 S9 = W7 - R3 * Z5
790 T8 = W6 + R3 * Z4
791 T9 = W7 + R3 * Z5
792 D = T8 * T8 + T9 * T9
793 H8 = (S8 * T8 + S9 * T9) / D - V8
794 H9 = (S9 * T8 - S8 * T9) / D - V9
795 REM ----- COMPUTE CONTRIBUTION TO SUM
796 S2 = W * (X(J) * R1 + Y(J) * R2 - (Z(J) - 2 * H(J2)) * R3)
797 S1 = COS(S2)
798 S2 = SIN(S2)
799 B1 = F3 * (S1 * CR(I) - S2 * CI(I))
800 B2 = F3 * (S1 * CI(I) + S2 * CR(I))
801 W6 = B1 * V8 - B2 * V9
802 W7 = B1 * V9 + B2 * V8
803 D = CA(L) * V1 + CB(L) * V2
804 Z6 = D * (B1 * H8 - B2 * H9)
805 Z7 = D * (B1 * H9 + B2 * H8)
806 X1 = X1 - (CA(L) * W6 + V1 * Z6)
807 X2 = X2 - (CA(L) * W7 + V1 * Z7)
808 Y1 = Y1 - (CB(L) * W6 + V2 * Z6)
809 Y2 = Y2 - (CB(L) * W7 + V2 * Z7)
810 Z1 = Z1 + CG(L) * W6
811 Z2 = Z2 + CG(L) * W7
812 NEXT F5
813 NEXT I
814 NEXT K
815 H2 = -(X1 * T1 + Y1 * T2 + Z1 * T3) * G0
816 H1 = (X2 * T1 + Y2 * T2 + Z2 * T3) * G0
817 X4 = -(X1 * V1 + Y1 * V2) * G0
818 X3 = (X2 * V1 + Y2 * V2) * G0
819 IF P$ = "D" THEN 827
820 IF RD = 0 THEN 842
821 H1 = H1 / RD
822 H2 = H2 / RD
823 X3 = X3 / RD
824 X4 = X4 / RD
825 GOTO 842
826 REM ----- PATTERN IN DB
827 P1 = -999
828 P2 = P1
829 P3 = P1
830 T1 = K9! * (H1 * H1 + H2 * H2)
831 T2 = K9! * (X3 * X3 + X4 * X4)
832 T3 = T1 + T2
833 REM ----- CALCULATE VALUES IN DB
834 IF T1 > 1E-30 THEN P1 = 4.343 * LOG(T1)
835 IF T2 > 1E-30 THEN P2 = 4.343 * LOG(T2)
836 IF T3 > 1E-30 THEN P3 = 4.343 * LOG(T3)
837 PRINT #3, Q2; TAB(15); Q1; TAB(29); P1; TAB(43); P2; TAB(57); P3
838 IF S$ = "Y" THEN PRINT #1, Q2; ","; Q1; ","; P1; ","; P2; ","; P3
839 GOTO 866
840 REM ----- PATTERN IN VOLTS/METER
841 REM ----- MAGNITUDE AND PHASE OF E(THETA)
842 S1 = 0
843 IF (H1 = 0 AND H2 = 0) THEN 845
844 S1 = SQR(H1 * H1 + H2 * H2)
845 IF H1 <> 0 THEN 848
846 S2 = 0
847 GOTO 851
848 S2 = ATN(H2 / H1) / P0
849 IF H1 < 0 THEN S2 = S2 + SGN(H2) * 180
850 REM ----- MAGNITUDE AND PHASE OF E(PHI)
851 S3 = 0
852 IF (X3 = 0 AND X4 = 0) THEN 854
853 S3 = SQR(X3 * X3 + X4 * X4)
854 IF X3 <> 0 THEN 857
855 S4 = 0
856 GOTO 859
857 S4 = ATN(X4 / X3) / P0
858 IF X3 < 0 THEN S4 = S4 + SGN(X4) * 180
859 PRINT #3, USING "###.## "; Q2, Q1;
860 PRINT #3, USING " ##.###^^^^"; S1 * F1;
861 PRINT #3, USING " ###.## "; S2;
862 PRINT #3, USING " ##.###^^^^"; S3 * F1;
863 PRINT #3, USING " ###.##"; S4
864 IF S$ = "Y" THEN PRINT #1, Q2; ","; Q1; ","; S1 * F1; ","; S2; ","; S3 * F1; ","; S4
865 REM ----- INCREMENT ZENITH ANGLE
866 Q2 = Q2 + ZC
867 NEXT I2
868 REM ----- INCREMENT AZIMUTH ANGLE
869 Q1 = Q1 + AC
870 NEXT I1
871 CLOSE #1
872 RETURN
873 REM ********** NEAR FIELD CALCULATION **********
874 REM ----- ENSURE CURRENTS HAVE BEEN CALCULATED
875 IF FLG < 2 THEN GOSUB 196
876 FVS = -1: O2 = PWR
877 PRINT #3, " "
878 PRINT #3, B$; " NEAR FIELDS "; B$
879 PRINT #3, " "
880 INPUT "ELECTRIC OR MAGNETIC NEAR FIELDS (E/H) "; N$
881 IF (N$ = "H" OR N$ = "E") GOTO 883
882 GOTO 880
883 PRINT
884 REM ----- INPUT VARIABLES FOR NEAR FIELD CALCULATION
885 PRINT "FIELD LOCATION(S):"
886 A$ = "-COORDINATE (M): INITIAL,INCREMENT,NUMBER "
887 PRINT " X"; A$;
888 INPUT XI, XC, NX
889 IF NX = 0 THEN NX = 1
890 IF O$ > "C" THEN PRINT #3, "X"; A$; ": "; XI; ","; XC; ","; NX
891 PRINT " Y"; A$;
892 INPUT YI, YC, NY
893 IF NY = 0 THEN NY = 1
894 IF O$ > "C" THEN PRINT #3, "Y"; A$; ": "; YI; ","; YC; ","; NY
895 PRINT " Z"; A$;
896 INPUT ZI, ZC, NZ
897 IF NZ = 0 THEN NZ = 1
898 IF O$ > "C" THEN PRINT #3, "Z"; A$; ": "; ZI; ","; ZC; ","; NZ
899 F1 = 1
900 PRINT
901 PRINT "PRESENT POWER LEVEL IS "; PWR; " WATTS"
902 INPUT "CHANGE POWER LEVEL (Y/N) "; A$
903 IF A$ = "N" THEN 908
904 IF A$ <> "Y" THEN 902
905 INPUT "NEW POWER LEVEL (WATTS) "; O2
906 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, "NEW POWER LEVEL (WATTS) = "; O2
907 GOTO 902
908 IF (O2 < 0 OR O2 = 0) THEN O2 = PWR
909 REM ----- RATIO OF POWER LEVELS
910 F1 = SQR(O2 / PWR)
911 IF N$ = "H" THEN F1 = F1 / S0 / 4 / P
912 PRINT
913 REM ----- DESIGNATION OF OUTPUT FILE FOR NEAR FIELD DATA
914 INPUT "SAVE TO A FILE (Y/N) "; S$
915 IF S$ = "N" THEN 923
916 IF S$ <> "Y" THEN 914
917 INPUT "FILENAME (NAME.OUT) "; F$
918 IF LEFT$(RIGHT$(F$, 4), 1) = "." THEN 919 ELSE F$ = F$ + ".OUT"
919 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, "FILENAME (NAME.OUT) "; F$
920 OPEN F$ FOR OUTPUT AS #2
921 PRINT #2, NX * NY * NZ; ","; O2; ","; N$
922 REM ----- LOOP OVER Z DIMENSION
923 FOR IZ = 1 TO NZ
924 ZZ = ZI + (IZ - 1) * ZC
925 REM ----- LOOP OVER Y DIMENSION
926 FOR IY = 1 TO NY
927 YY = YI + (IY - 1) * YC
928 REM ----- LOOP OVER X DIMENSION
929 FOR IX = 1 TO NX
930 XX = XI + (IX - 1) * XC
931 REM ----- NEAR FIELD HEADER
932 PRINT #3, " "
933 IF N$ = "E" THEN PRINT #3, B$; "NEAR ELECTRIC FIELDS"; B$
934 IF N$ = "H" THEN PRINT #3, B$; "NEAR MAGNETIC FIELDS"; B$
935 PRINT #3, TAB(10); "FIELD POINT: "; "X = "; XX; " Y = "; YY; " Z = "; ZZ
936 PRINT #3, " VECTOR", "REAL", "IMAGINARY", "MAGNITUDE", "PHASE"
937 IF N$ = "E" THEN A$ = " V/M "
938 IF N$ = "H" THEN A$ = " AMPS/M "
939 PRINT #3, " COMPONENT ", A$, A$, A$, " DEG"
940 A1 = 0
941 A3 = 0
942 A4 = 0
943 REM ----- LOOP OVER THREE VECTOR COMPONENTS
944 FOR I = 1 TO 3
945 X0 = XX
946 Y0 = YY
947 Z0 = ZZ
948 IF N$ = "H" THEN 958
949 T5 = 0
950 T6 = 0
951 T7 = 0
952 IF I = 1 THEN T5 = 2 * S0
953 IF I = 2 THEN T6 = 2 * S0
954 IF I = 3 THEN T7 = 2 * S0
955 U7 = 0
956 U8 = 0
957 GOTO 968
958 FOR J8 = 1 TO 6
959 K!(J8, 1) = 0
960 K!(J8, 2) = 0
961 NEXT J8
962 J9 = 1
963 J8 = -1
964 IF I = 1 THEN X0 = XX + J8 * S0 / 2
965 IF I = 2 THEN Y0 = YY + J8 * S0 / 2
966 IF I = 3 THEN Z0 = ZZ + J8 * S0 / 2
967 REM ----- LOOP OVER SOURCE SEGMENTS
968 FOR J = 1 TO N
969 J1 = ABS(C%(J, 1))
970 J2 = ABS(C%(J, 2))
971 J3 = J2
972 IF J1 > J2 THEN J3 = J1
973 F4 = SGN(C%(J, 1))
974 F5 = SGN(C%(J, 2))
975 F6 = 1
976 F7 = 1
977 U5 = 0
978 U6 = 0
979 REM ----- IMAGE LOOP
980 FOR K = 1 TO G STEP -2
981 IF C%(J, 1) <> -C%(J, 2) THEN 987
982 IF K < 0 THEN 1048
983 REM ----- COMPUTE VECTOR POTENTIAL A
984 F6 = F4
985 F7 = F5
986 REM ----- COMPUTE PSI(0,J,J+.5)
987 P1 = 0
988 P2 = 2 * J3 + J - 1
989 P3 = P2 + .5
990 P4 = J2
991 GOSUB 75
992 U1 = T1 * F5
993 U2 = T2 * F5
994 REM ----- COMPUTE PSI(0,J-.5,J)
995 P3 = P2
996 P2 = P2 - .5
997 P4 = J1
998 GOSUB 66
999 V1 = F4 * T1
1000 V2 = F4 * T2
1001 REM ----- REAL PART OF VECTOR POTENTIAL CONTRIBUTION
1002 X3 = U1 * CA(J2) + V1 * CA(J1)
1003 Y3 = U1 * CB(J2) + V1 * CB(J1)
1004 Z3 = (F7 * U1 * CG(J2) + F6 * V1 * CG(J1)) * K
1005 REM ----- IMAGINARY PART OF VECTOR POTENTIAL CONTRIBUTION
1006 X5 = U2 * CA(J2) + V2 * CA(J1)
1007 Y5 = U2 * CB(J2) + V2 * CB(J1)
1008 Z5 = (F7 * U2 * CG(J2) + F6 * V2 * CG(J1)) * K
1009 REM ----- MAGNETIC FIELD CALCULATION COMPLETED
1010 IF N$ = "H" THEN 1042
1011 D1 = (X3 * T5 + Y3 * T6 + Z3 * T7) * W2
1012 D2 = (X5 * T5 + Y5 * T6 + Z5 * T7) * W2
1013 REM ----- COMPUTE PSI(.5,J,J+1)
1014 P1 = .5
1015 P2 = P3
1016 P3 = P3 + 1
1017 P4 = J2
1018 GOSUB 56
1019 U1 = T1
1020 U2 = T2
1021 REM ----- COMPUTE PSI(-.5,J,J+1)
1022 P1 = -P1
1023 GOSUB 56
1024 U1 = (T1 - U1) / S(J2)
1025 U2 = (T2 - U2) / S(J2)
1026 REM ----- COMPUTE PSI(.5,J-1,J)
1027 P1 = -P1
1028 P3 = P2
1029 P2 = P2 - 1
1030 P4 = J1
1031 GOSUB 56
1032 U3 = T1
1033 U4 = T2
1034 REM ----- COMPUTE PSI(-.5,J-1,J)
1035 P1 = -P1
1036 GOSUB 56
1037 REM ----- GRADIENT OF SCALAR POTENTIAL
1038 U5 = (U1 + (U3 - T1) / S(J1) + D1) * K + U5
1039 U6 = (U2 + (U4 - T2) / S(J1) + D2) * K + U6
1040 GOTO 1048
1041 REM ----- COMPONENTS OF VECTOR POTENTIAL A
1042 K!(1, J9) = K!(1, J9) + (X3 * CR(J) - X5 * CI(J)) * K
1043 K!(2, J9) = K!(2, J9) + (X5 * CR(J) + X3 * CI(J)) * K
1044 K!(3, J9) = K!(3, J9) + (Y3 * CR(J) - Y5 * CI(J)) * K
1045 K!(4, J9) = K!(4, J9) + (Y5 * CR(J) + Y3 * CI(J)) * K
1046 K!(5, J9) = K!(5, J9) + (Z3 * CR(J) - Z5 * CI(J)) * K
1047 K!(6, J9) = K!(6, J9) + (Z5 * CR(J) + Z3 * CI(J)) * K
1048 NEXT K
1049 IF N$ = "H" THEN 1052
1050 U7 = U5 * CR(J) - U6 * CI(J) + U7
1051 U8 = U6 * CR(J) + U5 * CI(J) + U8
1052 NEXT J
1053 IF N$ = "E" THEN 1075
1054 REM ----- DIFFERENCES OF VECTOR POTENTIAL A
1055 J8 = 1
1056 J9 = J9 + 1
1057 IF J9 = 2 THEN 964
1058 ON I GOTO 1059, 1064, 1069
1059 H(3) = K!(5, 1) - K!(5, 2)
1060 H(4) = K!(6, 1) - K!(6, 2)
1061 H(5) = K!(3, 2) - K!(3, 1)
1062 H(6) = K!(4, 2) - K!(4, 1)
1063 GOTO 1097
1064 H(1) = K!(5, 2) - K!(5, 1)
1065 H(2) = K!(6, 2) - K!(6, 1)
1066 H(5) = H(5) - K!(1, 2) + K!(1, 1)
1067 H(6) = H(6) - K!(2, 2) + K!(2, 1)
1068 GOTO 1097
1069 H(1) = H(1) - K!(3, 2) + K!(3, 1)
1070 H(2) = H(2) - K!(4, 2) + K!(4, 1)
1071 H(3) = H(3) + K!(1, 2) - K!(1, 1)
1072 H(4) = H(4) + K!(2, 2) - K!(2, 1)
1073 GOTO 1097
1074 REM ----- IMAGINARY PART OF ELECTRIC FIELD
1075 U7 = -M * U7 / S0
1076 REM ----- REAL PART OF ELECTRIC FIELD
1077 U8 = M * U8 / S0
1078 REM ----- MAGNITUDE AND PHASE CALCULATION
1079 S1 = 0
1080 IF (U7 = 0 AND U8 = 0) THEN 1082
1081 S1 = SQR(U7 * U7 + U8 * U8)
1082 S2 = 0
1083 IF U8 <> 0 THEN S2 = ATN(U7 / U8) / P0
1084 IF U8 > 0 THEN 1086
1085 S2 = S2 + SGN(U7) * 180
1086 IF I = 1 THEN PRINT #3, " X ",
1087 IF I = 2 THEN PRINT #3, " Y ",
1088 IF I = 3 THEN PRINT #3, " Z ",
1089 PRINT #3, TAB(15); F1 * U8; TAB(29); F1 * U7; TAB(43); F1 * S1; TAB(57); S2
1090 IF S$ = "Y" THEN PRINT #2, F1 * U8; ","; F1 * U7; ","; F1 * S1; ","; S2
1091 REM ----- CALCULATION FOR PEAK ELECTRIC FIELD
1092 S1 = S1 * S1
1093 S2 = S2 * P0
1094 A1 = A1 + S1 * COS(2 * S2)
1095 A3 = A3 + S1 * SIN(2 * S2)
1096 A4 = A4 + S1
1097 NEXT I
1098 IF N$ = "E" THEN 1121
1099 REM ----- MAGNETIC FIELD MAGNITUDE AND PHASE CALCULATION
1100 FOR I = 1 TO 5 STEP 2
1101 S1 = 0
1102 IF (H(I) = 0 AND H(I + 1) = 0) THEN 1104
1103 S1 = SQR(H(I) * H(I) + H(I + 1) * H(I + 1))
1104 S2 = 0
1105 IF H(I) <> 0 THEN S2 = ATN(H(I + 1) / H(I)) / P0
1106 IF H(I) > 0 THEN 1108
1107 S2 = S2 + SGN(H(I + 1)) * 180
1108 IF I = 1 THEN PRINT #3, " X ",
1109 IF I = 3 THEN PRINT #3, " Y ",
1110 IF I = 5 THEN PRINT #3, " Z ",
1111 PRINT #3, TAB(15); F1 * H(I); TAB(29); F1 * H(I + 1); TAB(43); F1 * S1; TAB(57); S2
1112 IF S$ = "Y" THEN PRINT #2, F1 * H(I); ","; F1 * H(I + 1); ","; F1 * S1; ","; S2
1113 REM ----- CALCULATION FOR PEAK MAGNETIC FIELD
1114 S1 = S1 * S1
1115 S2 = S2 * P0
1116 A1 = A1 + S1 * COS(2 * S2)
1117 A3 = A3 + S1 * SIN(2 * S2)
1118 A4 = A4 + S1
1119 NEXT I
1120 REM ----- PEAK FIELD CALCULATION
1121 PK = SQR(A4 / 2 + SQR(A1 * A1 + A3 * A3) / 2)
1122 PRINT #3, " MAXIMUM OR PEAK FIELD = "; F1 * PK; A$
1123 IF (S$ = "Y" AND N$ = "E") THEN PRINT #2, F1 * PK; ","; O2
1124 IF (S$ = "Y" AND N$ = "H") THEN PRINT #2, F1 * PK; ","; O2
1125 IF S$ = "Y" THEN PRINT #2, XX; ","; YY; ","; ZZ
1126 NEXT IX
1127 NEXT IY
1128 NEXT IZ
1129 CLOSE #2
1130 RETURN
1131 REM ********** FREQUENCY INPUT **********
1132 REM ----- SET FLAG
1133 PRINT
1134 INPUT "FREQUENCY (MHZ)"; F
1135 IF F = 0 THEN F = 299.8
1136 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, "FREQUENCY (MHZ):"; F
1137 W = 299.8 / F
1138 REM -----VIRTUAL DIPOLE LENGTH FOR NEAR FIELD CALCULATION
1139 S0 = .001 * W
1140 REM ----- 1 / (4 * PI * OMEGA * EPSILON)
1141 M = 4.77783352# * W
1142 REM ----- SET SMALL RADIUS MODIFICATION CONDITION
1143 SRM = .0001 * W
1144 PRINT #3, " WAVE LENGTH = "; W; " METERS"
1145 REM ----- 2 PI / WAVELENGTH
1146 W = 2 * P / W
1147 W2 = W * W / 2
1148 FLG = 0
1149 RETURN
1150 REM ********** GEOMETRY INPUT **********
1151 REM ----- WHEN GEOMETRY IS CHANGED, ENVIRONMENT MUST BE CHECKED
1152 GOSUB 1369
1153 PRINT
1154 IF INFILE THEN 1160
1155 INPUT "NO. OF WIRES"; NW
1156 IF NW = 0 THEN RETURN
1157 IF NW <= MW THEN 1160
1158 PRINT "NUMBER OF WIRES EXCEEDS DIMENSION..."
1159 GOTO 1155
1160 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, "NO. OF WIRES:"; NW
1161 REM ----- INITIALIZE NUMBER OF PULSES TO ZERO
1162 N = 0
1163 FOR I = 1 TO NW
1164 IF INFILE THEN GOSUB 1557: GOTO 1190
1165 PRINT
1166 PRINT "WIRE NO."; I
1167 INPUT " NO. OF SEGMENTS"; S1
1168 IF S1 = 0 THEN 1153
1169 A$ = " END ONE COORDINATES (X,Y,Z)"
1170 PRINT A$;
1171 INPUT X1, Y1, Z1
1172 IF G < 0 AND Z1 < 0 THEN PRINT "Z CANNOT BE NEGATIVE": GOTO 1170
1173 A$ = " END TWO COORDINATES (X,Y,Z)"
1174 PRINT A$;
1175 INPUT X2, Y2, Z2
1176 IF G < 0 AND Z2 < 0 THEN PRINT "Z CANNOT BE NEGATIVE": GOTO 1174
1177 IF X1 = X2 AND Y1 = Y2 AND Z1 = Z2 THEN PRINT "ZERO LENGTH WIRE.": GOTO 1166
1178 A$ = " RADIUS"
1179 PRINT " "; A$;
1180 INPUT A(I)
1181 IF A(I) <= 0! THEN 1179
1182 REM ----- DETERMINE CONNECTIONS
1183 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, "WIRE NO."; I
1184 GOSUB 1299
1185 PRINT "CHANGE WIRE NO. "; I; " (Y/N) ";
1186 INPUT A$
1187 IF A$ = "Y" THEN 1165
1188 IF A$ <> "N" THEN 1185
1189 REM ----- COMPUTE DIRECTION COSINES
1190 X3 = X2 - X1
1191 Y3 = Y2 - Y1
1192 Z3 = Z2 - Z1
1193 D = SQR(X3 * X3 + Y3 * Y3 + Z3 * Z3)
1194 CA(I) = X3 / D
1195 CB(I) = Y3 / D
1196 CG(I) = Z3 / D
1197 S(I) = D / S1
1198 REM ----- COMPUTE CONNECTIVITY DATA (PULSES N1 TO N)
1199 N1 = N + 1
1200 N(I, 1) = N1
1201 IF (S1 = 1 AND I1 = 0) THEN N(I, 1) = 0
1202 N = N1 + S1
1203 IF I1 = 0 THEN N = N - 1
1204 IF I2 = 0 THEN N = N - 1
1205 IF N > MP THEN PRINT "PULSE NUMBER EXCEEDS DIMENSION": CLOSE : GOTO 1155
1206 N(I, 2) = N
1207 IF (S1 = 1 AND I2 = 0) THEN N(I, 2) = 0
1208 IF N < N1 THEN 1247
1209 FOR J = N1 TO N
1210 C%(J, 1) = I
1211 C%(J, 2) = I
1212 W%(J) = I
1213 NEXT J
1214 C%(N1, 1) = I1
1215 C%(N, 2) = I2
1216 REM ----- COMPUTE COORDINATES OF BREAK POINTS
1217 I1 = N1 + 2 * (I - 1)
1218 I3 = I1
1219 X(I1) = X1
1220 Y(I1) = Y1
1221 Z(I1) = Z1
1222 IF C%(N1, 1) = 0 THEN 1230
1223 I2 = ABS(C%(N1, 1))
1224 F3 = SGN(C%(N1, 1)) * S(I2)
1225 X(I1) = X(I1) - F3 * CA(I2)
1226 Y(I1) = Y(I1) - F3 * CB(I2)
1227 IF C%(N1, 1) = -I THEN F3 = -F3
1228 Z(I1) = Z(I1) - F3 * CG(I2)
1229 I3 = I3 + 1
1230 I6 = N + 2 * I
1231 FOR I4 = I1 + 1 TO I6
1232 J = I4 - I3
1233 X(I4) = X1 + J * X3 / S1
1234 Y(I4) = Y1 + J * Y3 / S1
1235 Z(I4) = Z1 + J * Z3 / S1
1236 NEXT I4
1237 IF C%(N, 2) = 0 THEN 1245
1238 I2 = ABS(C%(N, 2))
1239 F3 = SGN(C%(N, 2)) * S(I2)
1240 I3 = I6 - 1
1241 X(I6) = X(I3) + F3 * CA(I2)
1242 Y(I6) = Y(I3) + F3 * CB(I2)
1243 IF I = -C%(N, 2) THEN F3 = -F3
1244 Z(I6) = Z(I3) + F3 * CG(I2)
1245 GOTO 1255
1246 REM ---- SINGLE SEGMEN 0 PULSE CASE
1247 I1 = N1 + 2 * (I - 1)
1248 X(I1) = X1
1249 Y(I1) = Y1
1250 Z(I1) = Z1
1251 I1 = I1 + 1
1252 X(I1) = X2
1253 Y(I1) = Y2
1254 Z(I1) = Z2
1255 NEXT I
1256 REM ********** GEOMETRY OUTPUT **********
1257 PRINT #3, " "
1258 PRINT #3, " **** ANTENNA GEOMETRY ****"
1259 IF N > 0 THEN 1264
1260 PRINT
1261 PRINT "NUMBER OF PULSES IS ZERO....RE-ENTER GEOMETRY"
1262 PRINT
1263 GOTO 1155
1264 K = 1
1265 J = 0
1266 FOR I = 1 TO N
1267 I1 = 2 * W%(I) - 1 + I
1268 IF K > NW THEN 1279
1269 IF K = J THEN 1279
1270 J = K
1271 PRINT #3, " "
1272 PRINT #3, "WIRE NO. "; K; " COORDINATES", , , "CONNECTION PULSE"
1273 PRINT #3, "X", "Y", "Z", "RADIUS", "END1 END2 NO."
1274 IF (N(K, 1) <> 0 OR N(K, 2) <> 0) THEN 1279
1275 PRINT #3, "-", "-", "-", " -", " - - 0"
1276 K = K + 1
1277 IF K > NW THEN 1286
1278 GOTO 1270
1279 PRINT #3, X(I1); TAB(15); Y(I1); TAB(29); Z(I1); TAB(43); A(W%(I)); TAB(57);
1280 PRINT #3, USING "### ### ##"; C%(I, 1), C%(I, 2), I
1281 IF (I = N(K, 2) OR N(K, 1) = N(K, 2) OR C%(I, 2) = 0) THEN K = K + 1
1282 IF C%(I, 1) = 0 THEN C%(I, 1) = W%(I)
1283 IF C%(I, 2) = 0 THEN C%(I, 2) = W%(I)
1284 IF (K = NW AND N(K, 1) = 0 AND N(K, 2) = 0) THEN 1270
1285 IF (I = N AND K < NW) THEN 1270
1286 NEXT I
1287 PRINT
1288 CLOSE 1: IF INFILE THEN INFILE = 0: IF O$ > "C" THEN 1293
1289 INPUT " CHANGE GEOMETRY (Y/N) "; A$
1290 IF A$ = "Y" THEN 1153
1291 IF A$ <> "N" THEN 1289
1292 REM ----- EXCITATION INPUT
1293 GOSUB 1430
1294 REM ----- LOADS/NETWORKS INPUT
1295 GOSUB 1455
1296 FLG = 0
1297 RETURN
1298 REM ********** CONNECTIONS **********
1299 E(I) = X1
1300 L(I) = Y1
1301 M(I) = Z1
1302 E(I + NW) = X2
1303 L(I + NW) = Y2
1304 M(I + NW) = Z2
1305 G% = 0
1306 I1 = 0
1307 I2 = 0
1308 J1(I) = 0
1309 J2(I, 1) = -I
1310 J2(I, 2) = -I
1311 IF G = 1 THEN 1323
1312 REM ----- CHECK FOR GROUND CONNECTION
1313 IF Z1 = 0 THEN 1315
1314 GOTO 1318
1315 I1 = -I
1316 J1(I) = -1
1317 GOTO 1340
1318 IF Z2 = 0 THEN 1320
1319 GOTO 1323
1320 I2 = -I
1321 J1(I) = 1
1322 G% = 1
1323 IF I = 1 THEN 1358
1324 FOR J = 1 TO I - 1
1325 REM ----- CHECK FOR END1 TO END1
1326 IF (X1 = E(J) AND Y1 = L(J) AND Z1 = M(J)) THEN 1328
1327 GOTO 1333
1328 I1 = -J
1329 J2(I, 1) = J
1330 IF J2(J, 1) = -J THEN J2(J, 1) = J
1331 GOTO 1340
1332 REM ----- CHECK FOR END1 TO END2
1333 IF (X1 = E(J + NW) AND Y1 = L(J + NW) AND Z1 = M(J + NW)) THEN 1335
1334 GOTO 1339
1335 I1 = J
1336 J2(I, 1) = J
1337 IF J2(J, 2) = -J THEN J2(J, 2) = J
1338 GOTO 1340
1339 NEXT J
1340 IF G% = 1 THEN 1358
1341 IF I = 1 THEN 1358
1342 FOR J = 1 TO I - 1
1343 REM ----- CHECK END2 TO END2
1344 IF (X2 = E(J + NW) AND Y2 = L(J + NW) AND Z2 = M(J + NW)) THEN 1346
1345 GOTO 1351
1346 I2 = -J
1347 J2(I, 2) = J
1348 IF J2(J, 2) = -J THEN J2(J, 2) = J
1349 GOTO 1358
1350 REM ----- CHECK FOR END2 TO END1
1351 IF (X2 = E(J) AND Y2 = L(J) AND Z2 = M(J)) THEN 1353
1352 GOTO 1357
1353 I2 = J
1354 J2(I, 2) = J
1355 IF J2(J, 1) = -J THEN J2(J, 1) = J
1356 GOTO 1358
1357 NEXT J
1358 PRINT #3, " COORDINATES", " ", " ", "END NO. OF"
1359 PRINT #3, " X", " Y", " Z", "RADIUS CONNECTION SEGMENTS"
1360 PRINT #3, X1; TAB(15); Y1; TAB(29); Z1; TAB(57); I1
1361 PRINT #3, X2; TAB(15); Y2; TAB(29); Z2; TAB(43); A(I); TAB(57); I2; TAB(71); S1
1362 RETURN
1363 REM ********** ENVIROMENT INPUT **********
1364 PRINT
1365 PRINT " **** WARNING ****"
1366 PRINT "REDO GEOMETRY TO ENSURE PROPER GROUND CONNECTION/DISCONNECTION"
1367 PRINT
1368 REM ----- INITIALIZE NUMBER OF RADIAL WIRES TO ZERO
1369 NR = 0
1370 REM ----- SET ENVIRONMENT
1371 PRINT #3, " "
1372 A$ = "ENVIRONMENT (+1 FOR FREE SPACE, -1 FOR GROUND PLANE)"
1373 PRINT A$;
1374 INPUT G
1375 IF O$ > "C" THEN PRINT #3, A$; ": "; G
1376 IF G = 1 THEN 1428
1377 IF G <> -1 THEN 1373
1378 REM ----- NUMBER OF MEDIA
1379 A$ = " NUMBER OF MEDIA (0 FOR PERFECTLY CONDUCTING GROUND)"
1380 PRINT A$;
1381 INPUT NM
1382 IF NM <= MM THEN 1385
1383 PRINT "NUMBER OF MEDIA EXCEEDS DIMENSION..."
1384 GOTO 1380
1385 IF O$ > "C" THEN PRINT #3, A$; ": "; NM
1386 REM ----- INITIALIZE BOUNDARY TYPE
1387 TB = 1
1388 IF NM = 0 THEN 1428
1389 IF NM = 1 THEN 1396
1390 REM ----- TYPE OF BOUNDARY
1391 A$ = " TYPE OF BOUNDARY (1-LINEAR, 2-CIRCULAR)"
1392 PRINT " "; A$;
1393 INPUT TB
1394 IF O$ > "C" THEN PRINT #3, A$; ": "; TB
1395 REM ----- BOUNDARY CONDITIONS
1396 FOR I = 1 TO NM
1397 PRINT "MEDIA"; I
1398 A$ = " RELATIVE DIELECTRIC CONSTANT, CONDUCTIVITY"
1399 PRINT " "; A$;
1400 INPUT T(I), V(I)
1401 IF O$ > "C" THEN PRINT #3, A$; ": "; T(I); ","; V(I)
1402 IF I > 1 THEN 1414
1403 IF TB = 1 THEN 1414
1404 A$ = " NUMBER OF RADIAL WIRES IN GROUND SCREEN"
1405 PRINT " "; A$;
1406 INPUT NR
1407 IF O$ > "C" THEN PRINT #3, A$; ": "; NR
1408 IF NR = 0 THEN 1414
1409 A$ = " RADIUS OF RADIAL WIRES"
1410 PRINT " "; A$;
1411 INPUT RR
1412 IF O$ > "C" THEN PRINT #3, A$; ": "; RR
1413 REM ----- INITIALIZE COORDINATE OF MEDIA INTERFACE
1414 U(I) = 1000000!
1415 REM ----- INITIALIZE HEIGHT OF MEDIA
1416 H(I) = 0
1417 IF I = NM THEN 1422
1418 A$ = " X OR R COORDINATE OF NEXT MEDIA INTERFACE"
1419 PRINT " "; A$;
1420 INPUT U(I)
1421 IF O$ > "C" THEN PRINT #3, A$; ": "; U(I)
1422 IF I = 1 THEN 1427
1423 A$ = " HEIGHT OF MEDIA"
1424 PRINT " "; A$;
1425 INPUT H(I)
1426 IF O$ > "C" THEN PRINT #3, A$; ": "; H(I)
1427 NEXT I
1428 RETURN
1429 REM ********** EXCITATION INPUT **********
1430 PRINT
1431 A$ = "NO. OF SOURCES "
1432 PRINT A$;
1433 INPUT NS
1434 IF NS < 1 THEN NS = 1
1435 IF NS <= MP THEN 1438
1436 PRINT "NO. OF SOURCES EXCEEDS DIMENSION ..."
1437 GOTO 1432
1438 IF O$ > "C" THEN PRINT #3, " ": PRINT #3, A$; ": "; NS
1439 FOR I = 1 TO NS
1440 PRINT
1441 PRINT "SOURCE NO. "; I; ":"
1442 A$ = "PULSE NO., VOLTAGE MAGNITUDE, PHASE (DEGREES)"
1443 PRINT A$;
1444 INPUT E(I), VM, VP
1445 IF E(I) <= N THEN 1448
1446 PRINT "PULSE NUMBER EXCEEDS NUMBER OF PULSES..."
1447 GOTO 1443
1448 IF O$ > "C" THEN PRINT #3, A$; ": "; E(I); ","; VM; ","; VP
1449 L(I) = VM * COS(VP * P0)
1450 M(I) = VM * SIN(VP * P0)
1451 NEXT I
1452 IF FLG = 2 THEN FLG = 1
1453 RETURN
1454 REM ********** LOADS INPUT **********
1455 PRINT
1456 INPUT "NUMBER OF LOADS "; NL
1457 IF NL <= ML THEN 1460
1458 PRINT "NUMBER OF LOADS EXCEEDS DIMENSION..."
1459 GOTO 1456
1460 IF O$ > "C" THEN PRINT #3, "NUMBER OF LOADS"; NL
1461 IF NL < 1 THEN 1492
1462 INPUT "S-DOMAIN (S=jw) IMPEDANCE LOAD (Y/N)"; L$
1463 IF L$ <> "Y" AND L$ <> "N" THEN 1462
1464 A$ = "PULSE NO.,RESISTANCE,REACTANCE"
1465 IF L$ = "Y" THEN A$ = "PULSE NO., ORDER OF S-DOMAIN FUNCTION"
1466 FOR I = 1 TO NL
1467 PRINT
1468 PRINT "LOAD NO. "; I; ":"
1469 IF L$ = "Y" THEN 1476
1470 PRINT A$;
1471 INPUT LP(I), LA(1, I, 1), LA(2, I, 1)
1472 IF LP(I) > N THEN PRINT "PULSE NUMBER EXCEEDS NUMBER OF PULSES...": GOTO 1470
1473 IF O$ > "C" THEN PRINT #3, A$; ": "; LP(I); ","; LA(1, I, 1); ","; LA(2, I, 1)
1474 GOTO 1491
1475 REM ----- S-DOMAIN LOADS
1476 PRINT A$;
1477 INPUT LP(I), LS(I)
1478 IF LP(I) > N THEN PRINT "PULSE NUMBER EXCEEDS NUMBER OF PULSES...": GOTO 1476
1479 IF LS(I) > MA THEN PRINT "MAXIMUM DIMENSION IS 10": GOTO 1477
1480 IF O$ > "C" THEN PRINT #3, A$; ": "; LP(I); ","; LS(I)
1481 FOR J = 0 TO LS(I)
1482 A1$ = "NUMERATOR, DENOMINATOR COEFFICIENTS OF S^"
1483 PRINT A1$; J;
1484 INPUT LA(1, I, J), LA(2, I, J)
1485 IF O$ > "C" THEN PRINT #3, A1$; J; ":"; LA(1, I, J); ","; LA(2, I, J)
1486 NEXT J
1487 IF LS(I) > 0 THEN 1491
1488 LS(I) = 1
1489 LA(1, I, 1) = 0
1490 LA(2, I, 1) = 0
1491 NEXT I
1492 FLG = 0
1493 RETURN
1494 REM ********** MAIN PROGRAM **********
1495 REM ----- DATA INITIALIZATION
1496 REM ----- PI
1497 P = 4 * ATN(1)
1498 REM ----- CHANGES DEGREES TO RADIANS
1499 P0 = P / 180
1500 B$ = "********************"
1501 REM ----- INTRINSIC IMPEDANCE OF FREE SPACE DIVIDED BY 2 PI
1502 G0 = 29.979221#
1503 REM ---------- Q-VECTOR FOR GAUSSIAN QUADRATURE
1504 READ Q(1), Q(2), Q(3), Q(4), Q(5), Q(6), Q(7), Q(8), Q(9), Q(10), Q(11), Q(12)
1505 READ Q(13), Q(14)
1506 DATA .288675135,.5,.430568156,.173927423,.169990522,.326072577
1507 DATA .480144928,.050614268,.398333239,.111190517
1508 DATA .262766205,.156853323,.091717321,.181341892
1509 REM ---------- E-VECTOR FOR COEFFICIENTS OF ELLIPTIC INTEGRAL
1510 READ C0, C1, C2, C3, C4, C5, C6, C7, C8, C9
1511 DATA 1.38629436112,.09666344259,.03590092383,.03742563713,.01451196212
1512 DATA .5,.12498593397,.06880248576,.0332835346,.00441787012
1513 REM ----- IDENTIFY OUTPUT DEVICE
1514 GOSUB 1580
1515 PRINT #3, TAB(20); B$; B$
1516 PRINT #3, TAB(22); "MINI-NUMERICAL ELECTROMAGNETICS CODE"
1517 PRINT #3, TAB(34); "MININEC (3)"
1518 PRINT #3, TAB(24); DATE$; TAB(48); TIME$
1519 PRINT #3, TAB(20); B$; B$
1520 REM ----- FREQUENCY INPUT
1521 GOSUB 1133
1522 REM ----- ENVIRONMENT INPUT
1523 GOSUB 1369
1524 REM ----- CHECK FOR NEC-TYPE GEOMETRY INPUT
1525 GOSUB 1550
1526 REM ----- GEOMETRY INPUT
1527 GOSUB 1153
1528 REM ----- MENU
1529 PRINT
1530 PRINT B$; " MININEC(3) MENU "; B$
1531 PRINT " G - CHANGE GEOMETRY C - COMPUTE/DISPLAY CURRENTS"
1532 PRINT " E - CHANGE ENVIRONMENT P - COMPUTE FAR-FIELD PATTERNS"
1533 PRINT " X - CHANGE EXCITATION N - COMPUTE NEAR-FIELDS"
1534 PRINT " L - CHANGE LOADS"
1535 PRINT " F - CHANGE FREQUENCY Q - QUIT"
1536 PRINT B$; " {math co-processor} "; B$
1537 INPUT " COMMAND "; C$
1538 IF C$ = "F" THEN GOSUB 1133
1539 IF C$ = "P" THEN GOSUB 621
1540 IF C$ = "X" THEN GOSUB 1430
1541 IF C$ = "E" THEN GOSUB 1364
1542 IF C$ = "G" THEN GOSUB 1152
1543 IF C$ = "C" THEN GOSUB 497
1544 IF C$ = "L" THEN GOSUB 1455
1545 IF C$ = "N" THEN GOSUB 875
1546 IF C$ <> "Q" THEN 1529
1547 IF O$ = "P" THEN PRINT #3, CHR$(12) ELSE IF O$ = "C" THEN PRINT #3, " "
1548 CLOSE
1549 SYSTEM
1550 REM ********** NEC-TYPE GEOMETRY INPUT **********
1551 OPEN "MININEC.INP" FOR RANDOM AS #1 LEN = 30
1552 FIELD #1, 2 AS S$, 4 AS X1$, 4 AS Y1$, 4 AS Z1$, 4 AS X2$, 4 AS Y2$, 4 AS Z2$, 4 AS R$
1553 GET 1
1554 NW = CVI(S$)
1555 IF NW THEN INFILE = 1
1556 RETURN
1557 REM ---------- GET GEOMETRY DATA FROM MININEC.INP
1558 GET 1
1559 S1 = CVI(S$)
1560 X1 = CVS(X1$)
1561 Y1 = CVS(Y1$)
1562 Z1 = CVS(Z1$)
1563 X2 = CVS(X2$)
1564 Y2 = CVS(Y2$)
1565 Z2 = CVS(Z2$)
1566 A(I) = CVS(R$)
1567 IF G < 0 THEN IF Z1 < 0 OR Z2 < 0 THEN GOSUB 1572
1568 PRINT #3, " ": PRINT #3, "WIRE NO."; I
1569 IF X1 = X2 AND Y1 = Y2 AND Z1 = Z2 THEN PRINT "WIRE LENGTH IS ZERO.": GOTO 1547
1570 GOSUB 1299
1571 RETURN
1572 IF IZNEG THEN 1576
1573 PRINT "NEGATIVE Z VALUE ENCOUNTERED FOR GROUND PLANE."
1574 INPUT "ABORT OR CONVERT NEGATIVE Z VALUE TO ZERO (A/C)? "; A$
1575 IF A$ = "A" THEN 1547 ELSE IF A$ = "C" THEN IZNEG = 1 ELSE 1574
1576 IF Z1 < 0 THEN Z1 = -Z1
1577 IF Z2 < 0 THEN Z2 = -Z2
1578 RETURN
1579 REM ********** IDENTIFY OUTPUT DEVICE **********
1580 INPUT "OUTPUT TO CONSOLE, PRINTER, OR DISK (C/P/D)"; O$
1581 IF O$ = "C" THEN F$ = "SCRN:": GOTO 1586
1582 IF O$ = "P" THEN F$ = "LPT1:": GOTO 1586
1583 IF O$ <> "D" THEN 1580
1584 INPUT "FILENAME (NAME.OUT)"; F$
1585 IF LEFT$(RIGHT$(F$, 4), 1) = "." THEN 1586 ELSE F$ = F$ + ".OUT"
1586 OPEN F$ FOR OUTPUT AS #3
1587 CLS
1588 RETURN
1589 REM ********** CALCULATE ELAPSED TIME **********
1590 JH = VAL(MID$(T$, 1, 2)) - VAL(MID$(OT$, 1, 2))
1591 JM = VAL(MID$(T$, 4, 2)) - VAL(MID$(OT$, 4, 2))
1592 JS = VAL(MID$(T$, 7, 2)) - VAL(MID$(OT$, 7, 2))
1593 IF JS < 0 THEN JS = JS + 60: JM = JM - 1
1594 IF JM < 0 THEN JM = JM + 60: JH = JH - 1
1595 IF JH < 0 THEN JH = JH + 24
1596 T$ = ":" + MID$(STR$(JS + 100), 3)
1597 IF JH THEN T$ = MID$(STR$(JH), 2) + ":" + MID$(STR$(JM + 100), 3) + T$ ELSE T$ = MID$(STR$(JM), 2) + T$
1598 RETURN
1599 REM ********** CALCULATE APPROXIMATE TIME REMAINING **********
1600 IPCT = 100 * PCT
1601 T$ = TIME$
1602 JH = VAL(MID$(T$, 1, 2)) - VAL(MID$(OT$, 1, 2))
1603 IF JH < 0 THEN JH = JH + 24
1604 JM = VAL(MID$(T$, 4, 2)) - VAL(MID$(OT$, 4, 2))
1605 JS = VAL(MID$(T$, 7, 2)) - VAL(MID$(OT$, 7, 2))
1606 JS = JS + 60 * (JM + 60 * JH)
1607 JS = JS * (1 / PCT - 1)
1608 JM = INT(JS / 60)
1609 JS = JS MOD 60
1610 JH = INT(JM / 60)
1611 JM = JM MOD 60
1612 T$ = ":" + MID$(STR$(JS + 100), 3)
1613 IF JH THEN T$ = MID$(STR$(JH), 2) + ":" + MID$(STR$(JM + 100), 3) + T$ ELSE T$ = MID$(STR$(JM), 2) + T$
1614 LOCATE CSRLIN, 1
1615 PRINT Q$; IPCT; "% COMPLETE - APPROX TIME REMAINING "; T$; " ";
1616 RETURN
1617 END
REM $STATIC
DEFSNG I-K, N