\ FACS  - Forth Analog Computer Simulator  PropForth Version 1.1  NG Lordi 3/2011
\       PropForth version of a program originally written in F83 Forth (1986). 

\        Emulates a 48 Amplifier 24 Coefficient Attentuator Analog Computer

fl

fswrite propfacs.f

: FACS ;

decimal

\ Utilities  ( Words not available in current PropForth version. )

\ No operation.
: noop ;

\ depth ( -- n )  leaves stack depth on top of stack. 
: depth  _sttop _stptr COG@ - 3 -  ;

\ hold ( c -- ) places character (e.g., decimal point) in output number image
: hold  -1 >out W+! pad>out C! ;

\ The following functions are used to determine the time used by words in between.
\ pretime  ( -- ) stores the current clock count in CNT.
variable CNT
: pretime  cnt COG@ CNT L! ;
\ postime (  -- ) prints elapsed time in microseconds betweeen pretime and postime.
\ Value is adjusted ( 13 usec), assuming words are used in compile mode
: postime cnt COG@ CNT L@ - abs 80 / 13 - cr . cr ;   

\ The following words implement Bartholdi's "TO" solution.

wvariable %TO
0 %TO W!
: TO -1 %TO W! ;
: LTO?  %TO W@ if  L!  0 %TO W!  else  L@  then  ;
: WTO?  %TO W@ if  W!  0 %TO W!  else  W@  then  ;
: CTO?  %TO W@ if  C!  0 %TO W!  else  C@  then  ;

\ These words implement a sign check.

wvariable ?sign
0 ?sign W!
\ ( -x -- x )
: cksign  dup 0< if -1 ?sign W! else 0 ?sign W! then  abs  ; 

\ This word converts an odd integer to an even integer.
: oddtoeven  2 /mod + 2* ;

\ For convience.
: ON  -1  ;
: OFF  0  ;

\ Deferred Words

: model noop ;
: output noop ;
: INPUT noop ;
: smat noop ;
: slopes noop ;
: INTGR noop ;
wvariable modaddr    \ store deferred simulation model pfa
' model modaddr W!
wvariable outaddr    \ store deferred output mode pfa
' output outaddr W!
wvariable inaddr      \ store deferred input mode pfa
' INPUT inaddr W! 
wvariable mataddr    \ store deferred mat-vector pfa
' smat mataddr W! 
wvariable slopaddr   \ store deferred slopes-table pfa
' slopes slopaddr W!
wvariable intaddr    \ store deferred integrator pfa
' INTGR intaddr W!

\ Data Types

\ vector ( n -- ) <name>  defines <name> as a one dimensional array of n longs.
: vector lockdict create $C_a_dovarl w, 0 l, 4* allot 150 w, forthentry freedict ;

52 vector y()    \ state variable array - allows for 48 blocks 
                 \ (4 elements are reserved for system use)
12 vector ic()   \ initial Condition array - 2 elements used per each ic
24 vector k()    \ parameter array 

\ Actual words used to access IC(), Kx, and Y() vectors. 
: IC() ic() swap 4* + LTO?  ;
: Y() y() swap 4* + LTO? ;
: Kx  k() swap 4* + LTO? ;

\ matrix ( n1 n2 -- ) <name> defines<name> as a two-dimensional array of n1*n2 longs.
\ n1 is the no. rows and n2 is the no. columns.
: matrix lockdict create $C_a_dovarl w, dup w, swap dup w, * 4* 2+ allot 0 w, 150 w, 
      forthentry freedict ;

\ Used to access defered matrix smat.
: SMAT  smat rot 1- 4* rot + 4* + LTO? ;

\ array ( n -- ) <name>  defines <name> as a one dimensional array of n words.
: array  lockdict create $C_a_dovarw w, 0 w, 2* allot 150 w, forthentry freedict ;

10 array n()  \ System array which stores node numbers of state variables for
              \ which output is requested: default 10. 
20 array b()  \ System array which stores current values of all logic nodes.
 
\ Actual words used to access arays N() & b().
: N() n() swap 2* + WTO? ;
: B() b() swap 2* + WTO? ;

