[Prev][Next][Index][Thread]
Re: Magnifier system
[The following text is in the "ISO-8859-1" character set]
[Your display is set for the "US-ASCII" character set]
[Some characters may be displayed incorrectly]
Hi Malcolm,
> Subject: Re: Magnifier system
> > [ for anyone interested, I have a QuickBAsic program which will solve
the
> > coupled differential equations for up to 100 L's ( with user specified
> > mutual coupling) C's and R's in a tesla coil-like structure - the
voltage
> > or current at each node is plotted in time so that one can see how
> > discretized waves propegate on such a structure].
>
> Yes please! Anything that adds to my understanding of this
> fascinating structure is welcomed. Would you be able to post it as
> plain ASCII text please? I have endless difficulties with various
> coding schemes because I am working on a LAN system with all sorts of
> restrictions imposed.
Here it comes, let me know if it survived intact! I have done some
commenting
whihc will hopefully fill in the flavor of what I was trying to do so that
you can make any changes you would like. I'd be happy to answer any
questions about how the code works:
A single quote is the comment
++++++++++++++++clip here ++++++++++++++++++++++
' A discrete transmission line simulation program
' by Ed Harris
' a transmision line is simuilated by solving the circuit
' equations for n nodes of a discrete lumped LCR circuit with
' arbitrary mutual inductances
' its fun to run just to see the output -- forget about the physics!
' The basic circuit layout (may not come out lined up properly):
' Vin --R1---L1--.--R2---L2--.--- etc -------.--Rn---Ln--. Vn
' | | | |
' ____ C1 ____ C2 _____ Cn-1 _____ Cn
' ____ ____ _____ _____
' | | | |
' GND GND GND GND
' M12=M21 (not shown) is the mutual inductance between inductors
' L1 and L2
' the following line asks for the number of nodes you wish to
' use. I suggest starting with 20 or so just so that the
' program runs more swiftly. The only limitaion really
' is that QuickBasic has some memory limitaions -- the code does
' not have any inherent size limitaions however
INPUT "Numnber of nodes to run( 20-50 is good to start with)"; n
' set up for QuickBasic graphics
SCREEN 9
CLS
DIM m(n, n), y(n, n), vv(n), indx(n), b(n)
DIM oc(n), r(n), i(n), q(n), di(n), dq(n), v(n), f(n)
'I will fill the m matrix with pseudo inductances and mutual inductances
' in this case ---
' m( i,j) is the mutual inductance between nodes i and j
' m(i,i) is the self inductance of the node
' I just made up a function here for fun---
' you can set the m values to anything you like but remember that
' m(i,j)=m(j,i) by symmetery right?
' anyway for realistic computations I'd suggest using Grover's
' formulas to get the m matrix
' for untis consideration see below
FOR i% = 1 TO n
FOR j% = 1 TO n
m(i%, j%) = 1 / ((i% - j%) ^ 2 + 1)
NEXT j%
NEXT i%
' fill r with same r
' this is the nodal series resistance of the nodal inductance
' see note about units below
FOR i% = 1 TO n
r(i%) = .01
NEXT i%
'fill 1/c matrix with c's
' oc(node)=1/c(node)
' Ive got unitless variables in here now, but
' any units can be used as long as one is
' consistent: ie VOLTS/FARADS/AMPERES/HENRIES/COULOMBS etc...
FOR i% = 1 TO n
oc(i%) = 1
NEXT i%
' this is commented out - i used it to check the form of
' the mutual inductance matrix
'FOR i% = 1 TO n
'FOR j% = 1 TO n
'PRINT USING "#.## "; m(i%, j%);
'NEXT j%
'PRINT
'NEXT i%
' the following code up to the biterate lable is all just for
' inverting the mutual inductance matrix -- it's just the
' Numerical Recipes routine adapted for QUick Basic
FOR i% = 1 TO n
FOR j% = 1 TO n
y(i%, j%) = 0
NEXT j%
y(i%, i%) = 1
NEXT i%
d = 1
FOR i% = 1 TO n
aamax = 0
FOR j% = 1 TO n
IF (ABS(m(i%, j%)) > aamax) THEN
aamax = ABS(m(i%, j%))
END IF
NEXT j%
IF aamax = 0 THEN
PRINT "error-singular matrix"
END IF
vv(i%) = 1 / aamax
NEXT i%
FOR j% = 1 TO n
IF j% > 1 THEN
FOR i% = 1 TO j% - 1
sum = m(i%, j%)
IF i% > 1 THEN
FOR k% = 1 TO i% - 1
sum = sum - m(i%, k%) * m(k%, j%)
NEXT k%
m(i%, j%) = sum
END IF
NEXT i%
END IF
aamax = 0
FOR i% = j% TO n
sum = m(i%, j%)
IF j% > 1 THEN
FOR k% = 1 TO j% - 1
sum = sum - m(i%, k%) * m(k%, j%)
NEXT k%
m(i%, j%) = sum
END IF
dum = vv(i%) * ABS(sum)
IF dum > aamax THEN
imax = i%
aamax = dum
END IF
NEXT i%
IF j% <> imax THEN
FOR k% = 1 TO n
dum = m(imax, k%)
m(imax, k%) = m(j%, k%)
m(j%, k%) = dum
NEXT k%
d = -d
vv(imax) = vv(j%)
END IF
indx(j%) = imax
IF j% <> n THEN
IF m(j%, j%) = 0 THEN
m(j%, j%) = tiny
END IF
dum = 1 / m(j%, j%)
FOR i% = j% + 1 TO n
m(i%, j%) = m(i%, j%) * dum
NEXT i%
END IF
NEXT j%
IF m(n, n) = 0 THEN
m(n, n) = tiny
END IF
FOR p = 1 TO n
FOR q = 1 TO n
b(q) = y(q, p)
NEXT q
GOSUB lubksb
FOR q = 1 TO n
y(q, p) = b(q)
NEXT q
NEXT p
GOTO biterate
lubksb:
ii = 0
FOR i% = 1 TO n
ll = indx(i%)
sum = b(ll)
IF ii <> 0 THEN
FOR j% = ii TO i% - 1
sum = sum - m(i%, j%) * b(j%)
NEXT j%
ELSE
IF sum <> 0 THEN
ii = i%
END IF
END IF
b(i%) = sum
NEXT i%
FOR i% = n TO 1 STEP -1
sum = b(i%)
IF i% < n THEN
FOR j% = i% + 1 TO n
sum = sum - m(i%, j%) * b(j%)
NEXT j%
END IF
b(i%) = sum / m(i%, i%)
NEXT i%
RETURN
' the matrix inversion routine terminates here - just before
' the diff eq's are iterated
biterate:
'FOR i% = 1 TO n
'FOR j% = 1 TO n
'PRINT USING "##.## "; y(i%, j%);
'NEXT j%
'PRINT
'NEXT i%
CLS
' set the intial period equal to the number of nodes
p = n
' the next line is the main iteration loop lable
gloop:
t = t + dt ' increment the time by dt
'the next line is the pulsed voltage applied to the first node
' this can be cahnged to any function one would want -- another
' example is commented out and show below: drive with sine wave
Vin = 20 * INT(.5 * SIN(t / p) + .501)
'Vin = SIN(t / p)
a$ = INKEY$
'some interactive command keys:
' pressing u on the keyboard increments the period of the driving
' functions defined above
' d lowers the period
' q is used to quit in compiled verions of the program
' c prints then clears the vm (voltage maximum) variable which I used to
' track the max voltage on the last node
IF a$ = "u" THEN p = p + .1
IF a$ = "d" THEN p = p - .1
IF a$ = "q" THEN STOP
IF v(n) > vm THEN vm = v(n)
IF a$ = "c" THEN vm = 0
LOCATE 1, 1
PRINT p
PRINT vm
' time increment
dt = .5
' the following code does a fisrt order increment in the
' voltages and currents and charges
' q() is the charge on a gien nodal capacitor
' oc() is 1/capacitance of a given node
' i() is the loop current the mesh associated with a pair of nodes
' r() is the resistnace in series with the inductance at a node
' v() is the voltage across the capacitor in a node
' y( , ) matrix is the inverse mutual inductance matrix
f(1) = dt * (Vin - i(1) * r(1) - q(1) * oc(1))
FOR k% = 2 TO n
f(k%) = dt * (q(k% - 1) * oc(k%) - i(k%) * r(k%) - q(k%) * oc(k%))
NEXT k%
FOR k% = 1 TO n
sum = 0
FOR j% = 1 TO n
sum = sum + y(k%, j%) * f(j%)
NEXT j%
di(k%) = sum
i(k%) = i(k%) + di(k%)
NEXT k%
dq(n) = dt * i(n)
q(n) = q(n) + dq(n)
v(n) = q(n) * oc(n)
FOR k% = 1 TO n - 1
dq(k%) = dt * (i(k%) - i(k% + 1))
q(k%) = q(k%) + dq(k%)
v(k%) = q(k%) * oc(k%)
NEXT k%
' the following lines of code do the spatial waveform plotting
' node postion are incremented vertically with the voltage plotted
' ou horizontally from the node --- any quantity could abe added
'in for plotting
scale = 300 / n
LINE (0, 10)-(600, 10), 0
LINE (200, 10)-(200 + Vin * 10, 10), 4
FOR i% = 1 TO n
LINE (0, i% * scale + 10)-(600, i% * scale + 10), 0
LINE (200, i% * scale + 10)-(200 + v(i%) * 10, i% * scale + 10), 2
NEXT i%
' bad programming form eh?
GOTO gloop