DECLARE SUB invertH (nn%, y!()) DECLARE SUB dataerror () DECLARE SUB printfour (i%, ib%(), b!(), a$, g$) DECLARE SUB readnonlin (i%, b!(), nnl%, cn!()) DECLARE SUB examine (g$) DECLARE SUB modifyL (nbran%, nls%, ib%(), b!(), y!(), tag!()) DECLARE SUB nonlinZ (branch%, ib%(), b!(), e3!(), ci3!(), cn!(), t!, t1!, K9%) DECLARE SUB buildy (i%, nn%, nbran%, nsorce%, nidep%, ib%(), b!(), t1!, y!(), e1!()) DECLARE SUB waveshape (t, t3!(), ix3%, v4!, v5!, v6!, v7!, vmag!) DECLARE SUB addtoY (i%, ib%(), y1!, j1%, k1%, j2%, k2%, y!()) ' This program is TRANCALC , a program to do transient analysis ' of simple RLC circuits, including non-distorting transmission ' lines. It was developed by FA Fisher ' Lightning Technologies Inc. ' 10 Downing Parkway ' Pittsfield, Ma 01201 ' (413) 499-2135 ' ' ' ************************************ version$ = " this version dated 14 May 1995" 'update after changes ' *********************************** ' PRINT PRINT " This is 'TRANCALC' "; version$ ' PRINT : PRINT PRINT " This program comes to you from F. A. Fisher" PRINT " Lightning Technologies Inc." PRINT " 10 Downing Parkway " PRINT " Pittsfield, MA 01201" PRINT : PRINT ' ' ------------dimension statements---------" ' 100 '$DYNAMIC CLEAR DEFINT I-K, M-N ' DEFDBL B-C, E, G, L, T, W-Y 'use single precision for large ckts ' KEY(1) ON 'allows f(1) to interrupt execution ' maxnode = 100 'max number of nodes maxbran = 150 'max number of branches maxind = 40 'max number of inductors maxntl = 12 'max number of transmission lines nlinestep = 400 'max length of t-lines in delta t nmax = 1023 'max number of output points noutpmax = 10 'max number of output locations DIM ib(maxbran, 5) 'integer portion of branch table DIM b(maxbran, 6) 'floating point portion of branch table DIM y(maxnode, maxnode) 'branch admittance array DIM tag(maxind, maxind) 'tags inductor numbers for inversion DIM e1(16, 8) 'source table DIM cn(40, 16) 'non-linear data table DIM ci3(maxbran, 6) 'branch current table used for calculations DIM e3(maxbran, 8) 'branch I and v resulting from calculations DIM ci4(maxnode, 1) 'nodal current vector DIM e2(maxnode, 3) 'nodal voltage vector DIM l9(maxntl, 10) 'properties of transmission lines DIM q(maxntl) 'index to transmission line DIM q1(maxntl) 'length of transmission lines in delta t DIM g(6, 9) 'conventional switch table DIM g1(6, 7) 'semiconductor switch table DIM t3(20, 2) 'data for tabular source DIM descrip$(3) 'descriptive data for case DIM l1(nlinestep, maxntl) 'forward voltages for transmission lines DIM l2(nlinestep, maxntl) 'reverse voltages for transmission lines DIM l3(nlinestep, maxntl) 'forward currents for transmission lines DIM l4(nlinestep, maxntl) 'reverse currents for transmission lines DIM o$(noutpmax) 'names for type of output DIM outtype%(noutpmax) 'numbers corresponding to the names of type DIM outplace%(noutpmax) 'number of branch or node in output point ' ' pi = 3.14159 Pi2 = pi / 2 scrn = 12 'VGA monitor assumed colset = 14 vran = 360 godown = 13 WIDTH 80 CLS SCREEN 0 'initialize monitor COLOR colset ' FOR i = 2 TO 16 STEP 2 'set cn(2,i)=1 so that non -linear cn(2, i) = 1 'branches are all checked on start NEXT i ' 'jout = 0 'set first line of output to zero ' K9 = 1 'set flag for new topology - must formulate y matrix ' '-----------initialize counters for various types of elements---- ' ntl = 0 'counts transmission lines nsw = 0 'counts switches nsc = 0 'counts semiconductor switches nnl = 0 'counts non-linear branches nls = 0 'counts self inductors nlm = 0 'counts mutual inductors ' CLS ' 112 ON ERROR GOTO 114 ' PRINT PRINT " What file and drive (eg. A:filename.Dat ) (-99 to abort)" INPUT infile$ ' IF infile$ = "-99" THEN ' read name of data file END ELSE OPEN infile$ FOR INPUT AS #1 ' open file for input data GOTO 116 END IF ' 114 BEEP: BEEP: BEEP PRINT "Bad file name or disk descriptor----try again (-99 to quit)" RESUME 112 ' 116 outfile$ = LEFT$(infile$, INSTR(infile$, ".") - 1) + ".out" 117 PRINT "Output file name will be "; outfile$; " OK (y/n) "; ' 118 a$ = INKEY$: IF a$ = "" THEN 118 IF a$ = "n" OR a$ = "N" THEN INPUT "new filename "; outfile$ GOTO 117 END IF ' OPEN "tranpass" FOR OUTPUT AS #2 ' open file to pass filenames datafile$ = infile$ WRITE #2, datafile$, outfile$ CLOSE #2 OPEN outfile$ FOR OUTPUT AS #2 ' 120 ON ERROR GOTO 0 ' CLS PRINT "This run is made from the data in file "; infile$: PRINT ' LINE INPUT #1, descrip$(1) PRINT descrip$(1) a$ = LEFT$(descrip$(1), 1) IF a$ = "" THEN PRINT "You may have a blank line in your data file." PRINT "IF so, it will cause a run time error>" CALL dataerror END IF ' a$ = LEFT$(descrip$(1), 2) a1$ = LEFT$(a$, 1) IF a1$ = CHR$(34) THEN a$ = RIGHT$(a$, 1) ELSE a$ = a1$ IF a$ = "T" OR a$ = "t" THEN GOTO 130 ELSE PRINT "This is not marked as a transient data file" CALL dataerror END IF ' 130 LINE INPUT #1, descrip$(2) 'read case description lines LINE INPUT #1, descrip$(3) PRINT descrip$(2) PRINT descrip$(3) PRINT INPUT #1, nnode, nbran, t1, t2, t3, t4, increm 'nbran = number of branches PRINT PRINT nnode; " nodes and "; nbran; " branches" ' IF nnode > maxnode OR nbran > maxbran THEN PRINT "Only "; maxnode; " nodes and "; maxbran; " branches permitted" CALL dataerror END IF ' PRINT ' f$ = "& +#.###^^^^ & +#.###^^^^ " PRINT USING f$; "time step is "; t1; " and end time is "; t2 PRINT USING f$; " after "; t3; " the time step is "; t4 PRINT " 1 point out of every "; increm; " is saved" ' ' t1 = delta t, t2 = final time, t3=time to change delta t ' t4 = new delta t, t5 = increment for storing output ' tmax = t2 WRITE #2, tmax PRINT ' g$ = "n" 'suppresses stops to examine branch data ' ' --------read data and store in branch table--- ' PRINT "branch from to type value1 value2 " PRINT FOR i = 1 TO nbran 'end of loop at 150 'integer data stored in ib, floating point in b ' INPUT #1, ib(i, 1), ib(i, 2), ib(i, 3), ib(i, 4), b(i, 1) 'read first 5 IF ib(i, 1) > nbran THEN PRINT "There was an error in the branch number" PRINT "What was read was branch number = "; ib(i, 1) PRINT " from node = "; ib(i, 2) PRINT " to node = "; ib(i, 3) PRINT " branch type = "; ib(i, 4) PRINT " value = "; ib(i, 5) CALL dataerror ELSE IF ib(i, 4) < 1 OR ib(i, 4) > 16 THEN PRINT "Error in type when reading branch "; i; " Elements were:" FOR jj = 1 TO 5 PRINT "b("; i; ","; jj; ") = "; b(i, jj) NEXT jj PRINT "The error may actually be in the preceeding line of data." CALL dataerror END IF END IF ' IF ib(i, 4) = 8 THEN 'don't check nodes for mutual L GOTO 140 END IF ' IF ib(i, 2) > nnode OR ib(i, 3) > nnode THEN PRINT " node number higher than maximum specified" PRINT " first node was read as "; ib(i, 2) PRINT " second node was read as "; ib(i, 3) CALL dataerror END IF ' 140 f2$ = " +#.###^^^^ " SELECT CASE ib(i, 4) 'end of select is at 490 ' CASE 1 a$ = "1 -R" CALL printfour(i, ib(), b(), a$, g$) 'print first 5 elements of data ' CASE 2 a$ = "2 -L" CALL printfour(i, ib(), b(), a$, g$) b(i, 4) = b(i, 1) 'IF part of a L&M pair, this value is nls = nls + 1 'modified by matrix inversion ib(i, 5) = nls y(nls, nls) = b(i, 1) tag(nls, nls) = i ' CASE 3 a$ = "3 -C" CALL printfour(i, ib(), b(), a$, g$) ' CASE 4 a$ = "4 -SRL" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2) 'read and print l PRINT USING f2$; b(i, 2); ' CASE 5 a$ = "5 -SRC" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2) 'read and print c PRINT USING f2$; b(i, 2); ' CASE 6 a$ = "6 -PRL" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2) 'read and print l PRINT USING f2$; b(i, 2); ' CASE 7 a$ = "7 -PRC" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2) 'read and print c PRINT USING f2$; b(i, 2); ' CASE 8 a$ = "8 -M" CALL printfour(i, ib(), b(), a$, g$) nlm = nlm + 1 ' CASE 9 a$ = "9 -V dep R" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2) 'input symmetry flag for v dep r IF b(i, 2) = 1 THEN a$ = "symmetrical about zero" ELSE a$ = "not symmetrical" END IF PRINT : PRINT PRINT a$ PRINT " start on segment "; b(i, 1) PRINT : PRINT PRINT " volts amps" CALL readnonlin(i, b(), nnl, cn()) 'read v and i ' CASE 10 a$ = "10-T dep R" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2), b(i, 3), b(i, 4)'descriptors for tdep r PRINT : PRINT PRINT USING f$; " start on segment "; b(i, 1) PRINT USING f$; " time to start changing "; b(i, 2) PRINT USING f$; " voltage to start changing "; b(i, 3) PRINT " control branch "; b(i, 4) PRINT : PRINT PRINT " seconds ohms" CALL readnonlin(i, b(), nnl, cn()) 'read t and r ' CASE 11 a$ = "11-W dep R" CALL printfour(i, ib(), b(), a$, g$) PRINT : PRINT PRINT " start on segment "; b(i, 1) PRINT : PRINT PRINT " watt-sec ohms" CALL readnonlin(i, b(), nnl, cn()) ' CASE 12 a$ = "12-I dep L" CALL printfour(i, ib(), b(), a$, g$) PRINT : PRINT PRINT " start on segment "; b(i, 1) PRINT : PRINT PRINT " weber-turns amps" CALL readnonlin(i, b(), nnl, cn()) 'read flux and i ' CASE 13 a$ = "13-Switch" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2), b(i, 3) ' nsw = nsw + 1 IF nsw > 6 THEN PRINT "error - only 6 switches allowed" PRINT "value was read as "; nsw CALL dataerror END IF ' g(nsw, 1) = ib(i, 1) 'transfer branch number g(nsw, 2) = b(i, 3) 'transfer state g(nsw, 3) = b(i, 1) 'transfer ron g(nsw, 4) = b(i, 2) 'transfer roff ' IF b(i, 3) = 0 THEN b(i, 4) = b(i, 1) 'set initial state to ron a$ = " initially closed " ELSEIF b(i, 3) = 1 THEN b(i, 4) = b(i, 2) a$ = " initially open " ELSE PRINT "state may be defined wrong" PRINT "state was read as "; b(i, 3) CALL dataerror END IF ' PRINT FOR jj = 5 TO 9 'read remaining part of switch data INPUT #1, g(nsw, jj) NEXT jj PRINT f4$ = " & +#.###^^^^ " PRINT USING f4$; " r on "; g(nsw, 3) PRINT USING f4$; " r off "; g(nsw, 4) PRINT USING f4$; " t for close "; g(nsw, 5) PRINT USING f4$; " v for close "; g(nsw, 6) PRINT USING f4$; " t for open "; g(nsw, 7) PRINT USING f4$; " t for part "; g(nsw, 8) PRINT USING f4$; " I for opening "; g(nsw, 9) PRINT a$ ' CASE 14 a$ = "14-Diode" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2), b(i, 3), b(i, 4) ' nsc = nsc + 1 IF nsc > 6 THEN PRINT "error - only 6 semiconductor switches allowed "; "" PRINT "value was read as "; nsc CALL dataerror END IF ' g1(nsc, 1) = ib(i, 1) 'transfer branch number g1(nsc, 2) = b(i, 4) 'initial state g1(nsc, 3) = b(i, 1) 'transfer ron or forward bias r g1(nsc, 4) = b(i, 2) 'transfer roff or reverse bias r g1(nsc, 5) = b(i, 3) 'transfer voltage at intersection IF b(i, 4) = 0 THEN b(i, 6) = b(i, 1) a$ = " initially conducting " ELSEIF b(i, 4) = 1 THEN b(i, 6) = b(i, 2) a$ = " initially not conducting " ELSE PRINT "initial state may be defined wrong" PRINT "state was read as "; b(i, 4) CALL dataerror END IF ' PRINT f4$ = " & +#.###^^^^ " PRINT USING f4$; " r on "; g1(nsc, 3) PRINT USING f4$; " r off "; g1(nsc, 4) PRINT USING f4$; " voltage break "; g1(nsc, 5) PRINT a$ ' CASE 15 a$ = "15-Thyris" CALL printfour(i, ib(), b(), a$, g$) INPUT #1, b(i, 2), b(i, 3), b(i, 4) INPUT #1, b(i, 5) ' nsc = nsc + 1 IF nsc > 6 THEN PRINT "error - only 6 semiconductor switches allowed "; "" PRINT "value was read as "; nsc CALL dataerror END IF ' g1(nsc, 1) = ib(i, 1) 'transfer branch number g1(nsc, 2) = b(i, 5) 'transfer initial state g1(nsc, 3) = b(i, 1) 'transfer ron or conducting r g1(nsc, 4) = b(i, 2) 'transfer roff or non-conducting r g1(nsc, 5) = b(i, 3) 'transfer trigger and holding current g1(nsc, 6) = b(i, 4) 'transfer control branch IF b(i, 5) = 0 THEN b(i, 6) = b(i, 1) a$ = "initially conducting " ELSEIF b(i, 5) = 1 THEN b(i, 6) = b(i, 2) a$ = "initially not conducting " ELSE PRINT "initial state may be defined wrong" PRINT "initial state was read as "; b(i, 5) CALL dataerror END IF ' PRINT f4$ = " & +#.###^^^^ " PRINT USING f4$; " r on "; g1(nsc, 3) PRINT USING f4$; " r off "; g1(nsc, 4) PRINT USING f4$; " trigger I "; g1(nsc, 5) PRINT USING f4$; " control branch "; g1(nsc, 6) PRINT a$ ' CASE 16 a$ = "16-T line" CALL printfour(i, ib(), b(), a$, g$) ntl = ntl + 1 IF ntl > maxntl + 1 THEN PRINT "error - only"; maxntl; " transmission lines allowed" PRINT "value was read as "; ntl CALL dataerror END IF ' INPUT #1, b(i, 2), b(i, 3), b(i, 4) v1 = b(i, 1) 'surge impedance v2 = b(i, 2) 'propagation velocity v3 = b(i, 3) 'length v4 = .25 * b(i, 4) '1/4 of series resistance f4$ = " & +#.###^^^^ " PRINT PRINT USING f4$; " surge impedance "; b(i, 1) PRINT USING f4$; " propagation velocity "; b(i, 2) PRINT USING f4$; " length "; b(i, 3) PRINT USING f4$; " series resistance "; b(i, 4) l9(ntl, 0) = ib(i, 1) 'branch number l9(ntl, 1) = ib(i, 2) 'from node l9(ntl, 2) = ib(i, 3) 'to node l9(ntl, 3) = v1 l9(ntl, 4) = v3 * .3048 / (3E+08 * v2) 'length/velocity l9(ntl, 7) = v4 h = (v1 - v4) / (v1 + v4) l9(ntl, 8) = (1 + h) / 2 l9(ntl, 9) = (1 - h) / 2 l9(ntl, 10) = 1 / (v1 + v4) 'admittance of line q9 = INT((l9(ntl, 4) / t1) + .5) 'calculate required line storage IF q9 > nlinestep THEN PRINT : PRINT "caution--this line requires "; q9; "elements and the " PRINT "assigned space is only 200 elements. The line has been shortened" PRINT "to "; 200 / q9; " of the original length" q9 = 200 l9(ntl, 4) = q9 END IF q(ntl) = 1 q1(ntl) = q9 ' END SELECT ' PRINT IF g$ <> "g" THEN CALL examine(g$) END IF ' 150 NEXT i ' '--------- modify inductors if mutuals are involved ------- ' IF nlm > 0 THEN CALL modifyL(nbran, nls, ib(), b(), y(), tag()) END IF ' ' --------read source data----- ' e1$ = "### ### ### +#.###^^^^ " e2$ = "\ \ +#.###^^^^ " ' i = 0 'i counts all types of sources nitab = 0 'itab counts tabular sources nidep = 0 'idep counts dependent sources ' PRINT : PRINT "source data" PRINT PRINT "branch type shape start magnitude value 3 value 4" PRINT ' 160 i = i + 1 INPUT #1, e1(i, 1) IF e1(i, 1) = -99 THEN 200 'go on to initial conditions ' IF i > 16 THEN PRINT "too many sources (limit is 16 )" PRINT "value was read as "; i CALL dataerror END IF ' INPUT #1, e1(i, 2), e1(i, 3), e1(i, 4) 'rest of input is at 170 IF e1(i, 2) < 1 OR e1(i, 2) > 3 THEN PRINT "Error in type of source " PRINT "your entry was read as type "; e1(i, 2) CALL dataerror END IF ' SELECT CASE e1(i, 2) ' CASE 1 s1$ = "independent voltage" CASE 2 s1$ = "independent current" CASE 3 s1$ = "dependent current " nidep = nidep + 1 ' END SELECT ' IF e1(i, 2) < 3 THEN GOTO independent ELSE GOTO dependent END IF ' independent: ' IF e1(i, 3) < 1 OR e1(i, 3) > 10 THEN PRINT "Error in waveshape type " PRINT "your entry was read as type "; e1(i, 3) CALL dataerror END IF SELECT CASE e1(i, 3) ' CASE 1 s2$ = "trapezoidal " CASE 2 s2$ = "triangular " CASE 3 s2$ = "double exponential " CASE 4 s2$ = "modified exponential " CASE 5 s2$ = "sine squared " CASE 6 s2$ = "offset sine " CASE 7 s2$ = "exponentially increasing sine " CASE 8 s2$ = "decaying sinusoid " CASE 9 s2$ = "decaying cosinusoid " CASE 10 s2$ = "tabular " END SELECT ' PRINT "branch "; e1(i, 1); " "; s1$; " - "; s2$ ' 170 INPUT #1, e1(i, 5), e1(i, 6), e1(i, 7) ' PRINT USING e2$; " starting time "; e1(i, 4) PRINT USING e2$; " magnitude "; e1(i, 5) PRINT USING e2$; " defining time 1 "; e1(i, 6) PRINT USING e2$; " defining time 2 "; e1(i, 7) ' IF e1(i, 7) = 0 THEN e1(i, 7) = t1 'guard against zero front time IF e1(i, 3) = 10 THEN GOTO tabsorce ELSE GOTO examine END IF ' examine: ' IF g$ <> "g" THEN CALL examine(g$) END IF GOTO 160 'read next source ' tabsorce: ' ---------list data on tabular source------- ' PRINT : PRINT "data on tabular source": PRINT PRINT " "; "time", "amplitude" f$ = " +#.####^^^^ +#.####^^^^" ' nitab = nitab + 1 IF nitab > 1 THEN PRINT "error - only 1 tabular source allowed" PRINT "your entry was read as "; nitab CALL dataerror END IF j = 1 180 INPUT #1, t3(j, 1) 'check for end of tabular source IF t3(j, 1) = -99 THEN 160 'read next source INPUT #1, t3(j, 2) PRINT USING f$; t3(j, 1); t3(j, 2) j = j + 1 ' IF g$ <> "g" THEN CALL examine(g$) END IF ' GOTO 180 'read next line of tabular source ' ' dependent: ' PRINT "across branch "; e1(i, 1); " "; s1$ PRINT " control branch "; e1(i, 3) PRINT USING e2$; " transconductance "; e1(i, 4) ' IF g$ <> "g" THEN CALL examine(g$) END IF 190 GOTO 160 'read next source ' ' 200 nsorce = i - 1 '------------ tag E1 table to identify control branches for Idep sorce FOR i = 1 TO nsorce IF e1(i, 2) = 3 THEN FOR j = 1 TO nsorce IF e1(j, 1) = e1(i, 3) THEN e1(i, 5) = j 'j = line in table where control branch is located END IF NEXT j END IF NEXT i ' -------- read initial conditions -------- ' PRINT : PRINT "initial conditions": PRINT PRINT PRINT "branch", "type", "value": PRINT 205 INPUT #1, iloc 'read location for initial condition IF iloc = -99 THEN 210 INPUT #1, itype, value f$ = "#### #### +#.####^^^^" PRINT USING f$; iloc; itype; value IF itype < 1 OR itype > 2 THEN PRINT "Error in type - your data was read as "; itype CALL dataerror END IF IF ib(iloc, 4) = 16 THEN i = 1 'find location in L9 table, starting at top of table 206 IF l9(i, 0) = iloc THEN ilocline = 1 ELSE i = i + 1 GOTO 206 END IF 'end of locline routine ' q9 = q1(ilocline) 'length of line in delta t IF itype = 1 THEN FOR i = 1 TO q9 'spread voltage out along t line l1(i, ilocline) = value / 2 'load the forward voltages l2(i, ilocline) = value / 2 'load the reverse voltages NEXT i ELSE FOR i = 1 TO q9 'spread current out along t line l3(i, ilocline) = value / 2 'load the forward currents l4(i, ilocline) = value / 2 'load the reverse currents NEXT i END IF ELSE IF itype = 1 THEN e3(iloc, 2) = -value 'initial voltage into e3(-,2) ELSE e3(iloc, 4) = value 'initial current into e3(-,4) END IF END IF GOTO 205 'loop for next initial condition 210 PRINT ' '------- read locations for output ---------- ' PRINT : PRINT "Output will be at the following points" PRINT INPUT #1, noutp WRITE #2, noutp IF noutp > noutpmax THEN PRINT "Too many output points (max = "; noutpmax; ") or error in output type" PRINT "your entry was read as "; noutp CALL dataerror END IF ' FOR i = 1 TO noutp INPUT #1, o$(i) WRITE #2, o$(i) PRINT " point "; i; " will be "; o$(i) p$ = MID$(o$(i), 1, 2) q$ = MID$(o$(i), 3, 2) outplace%(i) = VAL(q$) SELECT CASE p$ ' CASE "NV", "nv" outtype%(i) = 1 CASE "BV", "bv" outtype%(i) = 2 CASE "BI", "bi" outtype%(i) = 3 CASE "BP", "bp" outtype%(i) = 4 CASE "BE", "be" outtype%(i) = 5 CASE "SN", "sn" outtype%(i) = 6 CASE "SB", "sb" outtype%(i) = 7 CASE "AI", "ai" outtype%(i) = 8 END SELECT ' NEXT i ' CLOSE #1 PRINT PRINT "Press any key to commence calculations" DO WHILE INKEY$ = "": LOOP ' '-------------------------------------------- ' 285 CLS PRINT : PRINT ' LOCATE 1, 1: PRINT " time " LOCATE 2, 1: PRINT " point " FOR i = 1 TO noutp LOCATE 2 + i, 1: PRINT i; LOCATE 2 + i, 5: PRINT o$(i); 'print names of output points NEXT i LOCATE 2, 40: PRINT " (use f1 to break)" ' LOCATE 23, 76: PRINT "done"; LOCATE 24, 80: PRINT "|"; ' 290 ' --------start calculation loop--------- ' ' -------initialize the non-linear branch admittances----- FOR i = 1 TO nbran IF ib(i, 4) < 9 OR ib(i, 4) > 12 THEN 292 CALL nonlinZ(i, ib(), b(), e3(), ci3(), cn(), t, t1, K9) 'to non-linear impedance routine 292 NEXT i ' ' t = -t1 'initialize time to -1 delta t because 't is incremented before loop ' ' 300 ' <<<---<<--- entry point for next time step if topology changes ' (end of loop is at 410) icount = 1 nn = nnode 'set dimensions for inversion CALL buildy(i, nn, nbran, nsorce, nidep, ib(), b(), t1, y(), e1()) 'goto "formulate matrix " routine" CALL invertH(nn, y()) 'invert Y matrix ' 310 ' <<<---<<<---<<<---<<<------- entry point if no changes to Y matrix ' K9 = 0 'k9 is a flag that is changed to k9=1 if circuit topology changes ' t = t + t1 'increment time IF icount = increm THEN jout = jout + 1 'increment position in ' output table every nth step ' 'check f1 key to see if break is required ON KEY(1) GOSUB breakit: ' ' ----- calculate instantaneous values of sources ' i = 1 WHILE i <= nsorce IF e1(i, 2) = 3 THEN 316 'bypass if dependent source j = e1(i, 1) 'location of source ci3(j, 1) = 0 ' ci3(j, 2) = 0 ' ix2 = e1(i, 2) 'type of source ix3 = e1(i, 3) 'waveshape v4 = e1(i, 4) 'start time v5 = e1(i, 5) 'magnitude v6 = e1(i, 6) 'defining value 1 v7 = e1(i, 7) 'defining value 2 CALL waveshape(t, t3(), ix3, v4, v5, v6, v7, vmag) e1(i, 8) = vmag 316 i = i + 1 WEND ' ' -- put source data into ci3 matrix and calc branch i FOR i = 1 TO nsorce j = e1(i, 1) 'branch where source is found SELECT CASE e1(i, 2) CASE 1 'independent branch voltage ci3(j, 2) = ci3(j, 2) + e1(i, 8) CASE 2 'independent branch current ci3(j, 1) = ci3(j, 1) + e1(i, 8) CASE 3 'dependent source treated later ' END SELECT NEXT i 320 FOR i = 1 TO nbran 'calculate total independent currents ci3(i, 5) = ci3(i, 1) - b(i, 5) * ci3(i, 2) + ci3(i, 3) ' source I = ind I - yb * eb + non-lin r inter i NEXT i ' IF nsorce = 0 OR nidep = 0 THEN 330 ' add the contribution of the independent voltage (Eb) in the ' control branch to the current in the controlled branch FOR i = 1 TO nsorce IF e1(i, 2) = 3 THEN j = e1(i, 1) 'branch across which the source is located k = e1(i, 5) 'line in E1 table where Eb is located IF e1(k, 2) = 1 THEN ci3(j, 5) = ci3(j, 5) - e1(i, 4) * e1(k, 5) ' source I = source I - gm * eb END IF END IF NEXT i ' 330 ' --------calculate history terms for branches---- ' FOR i = 1 TO nbran ' SELECT CASE ib(i, 4) ' CASE 2 'l ci3(i, 6) = b(i, 5) * e3(i, 2) + e3(i, 4) ' CASE 3 'c ci3(i, 6) = -b(i, 5) * e3(i, 2) - e3(i, 4) ' CASE 4 'srl ci3(i, 6) = b(i, 1) * .5 * t1 / b(i, 2) ci3(i, 6) = (1 - ci3(i, 6)) / (1 + ci3(i, 6)) ci3(i, 6) = b(i, 5) * e3(i, 2) + e3(i, 4) * ci3(i, 6) ' CASE 5 'src ci3(i, 6) = b(i, 1) * 2 * b(i, 2) / t1 ci3(i, 6) = (1 - ci3(i, 6)) / (1 + ci3(i, 6)) ci3(i, 6) = -b(i, 5) * e3(i, 2) - e3(i, 4) * ci3(i, 6) ' CASE 6 'prl parta = .5 * t1 / b(i, 2) - 1 / b(i, 1) ci3(i, 6) = parta * e3(i, 2) + e3(i, 4) ' CASE 7 'prc parta = 2 * b(i, 2) / t1 - 1 / b(i, 1) ci3(i, 6) = -parta * e3(i, 2) - e3(i, 4) ' CASE 8 'mut l j = ib(i, 2) 'branch number of Lp k = ib(i, 3) 'branch number of Ls ci3(j, 6) = ci3(j, 6) + b(i, 5) * e3(k, 2) ci3(k, 6) = ci3(k, 6) + b(i, 5) * e3(j, 2) ' CASE 12 'sat l bphik = t1 / (2 * b(i, 5)) aphik = ci3(i, 4) ci3(i, 6) = (e3(i, 8) - aphik) / bphik + b(i, 5) * e3(i, 2) ' CASE ELSE 'all resistive branches END SELECT ' NEXT i '--------------- calculate history terms for transmission lines -- i = 1 WHILE i <= ntl c1 = l2(q(i), i) * l9(i, 10) + l4(q(i), i) c2 = l1(q(i), i) * l9(i, 10) + l3(q(i), i) l9(i, 5) = l9(i, 8) * c1 + l9(i, 9) * c2 l9(i, 6) = l9(i, 8) * c2 + l9(i, 9) * c1 i = i + 1 WEND ' 340 ' ---form nodal current vector with r,l,c elements ' FOR i = 1 TO nnode: ci4(i, 1) = 0: NEXT i ' clear ci4 vector ' FOR i = 1 TO nbran ci9 = ci3(i, 5) - ci3(i, 6) j = ib(i, 2) 'from node k = ib(i, 3) 'to node ci4(j, 1) = ci4(j, 1) + ci9 ci4(k, 1) = ci4(k, 1) - ci9 NEXT i ' ' ---add transmission lines to nodal current vector---- ' i = 1 WHILE i <= ntl j = l9(i, 1) 'from node k = l9(i, 2) 'to node ci4(j, 1) = ci4(j, 1) + l9(i, 5) ci4(k, 1) = ci4(k, 1) + l9(i, 6) i = i + 1 WEND ' ' -------multiply Y and I vector and so get------ ' nodal voltage vector ' FOR i = 1 TO nnode e2(i, 1) = 0 FOR j = 1 TO nnode IF ci4(j, 1) = 0 THEN 341 e2(i, 1) = e2(i, 1) + y(i, j) * ci4(j, 1) 341 NEXT j NEXT i ' ' ----------- calculate branch voltages ----------- ' FOR i = 1 TO nbran j = ib(i, 2) 'from node k = ib(i, 3) 'to node IF j = 0 THEN e7 = 0 ELSE e7 = e2(j, 1) 'voltage at from node IF k = 0 THEN e8 = 0 ELSE e8 = e2(k, 1) 'voltage at to node e3(i, 1) = e7 - e8 + ci3(i, 2) ' ci3(i, 2) = indepen branch voltage NEXT i ' ' --- calculate branch currents -------- ' FOR i = 1 TO nbran e3(i, 3) = 0 ' SELECT CASE ib(i, 4) ' CASE 2, 3, 4, 5, 6, 7 ' e3(i, 3) = e3(i, 3) + e3(i, 1) * b(i, 5) + ci3(i, 6) ' i + eb * yb + iold ' CASE 1, 9, 10, 11, 12, 13, 14, 15 ' e3(i, 3) = e3(i, 3) + e3(i, 1) * b(i, 5) - ci3(i, 3) ' i + eb * yb - iinter ' CASE 8 j = ib(i, 2) 'add contribution of mutual inductance k = ib(i, 3) e3(j, 3) = e3(j, 3) + b(i, 5) * e3(k, 1) e3(k, 3) = e3(k, 3) + b(i, 5) * e3(j, 1) ' CASE 16 FOR j = 1 TO maxntl 'find appropriate line in l9 table IF l9(j, 0) = i THEN 350 NEXT j ' 350 ii = l9(j, 1) 'number of starting node e3(i, 3) = e2(ii, 1) * l9(j, 10) - l9(j, 5) ' e2(ii, 1) = node voltage at sending end of line ' l9(j, 5) = I-old at sending end ' e3(i, 3) = I at sending end 'use l9(j, 2) and l9(j, 6) for current at receiving end of line ' END SELECT ' NEXT i ' ' -------------- calculate V and I to go into transmission lines ---- ' i = 1 WHILE i <= ntl l1(q(i), i) = e2(l9(i, 1), 1) 'voltage into from node l2(q(i), i) = e2(l9(i, 2), 1) 'voltage into to node l3(q(i), i) = e2(l9(i, 1), 1) / l9(i, 3) - l9(i, 5) 'current into from node l4(q(i), i) = e2(l9(i, 2), 1) / l9(i, 3) - l9(i, 6) 'current into to node q(i) = q(i) + 1 IF q(i) < q1(i) + 1 THEN GOTO 358 'controls indexing of waves through lines END IF q(i) = 1 358 i = i + 1 WEND ' 360 FOR i = 1 TO nbran 'calculate branch powers ' SELECT CASE ib(i, 4) ' CASE 1, 9, 10, 11, 13, 14, 15, 16 'all resistive branches ' e3(i, 5) = e3(i, 1) * e3(i, 3) ' CASE 2, 8, 12 'all inductive branches ' parta = e3(i, 1) * b(i, 5) * 2 / t1 e3(i, 5) = e3(i, 1) * (parta + b(i, 4)) / 2 ' CASE 3 'c ' e3(i, 5) = e3(i, 3) * e3(i, 3) * b(i, 1) / 2 ' CASE 4 'srl ' e3(i, 5) = e3(i, 1) * e3(i, 1) * (b(i, 1) + b(i, 2) / 2) ' CASE 5 'src ' e3(i, 5) = e3(i, 1) * e3(i, 1) * b(i, 1) voncap = e3(i, 3) - e3(i, 1) * b(i, 1) e3(i, 5) = e3(i, 5) + b(i, 2) * voncap * voncap / 2 ' CASE 6 'prl ' e3(i, 5) = e3(i, 3) * e3(i, 3) / b(i, 1) curinl = e3(i, 1) - e3(i, 3) / b(i, 1) e3(i, 5) = e3(i, 5) + b(i, 2) * curinl * curinl / 2 ' CASE 7 'prc ' parta = 1 / b(i, 1) + b(i, 2) / 2 e3(i, 5) = e3(i, 3) * e3(i, 3) * parta ' END SELECT ' e3(i, 6) = e3(i, 6) + t1 * e3(i, 5) 'branch energy ' NEXT i ' ' ----calculate new integrals------- ' IF t = 0 THEN 370 'integrals = 0 at t = 0 FOR i = 1 TO nnode e2(i, 3) = e2(i, 3) + .5 * t1 * (e2(i, 1) + e2(i, 2)) NEXT i FOR i = 1 TO nbran e3(i, 7) = e3(i, 7) + .5 * t1 * (e3(i, 1) + e3(i, 2)) e3(i, 8) = e3(i, 8) + .5 * t1 * (e3(i, 3) * e3(i, 3) + e3(i, 4) * e3(i, 4)) NEXT i ' ' fout$ = " +#.###^^^^ " IF icount = increm THEN 370 ELSE 380 370 '-------- output from this time step ---- ' LOCATE 1, 20 PRINT USING " +#.###^^^^ "; t WRITE #2, t ' FOR i = 1 TO noutp 'PRINT "outtype and outplace are "; outtype%, outplace% 'INPUT a$ SELECT CASE outtype%(i) ' CASE 1 outdata(i) = e2(outplace%(i), 1) 'nv=node voltage CASE 2 outdata(i) = e3(outplace%(i), 1) 'bv=branch voltage CASE 3 outdata(i) = e3(outplace%(i), 3) 'bi=branch current CASE 4 outdata(i) = e3(outplace%(i), 5) 'bp=branch power CASE 5 outdata(i) = e3(outplace%(i), 6) 'be=branch energy CASE 6 outdata(i) = e2(outplace%(i), 3) 'sn=integral of node voltage CASE 7 outdata(i) = e3(outplace%(i), 7) 'sb=integral of branch voltage CASE 8 outdata(i) = e3(outplace%(i), 8) 'ai = action integral END SELECT LOCATE 2 + i, 20 PRINT USING " +#.###^^^^ "; outdata(i) WRITE #2, outdata(i) 'output to file ' NEXT i ' 380 IF icount = increm THEN icount = 1 ELSE icount = icount + 1 ' ' ----- update the "amount done" line putit = 1 + 80 * t / t2 IF putit > 80 THEN putit = 80 LOCATE 24, putit PRINT "-"; ' --- store newly calculated values in old columns --- ' FOR i = 1 TO nnode e2(i, 2) = e2(i, 1) NEXT i FOR i = 1 TO nbran e3(i, 2) = e3(i, 1) e3(i, 4) = e3(i, 3) ' ci3(i,7) = ci3(i,6) NEXT i ' ' --------check each switch to see if it should change state----- ' i = 1 WHILE i <= nsw ksw = g(i, 1) 'branch where switch is located state = g(i, 2) '0 = closed, 1 = open ron = g(i, 3) roff = g(i, 4) tclose = g(i, 5) vbo = g(i, 6) topen = g(i, 7) tpart = g(i, 8) cmar = g(i, 9): tot = topen + tpart volts = e3(ksw, 1) K7 = g(i, 2) 'note state of switch prior to checking for any change IF t >= topen AND t < tot THEN vb = vbo * (t - topen) / tpart ELSE vb = vbo END IF ' IF g(i, 3) = 0 THEN 'what is present state? GOTO state0 ELSE GOTO state1 END IF ' state0: '----- state 0 - initially closed ---- IF t < topen THEN GOTO closeit 'leave switch closed ELSEIF t < tclose THEN GOTO checkiv 'check for opening ELSE GOTO closeit 'close the switch END IF ' state1: '----- state 1 - initially open ---- IF t < tclose AND g(i, 2) = 1 THEN GOTO checkvbo 'check for closure on vbo ELSEIF t < topen THEN GOTO closeit 'close the switch ELSE GOTO checkiv 'check for opening END IF ' checkvbo: ' ---- check for breakdown on voltage ---- IF volts >= vbo THEN GOTO closeit 'close the switch ELSE GOTO openit 'leave switch open END IF ' checkiv: '---- check to see if switch can open --- IF g(i, 2) = 0 AND AMPS < cmar THEN GOTO openit 'open the switch ELSEIF g(i, 2) = 1 AND volts >= vb THEN GOTO closeit 'close the switch ELSE GOTO checkchange1 END IF ' closeit: b(ksw, 4) = ron 'switch is closed g(i, 2) = 0 a$ = "closed " GOTO checkchange1 ' openit: b(ksw, 4) = roff 'switch is open g(i, 2) = 1 a$ = "open " GOTO checkchange1 ' checkchange1: ' ----- check to see if switch has changed state ---- IF K7 <> g(i, 2) THEN K9 = 1 LOCATE 14 + i, 1: PRINT "branch "; ksw; " is "; a$ i = i + 1 'check next switch WEND ' ' ---check diodes and thyristors for changes--- i = 1 WHILE i <= nsc 'WEND is after checkchange2 ksc = g1(i, 1) 'branch where switch is located state = g1(i, 2) '0 =closed, 1 = open K7 = state 'state prior to checking volts = e3(ksc, 1) 'voltage across branch AMPS = e3(ksc, 3) 'current through branch ron = g1(i, 3) roff = g1(i, 4) vinter = g1(i, 5) trigi = e3(g1(i, 6), 3) 'current in control branch holdi = g1(i, 5) ' IF state = 0 THEN GOTO closed2 ELSE GOTO open2 END IF ' closed2: ' ----- conducting - check for opening ------- IF ib(ksc, 4) = 14 THEN IF volts < vinter THEN GOTO setoff ELSE GOTO checkchange2 'check diodes END IF ELSE IF volts < 0 OR AMPS < holdi THEN GOTO setoff ELSE GOTO checkchange2 'check thyristors END IF END IF open2: ' ----- not conducting - check for closing ------- IF ib(ksc, 4) = 14 THEN IF volts >= vinter THEN GOTO seton ELSE GOTO checkchange2 'check diodes END IF ELSE IF volts > 0 AND trigi >= holdi THEN GOTO seton ELSE GOTO checkchange2 'check thyristors END IF END IF ' seton: b(ksc, 6) = ron 'set to on resistance g1(i, 2) = 0 GOTO checkchange2 ' setoff: b(ksc, 6) = roff 'set to off resistance g1(i, 2) = 1 GOTO checkchange2 ' checkchange2: ' ----- check to see if semiconductor switch has changed state --- IF K7 <> g1(i, 2) THEN K9 = 1 'k9=1 IF state changed ' i = i + 1 WEND ' ---check each non-linear element for a change--- ' ' FOR i = 1 TO nbran IF ib(i, 4) < 9 OR ib(i, 4) > 12 THEN 390 'no gosub if not non-linear CALL nonlinZ(i, ib(), b(), e3(), ci3(), cn(), t, t1, K9) ' go to non-linear routine 390 NEXT i ' ---set k9=1 if any non-linear element has changed slope--- ' FOR i = 2 TO 16 STEP 2 IF cn(2, i) = 1 THEN 400 K9 = 1 400 NEXT i ' ' -----increment the time and repeat---- ' IF t - t1 > t2 THEN GOTO 420 'to output routine 'IF jout >= nmax THEN ' PRINT "output array is full " ' BEEP: BEEP: BEEP ' GOTO 420 'END IF ' IF k6 = 1 THEN 410 'don't check times after changing delta T IF t >= t3 THEN t1 = t4 k6 = 1 K9 = 1 END IF ' 410 IF K9 = 0 THEN 310 'do not reformulate y matrix IF K9 = 1 THEN 300 'do reformulate the y matrix ' 420 CLOSE #2 CHAIN "tranplot" END ' runerror: ' PRINT "A type "; ERR; " error has occurred" PRINT "somewhere after line "; ERL PRINT "The key to the errors is listed in your manual" PRINT "press any key to end the session" SLEEP (10) INPUT a$ 490 END breakit: GOTO 420 '*********************************************** ' '**************************************************** ' SUBROUTINES '***************************************************** SUB addtoY (i, ib(), y1, j1, k1, j2, k2, y()) ' '------subroutine to add to y matrix---- ' y(j1, j2) = y(j1, j2) + y1 y(k1, k2) = y(k1, k2) + y1 IF ib(i, 4) <> 16 THEN y(j1, k2) = y(j1, k2) - y1 y(k1, j2) = y(k1, j2) - y1 END IF END SUB ' '********************************************************* SUB buildy (i, nn, nbran, nsorce, nidep, ib(), b(), t1, y(), e1()) ' ------formulate y matrix----- ' FOR i = 1 TO nn 'clear Y matrix FOR j = 1 TO nn y(i, j) = 0 NEXT j NEXT i ' FOR i = 1 TO nbran 'calculate admittances ' SELECT CASE ib(i, 4) CASE 1 'resistor ' b(i, 5) = 1 / b(i, 1) ' CASE 2, 8 'inductor or mutual ' b(i, 5) = .5 * t1 / b(i, 4) ' CASE 3 'capacitor ' b(i, 5) = 2 * b(i, 1) / t1 ' CASE 4 'srl ' b(i, 5) = .5 * t1 / b(i, 2) b(i, 5) = b(i, 5) / (1 + b(i, 4) * b(i, 1)) ' CASE 5 'src ' b(i, 5) = 2 * b(i, 2) / t1 b(i, 5) = b(i, 5) / (1 + b(i, 5) * b(i, 1)) ' CASE 6 'prl ' b(i, 5) = .5 * t1 / b(i, 2) + 1 / b(i, 1) ' CASE 7 'prc ' b(i, 5) = 2 * b(i, 2) / t1 + 1 / b(i, 1) ' CASE 13 'switch ' b(i, 5) = 1 / b(i, 4) ' CASE 14, 15 'diode or scr ' b(i, 5) = 1 / b(i, 6) ' CASE 9, 10, 11, 12 'all four non-linear types ' 'non-linear impedances determined in "check for changes " routine" ' CASE 16 ' b(i, 5) = 1 / (b(i, 1) + b(i, 4) / 4) 't-line ' END SELECT NEXT i ' ' 'beeps every time the admittance array is calculated FOR i = 1 TO nbran y1 = b(i, 5) IF ib(i, 4) = 8 THEN jj = ib(i, 2) 'branch number of first L kk = ib(i, 3) 'branch number of second L ' j1 = ib(jj, 2) k1 = ib(jj, 3) 'nodes for mutual L j2 = ib(kk, 2) k2 = ib(kk, 3) CALL addtoY(i, ib(), y1, j1, k1, j2, k2, y()) ' j1 = ib(kk, 2) k1 = ib(kk, 3) j2 = ib(jj, 2) k2 = ib(jj, 3) CALL addtoY(i, ib(), y1, j1, k1, j2, k2, y()) ELSE j1 = ib(i, 2) k1 = ib(i, 3) 'nodes for normal branch j2 = j1 k2 = k1 CALL addtoY(i, ib(), y1, j1, k1, j2, k2, y()) END IF NEXT i ' ' ---------- add gm ----------- ' IF nidep > 0 THEN FOR i = 1 TO nsorce IF e1(i, 2) = 3 THEN j1 = ib(e1(i, 1), 2) k1 = ib(e1(i, 1), 3) j2 = ib(e1(i, 3), 2) k2 = ib(e1(i, 3), 3) y1 = e1(i, 4) 'y1 = gm (transconductance) ii = ib(e1(1, 1), 1) 'set i as used in addtoY routine to 'number of controlling branch. AddtoY uses i only to 'check whether it is getting a transmission line so lets 'send it a valid branch type (ii), NOT the source number CALL addtoY(ii, ib(), y1, j1, k1, j2, k2, y()) END IF NEXT i END IF ' FOR i = 1 TO nn 'check for disconnected nodes IF y(i, i) = 0 THEN PRINT "There is an error in the admittance matrix" PRINT "Perhaps there is a node not connected" PRINT "check node "; i PRINT "press any key to continue" INPUT a$ END IF NEXT i ' END SUB 'end of formulate y matrix routine '**************************************** ' SUB dataerror PRINT "enter any character to erase the screen" DO: LOOP WHILE INKEY$ = "" SLEEP (10) STOP END SUB '*********************************** SUB examine (g$) ' 1000 a$ = INKEY$: IF a$ = "" THEN 1000 'stop to examine the line IF a$ = "G" OR a$ = "g" THEN g$ = "g" 'suppress further stops END SUB '************************************************************ SUB invertH (nn, y()) 'matrix inversion routine based on Hornbeck ' "Numerical Methods, p 295" ' Quantum Publishers ' New York, 1975 ' 'This is a Gauss-Jordan elimination routine. The matrix is 'called in y(nn , nn) and is returned in y(nn , nn). The 'original matrix is consumed. ' 'DEFINT I-N DIM j(nn + 20) ' pd = 1 ' FOR ii = 1 TO nn dd = 0 ' FOR kk = 1 TO nn dd = dd + y(ii, kk) * y(ii, kk) NEXT kk ' dd = SQR(dd) ' PRINT ii, dd, pd: INPUT debug$ ' pd = pd * dd NEXT ii ' detm = 1 ' FOR ii = 1 TO nn j(ii + 20) = ii NEXT ii ' FOR ii = 1 TO nn cc = 0 mm = ii ' FOR kk = ii TO nn 'NOTE ii (lower case L), not 1 (one) IF (ABS(cc) - ABS(y(ii, kk))) >= 0 THEN 1200 mm = kk cc = y(ii, kk) IF cc = 0 THEN PRINT "Matrix error - perhaps uou have called for more nodes" PRINT "then you have actually used." PRINT "check the second entry in the fourth line of data" PRINT "input data, the one just after the two comment lines." END IF 1200 NEXT kk ' IF ii = mm THEN 1220 SWAP j(mm + 20), j(ii + 20) ' FOR kk = 1 TO nn SWAP y(kk, ii), y(kk, mm) NEXT kk ' 1220 y(ii, ii) = 1 'detm = detm * cc 'detm = value of the determinant ' 1230 FOR mm = 1 TO nn y(ii, mm) = y(ii, mm) / cc NEXT mm ' FOR mm = 1 TO nn IF ii = mm GOTO 1250 cc = y(mm, ii) IF cc = 0 GOTO 1250 1240 y(mm, ii) = 0 ' FOR kk = 1 TO nn y(mm, kk) = y(mm, kk) - cc * y(ii, kk) NEXT kk ' 1250 NEXT mm ' NEXT ii ' FOR ii = 1 TO nn IF j(ii + 20) = ii GOTO 1280 mm = ii 1260 mm = mm + 1 IF j(mm + 20) = ii GOTO 1270 IF nn > mm GOTO 1260 1270 j(mm + 20) = j(ii + 20) ' FOR kk = 1 TO nn SWAP y(ii, kk), y(mm, kk) NEXT kk ' j(ii + 20) = ii 1280 NEXT ii ' dtnrm = ABS(detm) / pd ' END SUB ' '******************************************** SUB modifyL (nbran, nls, ib(), b(), y(), tag()) ' FOR i = 1 TO nbran IF ib(i, 4) = 8 THEN j1 = ib(i, 2) 'branch number of first l j2 = ib(i, 3) 'branch number of second l j3 = ib(j1, 5) 'tag number of first l j4 = ib(j2, 5) 'tag number of second l y(j3, j4) = b(i, 1) y(j4, j3) = b(i, 1) tag(j3, j4) = i 'keep track of branch involved tag(j4, j3) = i END IF NEXT i ' nls = size for matrix inversion CALL invertH(nls, y()) ' FOR i = 1 TO nls FOR j = 1 TO nls k1 = tag(i, j) IF ABS(y(i, j)) < 1E-20 THEN b(k1, 4) = 1E+20 ELSE b(k1, 4) = 1 / y(i, j) END IF NEXT j NEXT i ' END SUB '**************************************************** SUB nonlinZ (branch%, ib(), b(), e3(), ci3(), cn(), t, t1, K9) ' ------- evaluate non-linear elements ----- i = 1 1010 IF cn(1, i) = ib(branch%, 1) THEN 1020 i = i + 2 IF i > 16 THEN PRINT "Error - a non-linear branch is defined wrong" PRINT "Your entry was "; ib(i, 4); " and i was "; i CALL dataerror ELSE GOTO 1010 END IF 1020 SELECT CASE ib(branch%, 4) ' CASE 9 ' ---- voltage dependent resistor ---- cn2 = (e3(branch%, 1) + e3(branch%, 2)) / 2 'search on average of old and new voltages IF b(branch%, 2) = 1 AND cn2 < 0 THEN j9 = -1 ELSE j9 = 1 END IF IF b(branch%, 2) = 1 THEN cn2 = ABS(cn2) END IF GOSUB 1050 cn8 = (cn(n4, i + 1) - cn(n5, i + 1)) / (cn(n4, i) - cn(n5, i)) 'returns 1/r ci3(branch%, 3) = j9 * (cn(n4, i) * cn8 - cn(n4, i + 1)) 'intercept current b(branch%, 5) = cn8 GOTO 1030 ' CASE 10 ' -- time dependent resistor -- ' control% = b(branch%, 4) IF e3(control%, 1) < e3(branch%, 3) THEN 1030 'vcontrl= cn3 THEN cn2 = t - cn3 ELSE GOTO 1030 END IF GOSUB 1050 cn8 = (cn(n4, i + 1) - cn(n5, i + 1)) / (cn(n4, i) - cn(n5, i)) cnr = (cn2 - cn(n4, i)) * cn8 + cn(n4, i + 1) b(branch%, 5) = 1 / cnr 'cnr=resistance GOTO 1030 ' CASE 11 ' -- energy dependent resistor -- ' energy = e3(branch%, 6) IF energy < b(branch%, 2) THEN 1030 ELSE cn2 = energy GOSUB 1050 cn8 = (cn(n5, i + 1) - cn(n4, i + 1)) / (cn(n5, i) - cn(n4, i)) cnr = (cn2 - cn(n4, i)) * cn8 + cn(n4, i + 1) b(branch%, 5) = 1 / cnr 'cnr=resistance GOTO 1030 ' CASE 12 ' -- flux dependent inductor -- ' cn2 = ABS(e3(branch%, 7)) GOSUB 1050 cn8 = (cn(n4, i + 1) - cn(n5, i + 1)) / (cn(n4, i) - cn(n5, i)) 'returns 1/l b(branch%, 5) = .5 * t1 * cn8 ci3(branch%, 4) = (cn(n4, i) - cn(n4, i + 1) / cn8) * SGN(e3(branch%, 7)) 'offset of flux linkages GOTO 1030 END SELECT ' ' ' ---------check for change of slope----------- ' 1030 IF cn8 = cn(2, i) THEN cn(2, i) = cn8 cn(2, i + 1) = 1 ELSE cn(2, i) = cn8 cn(2, i + 1) = 0 K9 = 1 'slope (and topology) have changed END IF GOTO 1072 '----------------------------------------------- 1050 n4 = cn(1, i + 1) 1052 n5 = n4 + 1 va = cn(n4, i) vb = cn(n5, i) 1054 IF vb = -99 THEN n4 = n4 - 1 ELSEIF cn2 >= vb THEN n4 = n4 + 1 GOTO 1052 ELSEIF cn2 < va THEN n4 = n4 - 1 IF n4 <= 3 THEN n4 = 3 GOTO 1060 END IF GOTO 1052 END IF 1060 IF cn2 >= 0 THEN n5 = n4 + 1 ELSEIF cn2 = va THEN n5 = n4 n4 = n4 - 1 ELSE n5 = n4 + 1 END IF cn(1, i + 1) = n4 1070 RETURN '-------------------------------------------- 1072 END SUB '************************************************ SUB printfour (i, ib(), b(), a$, g$) ' '------- print first four descriptors of branches ---- ' f1$ = "### ### ### \ \ +#.###^^^^ " f2$ = " +#.###^^^^ " PRINT USING f1$; ib(i, 1); ib(i, 2); ib(i, 3); a$; '; 'print elements 1-4 IF ib(i, 4) < 9 THEN PRINT USING f2$; b(i, 1); END SUB '************************************************** SUB readnonlin (i, b(), nnl, cn()) ' '----- read and print non-linear descriptors----- ' PRINT nnl = nnl + 1 IF nnl > 8 THEN PRINT "Error - only 8 non-linear elements allowed" END IF cn(1, 2 * nnl - 1) = i 'load branch number cn(1, 2 * nnl) = b(i, 1) + 2 'load starting segment k = 3 1300 INPUT #1, cn(k, 2 * nnl - 1) IF cn(k, 2 * nnl - 1) = -99 THEN GOTO 1310 INPUT #1, cn(k, 2 * nnl) IF k > 38 THEN PRINT "Error - too many elements in non-linear branch" PRINT "max is 38 and value read was number "; k CALL dataerror END IF nl$ = " +#.###^^^^ +#.###^^^^" PRINT USING nl$; cn(k, 2 * nnl - 1); cn(k, 2 * nnl) k = k + 1 IF g$ <> "g" THEN CALL examine(g$) GOTO 1300 1310 ' END SUB '******************************************** SUB waveshape (t, t3(), ix3, v4, v5, v6, v7, vmag) ' -- subroutines for instantaneous values of sources --------- SELECT CASE ix3 ' CASE 1 'trapezoidal (type 1) v8 = v4 + v6 v9 = v7 + v8 v10 = v9 + v6 IF t < v4 OR t > v10 THEN vmag = 0 ELSEIF t >= v8 AND t <= v9 THEN vmag = v5 ELSEIF t >= v4 AND t < v8 THEN vmag = v5 * (t - v4) / v6 ELSE vmag = v5 * ((v10 - t) / v6) END IF ' ------------------ CASE 2 'triangular (type 2) v8 = v4 + v6 v9 = v8 + v7 IF t < v4 OR t > v9 THEN vmag = 0 ELSEIF t >= v4 AND t < v8 THEN vmag = v5 * (t - v4) / v6 ELSE vmag = v5 * (v9 - t) / v7 END IF ' ------------------- CASE 3 'filtered double exponential (type 3) alpha = v6 beta = v7 gamma = 20 * v7 'filter time constant IF t < v4 THEN vmag = 0 ELSE IF ABS(alpha * (t - v4)) > 80 THEN v8 = 0 ELSE v8 = EXP(-alpha * (t - v4)) END IF IF ABS(beta * (t - v4)) > 80 THEN v9 = 0 ELSE v9 = EXP(-beta * (t - v4)) END IF IF ABS(gamma * (t - v4)) > 80 THEN v10 = 0 ELSE v10 = EXP(-gamma * (t - v4)) END IF fac1 = gamma / (gamma - alpha) fac2 = gamma / (gamma - beta) 'vmag = v5 * (v8 - v9) 'true double exponential vmag = v5 * (fac1 * (v8 - v10) - fac2 * (v9 - v10)) END IF ' ------------------ CASE 4 'modified exponential (type 4) v8 = (t * v7) ^ v4 v8 = v8 / (v8 + 1) v9 = t * v6 IF ABS(v9) > 80 THEN vmag = 0 ELSE vmag = v5 * v8 * EXP(-v9) END IF ' ------------------- CASE 5 'sine-squared (type 5) v8 = v4 + v6 IF t < v4 THEN vmag = 0 ELSEIF t > v8 THEN vmag = v5 * EXP(-((t - v8) * v7)) ELSE vmag = SIN((t - v4) * 1.5706 / v6) vmag = v5 * vmag * vmag END IF ' ------------------- CASE 6 'offset sinusoid (type 6) IF t < v4 THEN vmag = 0 ELSE vmag = v5 * 1.41421 * SIN(6.2832 * v6 * (t - v4) + v7 / 57.296) END IF ' ------------------- CASE 7 'increasing sinusoid (type 7) IF t < v4 THEN vmag = 0 ELSE vmag = v5 * 1.41421 * SIN(6.2832 * v6 * (t - v4)) * (1 - EXP(-(t - v4) * v7)) END IF ' ------------------ CASE 8 'decaying sinusoid (type 8) IF t < v4 THEN vmag = 0 ELSE vmag = v5 * SIN(6.2832 * v6 * (t - v4)) * EXP(-v7 * (t - v4)) END IF ' ------------------- CASE 9 'decaying cosinusoid (type 9) IF t < v4 THEN vmag = 0 ELSE vmag = v5 * COS(6.2832 * v6 * (t - v4)) * EXP(-v7 * (t - v4)) END IF ' ------------------- CASE 10 'tabular (type 10) itab = 1 1600 tt1 = t3(itab, 1) tt2 = t3(itab + 1, 1) IF t - v4 < tt1 OR tt2 = -99 THEN vmag = t3(itab, 2) GOTO 1602 ELSEIF t - v4 < tt2 THEN ttr = tt2 - tt1 tvr = t3(itab + 1, 2) - t3(itab, 2) vmag = ((t - v4) - tt1) * tvr / ttr + t3(itab, 2) GOTO 1602 ELSE itab = itab + 1 END IF IF itab > 20 THEN vmag = t3(itab, 2) GOTO 1602 ELSE GOTO 1600 END IF END SELECT ' 1602 END SUB