\ Scalars - System Parameters (Scaling factor is 10000 , ie, 1 = 10000.)

wvariable dt         \ The integration step interval. 
wvariable (dt)       \ Variable used to temporarily store dt.
wvariable comint     \ Communication interval, at least 2*dt. 

wvariable n          \ Total number of defined nodes.
wvariable n1         \ Counts communication intervals 
wvariable l0         \ Number of nodes+1 for which output is requested.

\ Actual words used to access system parameters.
: DT  dt WTO? ;
: (DT)  (dt) WTO? ;
: COMINT comint WTO? ;
: N n WTO? ;
: N1 n1 WTO? ;
: L0 l0 WTO? ;
\ default values
10 TO DT
0 TO (DT)
100 TO COMINT
10 TO N  
0 TO N1          
0 TO L0

\ Other system variables and constants.

variable temp       \ temporary storage of state variable
: TEMP temp LTO? ;  \ Word used to access temp.

decimal
3 wconstant FIX         \ Default no. decimal places alloted.
10000 wconstant scale   \ Parameter scale - 4 decimal precision.
65535 constant lscale   \ Table (log, antilog, sin) scale.
  
301030 constant Lg2     \ Scaled Log10[2] 6 dec. places.
693147 constant Ln2     \ Scaled Ln[2] 6 dec. places.  
31416 constant Pi       \ Scaled 4 dec. places.
27183 constant Ei       \ Scaled 4 dec. places.

decimal
1000000 constant Y0  \ State variable scaling factor (1= 1000000).
10000000 constant Ym  \ Maximum scaled state variable value allowed. 

\ Reference values for state variables, equivalent to 1 scaled units.
1 wconstant +REF         \ Assigned to node 1.
2 wconstant 0REF         \ Assigned to node 2.
3 wconstant -REF         \ Assigned to node 3. 

\ System Messages

: MESS1  cr cr ." SIMULATION RUN COMPLETED " cr ;
: MESS2  cr cr ." EXECUTION TERMINATED " cr reset ;
: +ERR   dup abs Ym > if ." OVER/UNDERFLOW " cr  reset then ;
: 0ERR   dup 0= if ." ZERO NOT ALLOWED " cr  reset then  ;
: -ERR   dup 0< if ." NEGATIVE NOT ALLOWED " cr reset then ;

\ Output

\ ( x -- 'x' )Basic output routine displays vector elements in decimal format, i.e.,
\  'x.xxx'.
: dout  cksign 1000 / <# FIX 0 do # loop 46 hold #s  ?sign W@ if 45 hold then  #> 
        .cstr 4 spaces ; 

\ ( -- ) Outputs time in decimal format (arbitrary units).
: tout  0 Y() 10 / <# FIX 0 do # loop 46 hold #s #> .cstr 2 spaces ; 

\ ( n -- ) Sets no. of nodes to n.
: NODES  TO N ; 

\ ( -- ) Clears array N() (sets to 0 ). 
: ZERO-N()  10 0 do 0 TO i N() loop ; 

\ ( n1 n2 ... -- ) Enter specific nodes for output.
: OUT-NODES  ZERO-N() depth 1+ dup TO L0 1  do TO i N() loop ; 

\ ( -- ) Tabulates output in columns:  time ...  output nodes .
: (display)  tout L0 1 do i N() Y() dout loop cr ; 

\ Store data in ram.
wvariable row
: ROW row WTO? ;
wvariable col
: COL col WTO? ;
: (mstore) L0 1 do i N() Y()  N1 i TO SMAT loop ; 
\ ( n1 n2 -- )  Reads n1 n2 SMAT and displays data previously stored in ram.
: MREAD TO COL 1+ 1 do  COL 1+ 1 do  j i SMAT . space loop cr loop ;

