%%HP:T(3)A(R)F(.);
DIR
@ --------------------------------------------------------
@ Title   : CSIM (a simple circuit simulator for the HP48)
@ Version : 2.61
@ Author  : Per Stenius
@ LastEdit: 16.12.91
@ Copyright Per Stenius (1991)
@ --------------------------------------------------------

CST
{Csim View node ymin ymax CLRSC
outp CIR\-> CIR \->CIR Setup dc
w ac A\->L t tstep tran
X G C Cc W Wlist
Euler iterdc}

Csim
\<<
  "  Csim_HP-48 2.61


 (c) Per Stenius 1991" CLLCD 2 DISP
  1 WAIT CLLCD
  "Setup?" "Y" INPUT
  IF
  "Y" SAME
  THEN
    "Wait..." CLLCD 1 DISP
    IF
    DEPTH 0 ==
    THEN
      CIR CIR\->
    END
    Setup
  END
  "Analysis? (D, A, T)" "" INPUT
  \-> analysis
  \<<
    CASE

      analysis "D" SAME
      THEN
        dc
      END

      analysis "A" SAME
      THEN
        "Sweep range?" {":wstart:
:wstop:" { 1 0 } V } INPUT
        OBJ\-> \-> wstart wstop
        \<<
          wstop wstart - 130 /
          'wstep' STO
          wstart 'w' STO
          'acplot' STEQ
          wstart wstop XRNG
          ymin ymax YRNG
          'w' INDEP
          DRAX                @ Add ERASE to clear PICT
          {(0,0) "jw" "f(jw)"} AXES LABEL
          DRAW GRAPH
        \>>
      END

      analysis "T" SAME
      THEN
        "Sweep range?" {":tstart:0
:tstep:0
:tstop:1" { 3 0 } V } INPUT
        OBJ\->
        \-> tstart ttstep tstop
        \<<
          IF
          ttstep 0 ==
          THEN
            tstop tstart - 130 /
            'tstep' STO
          ELSE
            ttstep 'tstep' STO
          END
          tstart tstep 130 * XRNG
          ymin ymax YRNG
          't' INDEP
          DRAX                @ Add ERASE to clear PICT
          {(0,0) "t" "f(t)"} AXES LABEL
          IF
          Euler NOT
          THEN
            tstep 2 / 'tstep' STO
            'tranTR' STEQ
          ELSE
            'tranBE' STEQ
          END
          G tstep * C + INV
          'iChG' STO
          DRAW GRAPH
        \>>
      END
    END
  \>>
\>>

outp                          @ Enables user defined
\<< node GET                  @ calculations
\>>

dc
\<<                           @ DC analysis
 Wlist\->W W checkCc
 Gdc / DUP 'X' STO            @ The result vector is
\>>                           @ returned to the stack

checkCc
\<<
 Cc C\->R DUP 
 IF
 CNRM NOT SWAP RNRM NOT AND
 THEN
   G + 'Gdc' STO
 ELSE
   "Im{Cc} \=/ 0 IN DC" DOERR 
 END
\>>

iterdc                        @ Iterative DC analysis, max 100
\<< 0 \-> i                   @ iterations
  \<<
    DO
      X dc
    UNTIL
    ==
    'i' INCR 100 > OR
    END
    IF
    i 100 >
    THEN
      "100 ITERATIONS
CHECK CONVERGENCE" 1 DISP 1 FREEZE
    ELSE
      dc
    END
  \>>
\>>

ac
\<<
  Wlist\->W W G C w * R\->C   @ AC analysis
  Cc + / DUP 'X' STO          @ The result vector is
\>>                           @ returned to the stack

tran
\<<                           @ Trapezoidal approx.
 iChG
 W Wlist\->W W + tstep 2 / *
 C G tstep 2 / * - X * + *
 DUP 'X' STO
 t tstep + 't' STO
\>>

acplot
\<<
  ac outp                     @ outp is always called last
  wstep w + 'w' STO           @ in a plotting program
\>>

tranBE
\<<
 iChG                         @ Inverse Euler approx.
 Wlist\->W W tstep *          @ Returns the next result to stack
 C X * + *                    @ Used as default when plotting
 DUP 'X' STO outp             @ outp is always called last
\>>

tranTR
\<<                           @ Trapezoidal approx.
 iChG
 W Wlist\->W W + tstep *
 C G tstep * - X * + *
 DUP 'X' STO outp             @ outp is always called last
\>>

Wlist\->W                     @ Functional values -> numerical
\<<
 Wlist LIST\-> 1 SWAP
 START
    \->NUM
    dim ROLL
 NEXT
 dim 1 getpos \->ARRY
 'W' STO
\>>

Setup
\<<
  -3 CF                       @ Set symbolic mode
  -17 SF -18 CF               @ and radian mode
  0 't' STO
  0 'ndim' STO
  0 'bdim' STO
  DEPTH 1 SWAP
  START
     1 GETI
     \-> cmptype
     \<<
       IF
       cmptype 'm' SAME NOT  @ Not a component!
       THEN
         cmptype
         incbdim GETI
         incndim GETI
         incndim
         IF
         cmptype 'O' SAME    @ Components with 4 nodes
         cmptype 'M' SAME OR @ New two-ports: add type here!
         cmptype 'T' SAME OR
         cmptype 'g' SAME OR
         cmptype 'r' SAME OR
         cmptype 'a' SAME OR
         cmptype 'u' SAME OR
         cmptype 'y' SAME OR
         cmptype 'z' SAME OR
         cmptype 'A' SAME OR
         THEN
           GETI incndim
           GETI incndim
         END
       END
       DROP DEPTH ROLL
     \>>
  NEXT
  ndim bdim + 'dim' STO
  [[ 0 ]] dim DUP getpos RDM DUP
  'G' STO 'C' STO
  [[ (0,0) ]] dim DUP getpos RDM
  'Cc' STO
  [[ 0 ]] dim 1 getpos RDM
  DUP 'X' STO 'W' STO
  1 dim
  START
     0
  NEXT
  dim \->LIST 'Wlist' STO
  DEPTH 1 SWAP
  START
     IFERR
        DUP 1 GET
        loadmatrix
        DEPTH ROLL
     THEN
        "SYNTAX ERROR" DOERR
     END
  NEXT
  DEPTH \->LIST 'CIR' STO
\>>

loadmatrix
\<< \-> cmptype
  \<<
    DUP 2 GET
    2 PICK 3 GET              @ cmp n1 n2
    CASE
       cmptype 'J' SAME       @ Ideal current source
       THEN
          getval
          putJ
       END

       cmptype 'E' SAME       @ Ideal voltage source
       THEN
          getval
          getbranch
          putE
       END

       cmptype 'G' SAME       @ Conductor and capacitor
       cmptype 'C' SAME OR
       THEN
          getval
          cmptype putGC
       END

       cmptype 'R' SAME 
       cmptype 'L' SAME OR    @ Resistor and inductor
       THEN
          getval
          getbranch
          IF
          cmptype 'R' SAME
          THEN
            'G'
            putRL
          ELSE
            putL
          END
       END

       cmptype 'Z' SAME       @ Constant valued impedance
       THEN
          getval INV
          putY
       END

       cmptype 'Y' SAME       @ Constant valued admittance
       THEN
          getval
          putY
       END

       cmptype 'S' SAME       @ Short-circuit
       THEN
          getval              @ n1 n2 branch          
          putS
       END

       cmptype 'O' SAME       @ Ideal opamp
       THEN
          getn34
          5 PICK 6 GET        @ n1 n2 n3 n4 branch
          putO
       END

       cmptype 'M' SAME       @ Transformer
       THEN
          getn34vb
          7 PICK 8 GET
          8 PICK 9 GET
          9 PICK 10 GET       @ n1 n2 n3 n4 l1 l2 m b1 b2
          putM
       END

       cmptype 'T' SAME       @ Lossless transmission line
       THEN
          getn34vb            @ n1 n2 n3 n4 ll Zo
          putT
       END

       cmptype 'm' SAME       @ Mutual inductance
       THEN
          getval
          putm                @ b1 b2 val
       END

       cmptype 'g' SAME       @ VCCS
       THEN
          getn34
          5 PICK 6 GET        @ n1 n2 n3 n4 val
          putg
       END

       cmptype 'r' SAME       @ CCVS
       THEN
          getn34vb
          7 PICK 8 GET        @ n1 n2 n3 n4 val b1 b2
          putr
       END

       cmptype 'p' SAME       @ CCVS version 2
       THEN
          getvb1b2
          putp                @ n3 n4 val b1 b2
       END

       cmptype 'a' SAME       @ CCCS
       THEN
          getn34vb            @ n1 n2 n3 n4 val branch
          puta
       END

       cmptype 'b' SAME       @ CCVS version 2
       THEN
          getn34
          putb                @ n3 n4 val b
       END

       cmptype 'u' SAME       @ VCVS
       THEN
          getn34vb            @ n1 n2 n3 n4 val branch
          putu
       END

       cmptype 'z' SAME       @ z-parameters (two-port)
       THEN
          getn34v1234         @ n1 n2 n3 n4 y11 y12 y21 y22
          {2 2} \->ARRY INV
          ARRY\-> DROP
          puty
       END

       cmptype 'y' SAME       @ y-parameters (two-port)
       THEN
          getn34v1234         @ n1 n2 n3 n4 y11 y12 y21 y22
          puty
       END

       cmptype 'A' SAME       @ ABCD-parameters (two-port)
       THEN
          getn34v1234         @ n1 n2 n3 n4 A B C D
          ABCDtoy
          puty
       END
                              @ Add new components here!
    END
  \>>
\>>

ABCDtoy
\<< \-> A B C D
  \<<
    D B /
    C D A * B / -
    B INV NEG
    A B /
  \>>
\>>

putGC                         @ Routines to load component
\<< \-> n1 n2 value type      @ stamp into matrix (or vector)
  \<<                         @ Add new stamps here!
    value n2 n1 checknodes
    type RCL
    n1 n2 value puty2
    type STO
  \>>
\>>

putRL
\<< \-> n1 n2 value branch matr
  \<<
    branch n2 n1 checknodes   @ Enables short-circuits
    n1 n2 branch putL2
    matr RCL
    branch DUP value NEG putmatrix
    matr STO
  \>>
\>>

putJ
\<< \-> n1 n2 value
  \<< 
    value n2 n1 checknodes
    Wlist DUP
    IF n1 0 >
    THEN
       n1 GET value - n1 SWAP
       PUT DUP
    END
    IF n2 0 >
    THEN
       n2 GET value + n2 SWAP PUT
    ELSE
       DROP
    END
    'Wlist' STO
  \>>
\>>

putE
\<< \-> n1 n2 value branch
  \<<
    value n2 n1 checknodes
    n1 n2 0 branch putL
    Wlist DUP
    branch GET value +
    branch SWAP PUT
    'Wlist' STO
  \>>
\>>

putM
\<< \-> n1 n2 n3 n4 l1 l2 m b1 b2
  \<<
    n1 n2 l1 b1 putL
    n3 n4 l2 b2 putL
    b1 b2 m putm
  \>>
\>>

putm
\<< \-> b1 b2 m
  \<<
    m b1 b2 checknodes
    C
    b1 b2 m NEG putmatrix
    b2 b1 m NEG putmatrix
    'C' STO
  \>>
\>>

putS
\<< \-> n1 n2 b
  \<<
    n1 n2 0 b putL
  \>>
\>>

putT
\<< \-> n1 n2 n3 n4 ll Zo
  \<<
    ll 2 \135 * * \->NUM \-> gamma
    \<<
      ll n1 n3 checknodes
      IF
      n2 n4 \139
      THEN
        "n2 MUST EQUAL n4 IN T" DOERR
      ELSE
        'INV(i*Zo*SIN(gamma))' \->NUM
        Cc
        n1 n3 4 PICK puty2
        SWAP
        'COS(gamma)-1' \->NUM *
        SWAP
        n1 n2 4 PICK puty2
        n3 n4 4 ROLL puty2
        'Cc' STO
      END
    \>>
  \>>
\>>

putg
\<< \-> n1 n2 n3 n4 value
  \<<
    value n2 n1 checknodes
    value n3 n4 checknodes
    G
    n1 n2 n3 n4 value putg2
    'G' STO
  \>>
\>>

putr
\<< \-> n1 n2 n3 n4 val b1 b2
  \<<
    b1 n2 n1 checknodes
    n1 n2 b1 putS             @ Short circuit
    n3 n4 val b1 b2 putp
  \>>
\>>

putp
\<< \-> n3 n4 val b1 b2
  \<<
    val n3 n4 checknodes
    G
    b2 n3 1 putmatrix
    b2 n4 -1 putmatrix
    b2 b1 val NEG putmatrix
    n3 b2 1 putmatrix
    n4 b2 -1 putmatrix
    'G' STO
  \>>
\>>

putu
\<< \-> n1 n2 n3 n4 value branch
  \<<
    value n2 n1 checknodes
    branch n3 n4 checknodes
    G
    branch n1 value NEG putmatrix
    branch n2 value putmatrix
    branch n3 1 putmatrix
    branch n4 -1 putmatrix
    n3 branch 1 putmatrix
    n4 branch -1 putmatrix
    'G' STO    
  \>>
\>>

puta
\<< \-> n1 n2 n3 n4 val branch
  \<<
    val n2 n1 checknodes
    n1 n2 branch putS         @ Short circuit
    n3 n4 val branch putb
  \>>
\>>

putb
\<< \-> n3 n4 val branch
  \<<
    val n3 n4 checknodes
    G
    n3 branch val putmatrix
    n4 branch val NEG putmatrix    
    'G' STO    
  \>>
\>>

putO
\<< \-> n1 n2 n3 n4 branch
  \<<
    1 n2 n1 checknodes
    1 n3 n4 checknodes
    G
    branch n1 1 putmatrix
    branch n2 -1 putmatrix
    n3 branch 1 putmatrix
    n4 branch -1 putmatrix
    'G' STO
  \>>
\>>

putY
\<< \-> n1 n2 value
  \<<
    value n2 n1 checknodes
    Cc
    n1 n2 value puty2
    'Cc' STO
  \>>
\>>

puty
\<< \-> n1 n2 n3 n4 y11 y12 y21 y22
  \<<
    y11 n2 n1 checknodes
    y22 n2 n1 checknodes
    Cc
    n1 n2 y11 puty2
    n3 n4 y22 puty2
    n1 n2 n3 n4 y21 putg2
    n3 n4 n1 n2 y12 putg2
    'Cc' STO
  \>>
\>>

putL
\<<
 'C' putRL
\>>

putL2
\<< \-> n1 n2 branch
  \<<
    G
    n1 branch 1 putmatrix
    n2 branch -1 putmatrix
    branch n1 1 putmatrix
    branch n2 -1 putmatrix
    'G' STO
  \>>
\>>

putg2
\<< \-> n1 n2 n3 n4 value
  \<<
    n3 n1 value putmatrix
    n4 n2 value putmatrix
    n3 n2 value NEG putmatrix
    n4 n1 value NEG putmatrix
  \>>
\>>

puty2
\<< \-> n1 n2 value
  \<<
    n1 n1 value putmatrix
    n2 n2 value putmatrix
    n1 n2 value NEG putmatrix
    n2 n1 value NEG putmatrix
  \>>
\>>

putmatrix
\<< \-> row col val
  \<<
    IF
    row col AND
    THEN
       row col getpos
       DUP2                   @ matrix in level two
       GET val +
       PUT
    END
  \>>
\>>

incbdim                       @ Increase matrix dimension
\<< \-> cmptype               @ (branch)
  \<<
    IF
    cmptype 'E' SAME
    cmptype 'R' SAME OR
    cmptype 'L' SAME OR
    cmptype 'S' SAME OR
    cmptype 'O' SAME OR
    cmptype 'u' SAME OR
    cmptype 'a' SAME OR
    cmptype 'p' SAME OR
    THEN
       bdim 1 + 'bdim' STO
    ELSE
       IF
       cmptype 'M' SAME
       cmptype 'r' SAME OR
       THEN
          bdim 2 + 'bdim' STO
       END
    END
  \>>
\>>

incndim                       @ Increase matrix dimension
\<< \-> x                     @ (node)
  \<<
    IF
    x ndim >
    THEN
       x 'ndim' STO
    END
  \>>
\>>

checknodes
\<< \-> value n2 n1
  \<<
    CASE
       n1 0 <
       n2 0 < OR
       THEN
          "NEGATIVE NODE NO." DOERR
       END

       n1 n2 ==
       THEN
          "BOTH NODES SAME" DOERR
       END

       value 0 SAME
       THEN
          "ZERO VALUE OR BRANCH" DOERR
       END
    END
  \>>
\>>

getn34
\<<
 3 PICK 4 GET
 4 PICK 5 GET
\>>

getval
\<<
 3 PICK 4 GET
\>>

getvb1b2
\<<
 getn34
 5 PICK 6 GET
\>>

getn34vb
\<<
 getvb1b2
 6 PICK 7 GET
\>>

getn34v1234
\<<
 getn34vb
 7 PICK 8 GET
 8 PICK 9 GET
\>>

getbranch
\<<
 4 PICK 5 GET
\>>

getpos
\<<
 2 \->LIST
\>>

View                          @ Stack-View application
\<<
 PICT RCL \-> pict
  \<<
    PICT PURGE 1
    DEPTH 1 - 10 MIN
    DUP
    IF
      8 >
    THEN # 6d 1
    ELSE # 8d 2
    END \-> rowht tsize
    \<<
      FOR I
        PICT # 0d 65 I rowht * -
        2 \->LIST I \->STR ": "
        + I 3 + PICK \->STR +
        tsize \->GROB GOR
      NEXT
      { } PVIEW pict
      PICT STO
    \>>
  \>>
\>>

CV
\<< \-> node
  \<<
    X node GET
  \>>
\>>

CI
\<< \-> branch
  \<<
    X branch GET
  \>>
\>>

CLRSC                         @ Clears the screen
\<<
 ERASE
\>>

A\->L
\<<
 ARRY\-> DROP 4 \->LIST
\>>

CIR\->
\<<
 LIST\-> DROP
\>>

\->CIR
\<<
 DEPTH \->LIST
\>>

MTXSLV                        @ Solves a matrix equation Ax=B
\<< \-> B A                   @ Increased accuracy (iteration)
  \<<
    B A / B A
    3 PICK RSD A / +
  \>>
\>>

@ -------------------------------------------------------------
@               Default values and a sample circuit

CIR                           @ Sample circuit
{
 {E 1 0 'IFTE(t MOD 2 > 1,-1,1)' 4}
 {E 3 0 'IFTE(CV(2) > 0,-1,1)' 5}
 {G 1 2 10}
 {C 2 0 2}
}

ymin
-1

ymax
1

node
2

w                             @ Angular frequency (omega)
0

t                             @ Time
0

Euler
1

END