\ XY-Plotter routines for TV system "GreyTV". 
\ Variables setting plotter ranges.   
variable _xscale
variable _yscale
: XSCALE _xscale LTO? ;
: YSCALE _yscale LTO? ;
wvariable _wait
: WAIT _wait WTO? ;
variable _xlevel
variable _ylevel
: XLEVEL _xlevel LTO? ;
: YLEVEL _ylevel LTO? ;
100 TO WAIT
0 TO XLEVEL
0 TO YLEVEL
\ defaults are maximum values for 130X192 tv screen.
120 TO XSCALE
100 TO YSCALE
\ Sends y-value to plotter.
\ ( n --- ) Basic plot routine where Y(n) is the vector to be plotted.
: yplot  Y() 5000 + YSCALE Y0 */ YLEVEL + WAIT delms .byte ;
: xplot  Y() 5000 + XSCALE Y0 */ XLEVEL + WAIT delms .byte ;
\ Enables  vector-time plots where 1 N() contains vector node.
: (typlot) 1 N() yplot  ;
\ Plots two vectors identified by 1 N() and 2 N().
: (typlot2) 1 N() yplot  2 N() yplot ;  
\ Sends x- and y-values to plotter where x is 1 N()Y() and and y is 2 N() Y(). 
: (xyplot) 1 N() xplot 2 N() yplot  ;

\ System Control

\ ( -- ) Alias of ':' for clarity in defining models.
: MODEL:  :  ;  

\ ( -- ) OUTPUT-IS <output> activates <output> as the currect output mode.
: OUTPUT-IS ' outaddr W@ W! ; 

\ ( -- ) INPUT-IS  <input> activates <input> as the current input mode.
: INPUT-IS  ' inaddr W@ W! ; 

\ ( -- ) SIMULATE <model> activates <model> as the current simulation objective.
: SIMULATE ' modaddr W@ W! ;  

\ SMAT-IS <matrix name> activates deferred matrix SMAT.
: SMAT-IS ' mataddr W@ W! ;  

\ (  -- ) Increments time DT units at each execution.
: TIME   0 Y() DT + TO 0 Y() ;  

\ ( -- ) Sets array B() to 0.
: ZERO-B() 20 0 do 0 TO i B()  loop ; 

\ ( -- ) Sets all Y() vector values to 0.
 : ZERO-Y()  N 1+ 0 do 0 TO i Y()  loop  ; 

\ ( -- ) Sets all initial conditions to 0.
: ZERO-IC   12 0 do 0 TO i IC() loop ; 

\ (... y3 n3 y2 y1 n1 -- ) Assigns initial conditions, e.g., 
\ IC(1) = n1, IC(2) = y1, IC(3) = n2, etc. 
: ASSIGN-IC ZERO-IC depth dup TO 0 IC() 1+ 1 do TO i IC() loop ; 

\ ( -- )  Assigns reference values to nodes 1, 2, & 3.
: ASSIGN-REF  Y0 TO 1 Y()  0 TO 2 Y()  Y0 negate TO 3 Y()  ;  

\ ( -- ) Resets all initial conditions and resets time to 0.
: RESET   ZERO-B() ZERO-Y() ASSIGN-REF  0 IC() 1 do i 1+ IC() i IC() TO Y() 2 +loop 
         0 TO N1 ;

\ ( -- ) Initializes model block outputs to starting values.
: INITIALIZE  DT TO (DT) 0 TO DT  model N1 1+ (DT) TO DT output ;  

\ ( n1 -- ) Debugger to static check correctness of block connections. 
\ Prints outputs from node 4 to node n1. 
\ Variable used by RUN to temporarily store no. runs.
wvariable _runo
: RUNO  _runo WTO? ; 

: CHECK TO RUNO RESET  INITIALIZE RUNO 1+ 4 do  i . space i Y() . cr loop ; 

\ Simulation run will pause when any key is pressed , then continue when any key
\ is pressed.
: HALT  key? if key drop 10 delms begin key? until key drop 10 delms then ;
   
\ (n -- ) Simulation will continue for n communication intervals.
: CONTINUE  0 do HALT  COMINT 0 do TIME model  loop N1 1+ TO N1 output  loop   ; 

\ Variable used by RUN to temporarily store no. runs.
wvariable _runo
: RUNO  _runo WTO? ; 
\ ( n -- ) Simulation will run for n communication intervals.
: RUN  TO RUNO  RESET  INITIALIZE  RUNO CONTINUE  MESS1  ; 

\ ( n -- ) Dynamic chack of block connections. Requires EULER integrator.
: CHECK-RUN TO RUNO RESET INITIALIZE RUNO 0 do  COMINT 0 do TIME model  
      DT TO (DT) 0 TO DT  model (DT) TO DT loop output loop ; 

\ PRUN is a version of RUN used for the xy-plotter configuration.
: PRUN TO RUNO 10 delsec RESET DT TO (DT) model (DT) TO DT RUNO 0 do  COMINT 0 do TIME model 
     loop  output  loop  ;

\ Block Definitions

\ Input Functions 

\ ( n1 -- ) Y(n1) = random number 0 - Y0.
: RND  rnd Y0 255 */ TO swap Y() ; 

\ Auxiliary Words
variable _period
: PERIOD _period LTO? ;
variable _xt
: XT _xt LTO? ;
0 TO XT
variable _xo
: XO _xo LTO? ;
\ XE may be + plus or minus.  Sets  zero-level for TRIWAVE and SQWAVE.
variable _xe
: XE _xe LTO? ;
variable _slope
: SLOPE _slope LTO? ;
\ Note that +SLOPE & -SLOPE must be positive values.
variable _+slope
: +SLOPE _+slope LTO? ;
variable _-slope
: -SLOPE _-slope LTO? ;

wvariable _flg
: FLG _flg WTO? ;

\ Ramp
\ ( n1 -- ) Y(n1) is linear from XO to Y0, if > Y0 is 0. 
: RAMP >r 0 Y() SLOPE scale */ XO + dup Y0 >= if 0 TO r> Y() drop 
       else TO r> Y() then ;
 
\ ( n1 -- )  Y(n1) is a step function with stepsize XO and width PERIOD.
: STEP  0 Y() PERIOD /mod drop 0= if dup Y() XO + TO swap Y() else drop then ;

\ Square Wave

\ ( n1 -- ) Y(n1) is a repetitive square wave with amplitude XO.
: IMPULSE  >r 0 Y() PERIOD /mod drop 0= if XO TO r> Y() else 0 TO r> Y() then ;

: SQWAVE >r XT PERIOD /mod swap drop 0 = if XO else 0 then XE + TO r> Y()
     XT PERIOD 2* DT - = if 0 TO XT else XT DT + TO XT then  ;

\ ( n1 -- ) Triangle Wave allows diferent + & - slopes. 
\ Y(n1) is a repetitive triangle wave with a period dependent on the SLOPE value.
: TRIWAVE  dup >r Y() 0 <= if DT TO XT  1 TO FLG else r@ Y() Y0 >= if DT TO XT 2 TO FLG 
    thens FLG 1 = if XT +SLOPE scale */ XE + TO r> Y() else Y0 XT -SLOPE scale */ - 
    XE + TO r> Y() then XT DT + TO XT ;  

\ ( n1 --- ) n1 INPUT  sets the input node for selected input mode, all the above.

\ Function Generator simulates the classic diode function generator.
\ It requires a table of slopes. Break points are set at equal intervals determined
\ by FPERIOD. Y(n1) is output and Y(n2) input. If n2 is 0, the  generator will 
\ produce a time-dependent wave form, otherwise, a transformation of Y(n2). 
\ The definition of a table with 11 elements follows. The last element is a
\ terminating character. Note that slopes are scaled using 10000. slopes is a
\ deferred word, so the user must define this table as part of a model definition. 
\ variable <name> 1000 l, 5000 l, 2000 l, 500 l, -500 l, -700 l, -800 l,...
\ Slope values are calculated assuming DT = 1.  FGAIN is calculated as 10000/DT and
\ adjusted up or down to adjust the output magnitude if necessary. 
: SLOPES 4* slopes + LTO? ;  \ ( n -- )
: SLOPES-IS ' slopaddr W@ W! ; \ Slopes-IS <name>
\ Number of break points.
wvariable _nbpt
: NBPT  _nbpt WTO? ;

variable _fgain
: FGAIN _fgain LTO? ;    
variable _fperiod
: FPERIOD _fperiod LTO? ;

\ ( n1 n2 -- ) Y(n1) = Y(n2)+ Y(n2)*[ n slopes] where n is Y(n1)/FPERIOD.
: FUNGEN Y() FPERIOD 1+ / 1+ dup NBPT > if drop 0 else SLOPES then 
      DT *  FGAIN scale */ over Y() +  TO swap Y() ;

\ Integrators

\ ( n1 n2  -- ) Euler Integrator - Y(n1) = Integral[Y(n2)]
: EULER   Y() DT scale */ over Y() + TO swap Y() ;  

\ ( n1 n2  -- ) 2nd order Rung Kutta Integrator  Y(n1) = Integral(Y(n2). 
: RK2  Y() DT scale */ dup 2000000 + Y0 */ 2/ over Y() + TO swap Y() ;  

\ INTGR-IS <name>  Set <name> as the active integrator, where <name> is EULER or RK2. 
\ Use n1 n2 INTGR in model definitions instead of specifing a specific integrator.   
: INTGR-IS ' intaddr W@ W!  ;

{
\ ( n1 n2 n3 .... -- )  Y(n1) = Integral[Y(n2)+Y(n3)+ ... ]
\ RK2SUM is a summing integrator that combines the functions of SUM and RK2 in one
\ block.  It cannot be used with INTGR.
: RK2SUM  depth 2- swap Y() swap 0 do swap Y() + loop  DT scale */ dup 2000000 +
    Y0 */ over Y() + TO swap Y()  ;  
}
\ Mathematical Operations

\ ( n1 n2 -- )Y(n1) = -Y(n2)
: INV  Y() negate TO swap Y() ; 

\ ( k n1 n2 -- )Y(n1) = k*Y(n2) where k = 0-1  Note: if k is - combines POT with INV.
\ If k > 1  max 10, block is treated as a pot with a Gain of 10.
: POT  Y() rot scale */ TO swap Y() ; 

\ ( n1 n2 n3... -- ) Y(n1) = Y(n2)+Y(n3)+ ....
: SUM  depth 2- swap Y() swap 0 do swap Y() + loop TO swap Y() ; 

\ ( n1 n2 -- )Y(n1) = absolute value of Y(n2)
: ABS  Y() abs TO swap Y()  ; 

\ ( k n1 n2 -- ) Y(n1) = k+Y(n2) where k is a long.
: OFFSET  Y() rot + +ERR TO swap Y() ; 

\ ( n1 n2 n3 -- ) Y(n1) = Y(n2)*Y(n3)
: MULT  Y() swap Y() Y0 */ +ERR TO swap Y()  ; 

\ ( n1 n2 -- )Y(n1) = Y(n2)*Y(n2)
: SQR  Y() dup Y0 */ TO swap Y()  ; 

\ Variables required by square root (Newton's) method.
variable xs
: Xs xs LTO? ;
variable xn
: Xn xn LTO? ;
variable xi
: Xi xi LTO? ;
\ ( n1 n2 -- ) Y(n1) = Square Root of Y(n2) (scaled)
: SQRT  Y() TO Xs Xs TO Xi begin Xi TO Xn  Xn Xs Y0 Xn */ + 2/ dup TO Xi 
        Xn - 0= until Xn TO swap Y() ;

\ ( n1 n2 -- ) Y(n1) = Y(n3)/ Y(n2) - Zero-Error if Y(n2) <= 0.
: DIV  Y() 0REF rot Y() */  TO swap Y() ;

\ (n1 n2 -- ) Y(n1) = Sin[Y(n2)]  Uses prop sine table. Y() assumed to be in radians 
\ limited to angles 0 - 90 degrees.
: SIN  Y() 2000 766616  */ oddtoeven 57344 + W@ Y0 lscale */ TO swap Y() ; 

\ ( n1 n2 -- ) Y(n1) = LOG(base 2)[1 + Y(n2)] or LOG2[Y(n2]+1 in range Y(n1) = 0.5-1.
: LOG2 Y() scale 2442598 */ oddtoeven 49152 + W@ Y0 lscale */ TO swap Y() ; 

\ (n1 n2 -- ) Y(n1) = Antilog2[1 + Y(n2)] - 1  
: ALOG2  Y() scale 2442598 */ oddtoeven 53248 + W@ Y0 lscale */ TO swap Y() ; 

\ Logic and Logic-to-Analog Interfaces 

\ ( n1 n2 b1 -- )Track/Hold block precedes an integrator. If B(b1) is true, 
\     Y(n1)  = 0 else  Y(n1) = Y(n2).
: T/H  if 0 TO swap Y() else Y() TO swap Y() then  ;  

\ ( n1 n2 n3 -- )  B(n3)= -1 if Y(n2) > Y(n3) else B(n3) = 0
: CMP   Y() swap Y() > TO swap B()  ;  

\ ( n1 n2 n3 n4 -- )If B(n4) = -1, Y(n1) = Y(n3) else Y(n1) = Y(n2)
: SWITCH  B() if  Y() TO rot Y() drop  else drop Y() TO swap Y() then  ;

\ ( b1 b2 b3 -- )B(b1) = B(b2) or B(b3)
: OR B() swap B() or TO swap B() ;  

\ ( b1 b2 b3 -- ) B(b1) = B(b2) and B(b3)
: AND B() swap B() and TO swap B() ; 

\ ( b1 b2 -- ) B(b2) = invert B(b1)
: NOT  B() invert TO swap B() ;  

\ ( --- b ) B(b) is random bits on/off (true/false).
: RNDB   rndtf TO swap B() ;
      
\ A Pulse Generator which produces on/off bits at user selected widths and intervals.
\ Generally used with switches to control aspects of specific simulations. When combined
\ with logic elements it emulates a hybrid analog computer.
 
variable _start
wvariable _width
wvariable _interval

: START _start LTO? ;
: WIDTH _width WTO? ;
: INTERVAL _interval WTO? ; 

: PULSE  TO 0 B() START TO TEMP  0 Y() TEMP = if  -1 0 B() TO B() then  
       0 Y() TEMP WIDTH + = if 0 0 B() TO B() TEMP INTERVAL + TO TEMP then ;

\ Nonlinear blocks

\ ( n1 n2 -- ) Y(n1) = only positive values of Y(n2)
: -CLIP  Y() dup 0< if drop 0 then TO swap Y()  ;  

\ ( n1 n2 -- ) Y(n1) = only negative values of Y(n2) 
: +CLIP  Y() dup 0> if drop 0 then TO swap Y()  ;  

\ ( n1 n2 -- ) If Y(n2) is positive Y(n1) = Y0  else if Y(n2) is negative Y(n1) = -Y0.
: BANG-BANG  Y() 0> if +REF Y() else  -REF Y()  then TO swap Y()  ; 

\ ( k1 k2 n1 n2 -- ) Dead Zone block sets Y(n1) = k1 if Y(n2) is >k1 & <k2.
: DEAD  Y() TO TEMP swap >r rot dup TEMP > swap TEMP < and if  
      else drop TEMP then TO r> Y() ;  

\ ( k1 k2 n1 n2 -- ) Limit block sets Y(n1) = k1 if Y(n2) < k1 or k2 if Y(n2) > k2 
\ otherwise, Y(n1) = Y(n2).
: LIMIT  Y() swap >r scale / dup < if drop else swap drop dup > 
     if swap thens scale * TO r> Y() drop ;  

\ ( n1 n2 -- ) Terminate run if Y(n2) is greater than Y(n1).
: STOP  Y() swap Y() > if MESS2 then ; 

\ The Delay Line block is emulated as a FIFO shift register with each cell delayed
\ DT time units.  The shift register is defined as a vector with n cells (longs).
\ The delay time is n*DT.  

: dly() noop ;  \ Deferred vector.
wvariable dlyaddr
' dly() dlyaddr W!
: DLY()  dly() swap 4* +  LTO? ;
: DELAY-IS ' dlyaddr W@ W! ; 
wvariable _dn
: DN _dn WTO? ;
: CLEAR-DLY  DN 1+ 1 do 0 TO i DLY() loop  ;

\ ( n1 n2  -- ) Y(n1) = Y(n2) delayed  DN*DT time units.
: DELAY DN DLY() TO rot Y()  DN 1 do  DN i - DLY()  DN i 1- - TO DLY()
   loop Y() TO 1 DLY()  ; 

...