path: root/benchmarks/dc68.dat

BEGIN NEW DATA CASE
C     BENCHMARK DC-68                                                           
C         Illustration of MODELS using six of Laurent Dube's simple test cases           
C         1st of 16 data subcases.  This was taken from Laurent's file  TC1.DAT.  
C         July of 1995,  this is converted to the new  MODELS STAND ALONE  by
C         adding  STAND ALONE  to the  MODELS  declaration.  Also,  all data
C         between  ENDMODELS  & batch-mode plot cards (just BLANK) was deleted.
     0.1    1.00
       1       1       1       1       1      
MODELS STAND ALONE
STORAGE integer pages:0 real pages:0 ENDSTORAGE
MODEL tc1
 VAR a
 INPUT i
 EXEC
  a:=t
 ENDEXEC
ENDMODEL
USE tc1 AS tc1a
 INPUT i:=0
ENDUSE
RECORD tc1a.a AS a
ENDMODELS
C   Step      Time      MODELS
C                       A
C      0       0.0         0.0
C      1       0.1         0.1
C      2       0.2         0.2
C      3       0.3         0.3
C      4       0.4         0.4
C      5       0.5         0.5
C      6       0.6         0.6
C      7       0.7         0.7
C      8       0.8         0.8
C      9       0.9         0.9
C     10       1.0         1.0
C  Variable maxima :       1.0
C  Times of maxima :       1.0
C  Variable minima :       0.0
C  Times of minima :       0.0
BLANK card terminating plotting
BEGIN NEW DATA CASE
C         Illustration of MODELS using six of Laurent Dube's simple test cases           
C         2nd of 16 data subcases.  This was taken from Laurent's file  TC2.DAT.  
   0.001    0.01
       1       1       1       1       1      
MODELS
INPUT elec90 {v(NOD1)}
      elec91 {i(NOD2)}
      elec93 {switch(NOD2)}
      elec94 {imssv(NOD1)}
      elec95 {imssi(NOD2)}
MODEL TC2
 INPUT i[90..95]
 VAR a[90..95]
 EXEC a[90..95]:=i[90..95] ENDEXEC
ENDMODEL
USE tc2 AS tc2
 INPUT i[90]:=elec90
       i[91]:=elec91
       i[92]:=0
       i[93]:=elec93
       i[94]:=elec94
       i[95]:=elec95
ENDUSE
RECORD tc2.a[90] AS a90
       tc2.a[91] AS a91
       tc2.a[92] AS a92
       tc2.a[93] AS a93
       tc2.a[94] AS a94
       tc2.a[95] AS a95
ENDMODELS
  NOD1  NOD2               1.00
BLANK card after last electric network branch
  NOD2                                                MEASURING
BLANK card ends switches
14NOD1        100.0       50.0      -30.                         -1.0
BLANK card after last electric network source
  NOD1
C    Step      Time      NOD1       MODELS     MODELS     MODELS     MODELS     MODELS     MODELS
C                                   A90        A91        A92        A93        A94        A95
C ***  Phasor I(0) =  8.6602540E+01    Switch  "NOD2  "  to  "      "  closed in the steady-state.
C       0       0.0  86.6025404 86.6025404 86.6025404        0.0        1.0       -50.       -50.
C       1     .1E-2  97.8147601 97.8147601 97.8147601        0.0        1.0        0.0        0.0
C       2      .002  99.4521895 99.4521895 99.4521895        0.0        1.0        0.0        0.0
C       3      .003  91.3545458 91.3545458 91.3545458        0.0        1.0        0.0        0.0
BLANK card ends selective node voltage outputs
C       8      .008  -40.673664 -40.673664 -40.673664        0.0        1.0        0.0        0.0
C       9      .009  -66.913061 -66.913061 -66.913061        0.0        1.0        0.0        0.0
C      10       .01   -86.60254  -86.60254  -86.60254        0.0        1.0        0.0        0.0
C Variable maxima :  99.4521895 99.4521895 99.4521895        0.0        1.0        0.0        0.0
C Times of maxima :        .002       .002       .002        0.0        0.0      .1E-2      .1E-2
C Variable minima :   -86.60254  -86.60254  -86.60254        0.0        1.0       -50.       -50.
C Times of minima :         .01        .01        .01        0.0        0.0        0.0        0.0
BLANK card terminating plotting
BEGIN NEW DATA CASE
PRINTED NUMBER WIDTH, 12, 2, { Keep dT loop precision the same, but 2 blank separators
C         Illustration of MODELS using six of Laurent Dube's simple test cases           
C         3rd of 16 data subcases.  This was taken from Laurent's file  TC9.DAT.  
   0.001    0.01
       1       1       1       1       0      
MODELS
OUTPUT tacs1, tacs2, arcsh, arcch, arcth
MODEL tc9
 INPUT i {dflt:0}
 VAR tacs[1..2]
 OUTPUT tacs[1..2]
 EXEC
  tacs[1]:=cos(2*pi*1000/6*t)
  tacs[2]:=sin(2*pi*1000/6*(t+startstep))
 ENDEXEC
ENDMODEL
USE tc9 AS tc9
 OUTPUT tacs1:=tacs[1]
        tacs2:=tacs[2]
ENDUSE
C  21 September 2002,  add HYPER to demonstrate corrected hyperbolic functions.
C  Previously, the inverse sinh, cosh, and tanh were missing.  Dube did not halt
C  execution, however.  Instead, he returned impossible function valu 88888.88..
C  as first reported by Orlando Hevia of Santa Fe, Argentina.  WSM corrects XPR2
C  Note the internal names such as  arcsinhyperbolic  are impractically long.
C  This is done on purpose to demonstrate Dube's flexibility in this regard.
MODEL HYPER  --  module that serves only to verify 3 inverse hyperbolic function
 VAR    arcsinhyperbolic, arccoshyperbolic, arctanhyperbolic, arg
 OUTPUT arcsinhyperbolic, arccoshyperbolic, arctanhyperbolic
EXEC
  arg := 200 * T  --  argument of 3 function usages to follow
  arcsinhyperbolic := ASINH( arg )  --  test hyperbolic arc sine over [0, 2]
    IF  arg >= 1.0    -- if arg of arc cosh is possible & arc tanh is impossible
    THEN
  arccoshyperbolic := ACOSH( arg )  --  test hyperbolic arc cosine over ...
  arctanhyperbolic := -9999.  
    ELSE  -- Alternate, if arg of arc cosh is impossible & arc tanh is possible:
  arccoshyperbolic := -9999.  
  arctanhyperbolic := ATANH( arg )  --  test hyperbolic arc tangent over ...
    ENDIF
ENDEXEC
ENDMODEL
USE HYPER AS HYPER
 OUTPUT ARCSH:=arcsinhyperbolic  --  ATP output name ARCSH is connected to internal model name
 OUTPUT ARCCH:=arccoshyperbolic  --  ATP output name ARCCH is connected  ...
 OUTPUT ARCTH:=arctanhyperbolic  --  ATP output name ARCTH is connected  ...
ENDUSE
RECORD tacs1 AS tacs1
       tacs2 AS tacs2
       arcsh AS arcsh
       arcch AS arcch
       arcth AS arcth
ENDMODELS
  NOD1  NOD2               1.00
  TACS2 NOD2               1.00
  NOD2  NOD3               1.00
BLANK card after last electric network branch
  NOD2                                                MEASURING
13NOD3                                                                TACS1   1
BLANK card ends switches
17TACS2
14NOD1        100.0       50.0                                   -1.0
60TACS2
BLANK card after last electric network source
  TACS2
C  First  1     output variables are electric-network voltage differences (upper voltage minus lower voltage);
C  Next   5     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C   Step      Time       TACS2       MODELS      MODELS      MODELS      MODELS      MODELS
C                                    TACS1       TACS2       ARCSH       ARCCH       ARCTH
C ***      Phasor I(0) =  1.0000000E+02                     Switch  "NOD2  "  to  "      "  closed in the steady-state.
C      0       0.0          0.0         1.0  .866025404         0.0      -9999.         0.0
C  Switch  "NOD3  "  to  "      "  closing after  1.00000000E-03  sec.
C      1     .1E-2   .866025404         0.5  .866025404   .19869011      -9999.  .202732554
C                                                Switch  "NOD3  "  to  "      "  opening after  2.00000000E-03  sec.
C      2      .002   .866025404         -.5  .56655E-15   .39003532      -9999.   .42364893
C      3      .003   .56655E-15         -1.   -.8660254  .568824899      -9999.  .693147181
BLANK card ends selective node voltage outputs
C      4      .004    -.8660254         -.5   -.8660254  .732668256      -9999.  1.09861229
C  Switch  "NOD3  "  to  "      "  closing after  5.00000000E-03  sec.
C      5      .005    -.8660254         0.5  -.1133E-14  .881373587         0.0      -9999.
C      6      .006   -.1133E-14         1.0  .866025404  1.01597313  .622362504      -9999.
C      7      .007   .866025404         0.5  .866025404  1.13798205  .867014726      -9999.
C                                                Switch  "NOD3  "  to  "      "  opening after  8.00000000E-03  sec.
C      8      .008   .866025404         -.5  .36738E-15  1.24898333  1.04696792      -9999.
C      9      .009   .36738E-15         -1.   -.8660254  1.35044074  1.19291073      -9999.
C     10       .01    -.8660254         -.5   -.8660254  1.44363548   1.3169579      -9999.
BLANK card terminating plotting
BEGIN NEW DATA CASE
PRINTED NUMBER WIDTH, 11, 1,        { Cancel preceding declaration (extra blank)
C         Illustration of MODELS using six of Laurent Dube's simple test cases
C         4th of 16 data subcases.  This was taken from Laurent's file  TC10.DAT  
     1.0     1.0
       1       1       1       1       1      
MODELS
MODEL test
COMMENT
 Provides a quick check on the use of nested pointlist functions
 for the representation of parametric curves.
 Here the parameter is the variable 'r',
 and the main functin is used as y=f(x,r)
 using a=mainfun(x) at various values of r.
ENDCOMMENT
 FUNCTION fun1 POINTLIST (-inf,inf), (0,0), (inf,-inf)
 FUNCTION fun2 POINTLIST (-inf,2*inf), (0,0), (inf,-2*inf)
 FUNCTION fun3 POINTLIST (-inf,3*inf), (0,0), (inf,-3*inf)
 FUNCTION fun4 POINTLIST (-inf,4*inf), (0,0), (inf,-4*inf)
 FUNCTION mainfun POINTLIST 
  (0,0),(1,fun1(r)),(2,fun2(r)),(3,fun3(r)),(4,fun4(r))
 FUNCTION ref(x,r):=-r*x  --Same function in its analytic form
 INPUT i
 VAR a,r
 INIT
  r:=-1
  WHILE r<5 DO
   r:=r+1
   write(' ')
   write('r:',r)
   FOR k:=0 to 3.5 by 0.5 DO
    a:=mainfun(k)
    write('  ',k,': ',a, '   ', ref(r,k))
   ENDFOR
  ENDWHILE
 ENDINIT
 EXEC
 ENDEXEC
ENDMODEL
USE test AS test
 INPUT i:=0
ENDUSE
ENDMODELS
  NOD1                     1.00
BLANK card after last electric network branch
BLANK card ends switches
14NOD1        100.0       50.0                                   -1.0
BLANK card after last electric network source
  NOD1 
BLANK card ends selective node voltage outputs
C       Step      Time       NOD1                                                                                                       
C          0       0.0        100.                                                                                                      
C          1       1.0        100.                                                                                                      
C    Variable maxima :        100.                                                                                                   
C    Times of maxima :         0.0                                                                                                   
C    Variable minima :        100.                                                                                                   
C    Times of minima :         0.0                                                                                                   
BLANK card terminating plotting
BEGIN NEW DATA CASE
C         Illustration of MODELS using six of Laurent Dube's simple test cases           
C         5th of 16 data subcases.  This was taken from Laurent's file  TC11.DAT  
     0.1    1.00
       1       1       1       1       1      
MODELS
MODEL tc11
 VAR a
 INPUT i
 EXEC
  a:=t
  write ( 'Writing a line...' )
 ENDEXEC
ENDMODEL
USE tc11 AS tc11
 INPUT i:=0
ENDUSE
RECORD tc11.a AS a
ENDMODELS
  NOD1                     1.00
BLANK card after last electric network branch
BLANK card ends switches
14NOD1        100.0       50.0                                   -1.0
BLANK card after last electric network source
  NOD1
C          Step      Time      NOD1        TACS                                                                                            
C                                            A                                                                                               
C            0       0.0        100.        0.0                                                                                           
C            1       0.1        100.        0.1                                                                                           
C            2       0.2        100.        0.2                                                                                           
C            3       0.3        100.        0.3                                                                                           
BLANK card ends selective node voltage outputs
C            9       0.9        100.        0.9                                                                                           
C           10       1.0        100.        1.0                                                                                           
C      Variable maxima :        100.        1.0                                                                                        
C      Times of maxima :         0.0        1.0                                                                                        
C      Variable minima :        100.        0.0                                                                                        
C      Times of minima :         0.0        0.0                                                                                        
BLANK card terminating plotting
BEGIN NEW DATA CASE
C         Illustration of MODELS using six of Laurent Dube's simple test cases           
C         6th of 16 data subcases.  This was taken from Laurent's file  TC13.DAT  
     0.1    1.00
       1       1       1       1       1      
MODELS
MODEL test
 INPUT i
 VAR a, b
 HISTORY b
 EXEC
  a:=b {MAX:0.7}
  b:=t
 ENDEXEC
ENDMODEL
USE test AS test
 INPUT i:=t
 HISTORY b:=t
ENDUSE
RECORD test.a AS a
ENDMODELS
  NOD1                     1.00
BLANK card after last electric network branch
BLANK card ends switches
14NOD1        100.0       50.0                                    0.0
BLANK card after last electric network source
  NOD1
C        Step      Time        NOD1       TACS                                                                                            
C                                           A                                                                                               
C           0       0.0         0.0        -.1                                                                                           
C           1       0.1        100.        0.0                                                                                           
C           2       0.2        100.        0.1                                                                                           
C           3       0.3        100.        0.2                                                                                           
C           4       0.4        100.        0.3                                                                                           
BLANK card ends selective node voltage outputs
C           9       0.9        100.        0.7                                                                                           
C          10       1.0        100.        0.7                                                                                           
C     Variable maxima :        100.        0.7                                                                                        
C     Times of maxima :         0.1        0.9                                                                                        
C     Variable minima :         0.0        -.1                                                                                        
C     Times of minima :         0.0        0.0                                                                                        
BLANK card terminating plotting
BEGIN NEW DATA CASE
C    7th of 16 subcases was added 22 August 1998 to illustrate a supplemental
C    variable that extends through column 80.  Furthermore, it is unusual in
C    that it has no blank byte to the right of the equal sign in column 11.
C    The result is easily checked since, by inspection, it has value  9 + T.
C    For background text, see the July, 1998, newsletter.  Dr. Yuan Bin called
C    our attention to this problem.  But that was not the end.  It was about a
C    year later that Carlos Mata pointed to the weakness of the correction: it
C    failed if column 80 contained a closed parenthesis.  The correction was
C    fortified March 18, 1998,  and new variable  CARLOS  confirms correct
C    operation.  Note that the values of  I1FG50  and  CARLOS  are equal.
ABSOLUTE TACS DIMENSIONS
      10      90     100      20      30     400     350      60
     .02     .10
       1       1       0       0       
TACS STAND ALONE
99I1FG00  = 1.
99I1FG02  = TIMEX
99I1FG04  = 1.
99I1FG06  = 1.
99I1FG08  = 1.
99I1FG10  = 1.
99I1FG12  = 1.
99I1FG14  = 1.
99I1FG16  = 1.
99I1FG18  = 1.
98YUAN    =I1FG00+I1FG02+I1FG04+I1FG06+I1FG08+I1FG10+I1FG12+I1FG14+I1FG16+I1FG18
98CARLOS  =I1FG00+I1FG02+I1FG04+I1FG06+I1FG08+I1FG10+I1FG12+ 1.5*(I1FG14+I1FG16)
33YUAN  CARLOS
BLANK card ending all TACS data cards
BLANK card ending plot cards
BEGIN NEW DATA CASE
C    8th of 16 subcases was added 18 August 1998 to illustrate new request words
C    that are related to repeatability of random numbers,  or lack of it.  Here
C    we illustrate repeatability, but the non-repeatable case really is the more
C    important.  If the  REPEATABLE  card below is replaced by the comment card
C    that immediately precedes it,  the new initialization of OVER8 (see the
C    October newsletter) is demonstrated.  Previously,  ATP GNU answers differed
C    from Salford and Watcom answers.  Also illustrated is no output variable,
C    with  MODELS  doing all output.  Of course,  with only 2 energizations and
C    a dummy electric network,  this is just a demonstrative shadow of the real
C    data case of interest (Gabor Furst's BACKFL.DAT).  The dominant structural
C    feature is use of  MODELS  for all output of a SYSTEMATIC data case.  There
C    _is_  a base case,  although the engineering significance is not known.
TRULY RANDOM NUMBERS     { Alternative to following card removes repeatability
REPEATABLE RANDOM NUMBERS     { English-language override of NSEED 5 lines below
C    Note about preceding:  Since the TRULY declaration is overridden by the
C    REPEATABLE  declaration,  in fact the  TRULY declaration does nothing.  But
C    it does illustrate acceptance by ATP.
PRINTED NUMBER WIDTH, 13, 2,  { Request maximum precision (for 8 output columns)
    .060    .300   { Time step is artificially large (5 steps to completion)
       1       1       1       2       1      -1                      -2
C    ISW   ITEST   IDIST    IMAX   IDICE  KSTOUT                           NSEED
       0       1       0       0       1                                       1
       5       5      20      20      { Printout frequency changes for base case
MODELS HYBRID
STORAGE integer pages:0 real pages:0 ENDSTORAGE
MODEL tc1
 VAR angle
 INPUT i
   INIT angle:=0
   ENDINIT
 EXEC
    angle:= random()
C   write('   T =',  t,  '  MODELS random() yields =', angle )  
--- Only program output within dT loop is preceding output of variable ANGLE ---
C     To preceding output,  add the energization number  KNT  to demonstrate
C     that a 4-byte integer is handled properly by XPR1 for  WRITE(.  GNU ATP
C     uses different code of IEIGHT for this than for ISTEP.  WSM, 28 Aug 2010
    write('   T =',  t,  'KNT =', atp(knt),  'MODELS random() ==>', angle )
 ENDEXEC
ENDMODEL
USE tc1 AS tc1a
 INPUT i:=0
ENDUSE
C RECORD tc1a.angle AS angle  { If activates,  variable angle becomes ATP output
ENDMODELS
  GENA  A1                    1. 
  A1                          1.                                             
  ASW1                        1.                                             
BLANK card ending branch cards
  A1    ASW1       4.E-3     .2E-3         2          SYSTEMATIC
BLANK card ending switch cards
14GENA     303.       60.         0.0                        -1.
BLANK card ending source cards
BLANK card ending output variable requests (node voltages only, here)
  PRINTER PLOT
BLANK card ending plot cards (for base case only)
BLANK card ending statistical tabulate (none possible, since no output variable)
BEGIN NEW DATA CASE
C   9th of 16 subcases is a meaningful case from Trainee Rod Price of BPA, who
C   assembled and tested the data during the summer  and  fall  of 1989.  The 
C   study  involves  the  500-kV Coulee-Raver overhead line,  and  the  model
C   includes series capacitors at Columbia.  These  capacitors  are protected 
C   by  nonlinear resistors  (ZnO modeling) that are monitored  and  bypassed 
C   using  MODELS.   To see more interesting simulation,  the  end-time  TMAX
C   should be increased to 100 msec.   Reformulated by Laurent Dube  8 Aug 91
C   following his correction of a MODELS bug for Bob Wilson at Univ. of Idaho
C   in Moscow.   It was observed that the answer of this 7-th subcase changed,
C   too (i.e., before it was wrong,  although perhaps not by much)!       
C   14 Sept 1993,  WSM  adds  $BLANK DATA  below to protect the user who has 
C   NOBLAN  of  STARTUP  equal to unity (equivalent to  $BLANK HALT  to halt if
C   a true blank card is found.   This  MODELS  data does use blank lines as 
C   comments,  and  MODELS  ignores the blanks.   The important thing is to
C   prevent EMTP data input from halting before  MODELS  receives the blanks!
C      14 October 2000,  the addition of one IF-statement to OVER16 (see 8th
C      subcase of DCNEW-16 and January, 2001, newsletter) changes the answers
C      of this subcase slightly.  To document this,  output variables are have
C      been made more relevant.  Differences follow the fault to ground.  From
C      the  .DBG  file,  the new logic is used 3 times:
C          OVER16 bypasses compensation.  J, T =  43   2.000000000000E-02
C          OVER16 bypasses compensation.  J, T =  43   3.150000000000E-02
C          OVER16 bypasses compensation.  J, T =  43   4.990000000000E-02
C      It is not surprising that the difference shows up after step 400, which
C      corresponds to that switch closure at T = 20 msec (1st of 3 uses).
C      The following comparison shows how extrema are affected only slightly:
C      Next   3     output variables are branch currents (flowing from the upper node to the lower node);
C      Next   5     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C      Step    Time       CR25A       CR20A       CRZ1A       MODELS      MODELS      MODELS      MODELS      MODELS
C                         CR30A       TERRA       CR20A       VCAPA       IZNA        PZNA        EZNA        GAPA
C      Old maxima :   27237.7816  15421.0962  .445513E-4  39549.6402  .445515E-4  .139454E10  .500879 E7         1.0
C      New maxima :   26931.1785  15312.52    .445513E-4  39549.6402  .445515E-4  .135464E10  .4929451E7         1.0
C      Times of both :    .03175       .0375       .0125       .0125       .0125        .028         .05       .0315
C      Old minima :   -7552.2441  -28033.65   -7291.149   -191594.73  -7291.149   .432745E-9         0.0         0.0
C      New minima :   -7466.2202  -28099.511  -7350.3614  -184341.43  -7350.3614  .432745E-9         0.0         0.0
C      Times of old :      .0378      .02865      .02815      .02815      .02815         0.0         0.0         0.0
C      Times of new :      .0378      .02865      .02805      .02805      .02805         0.0         0.0         0.0
C      About switching times, both have: Gap "CR25A "  to  "CR30A "  closing
C      after  3.15000000E-02  sec.  But openings differ by one time step:
C      Old:  Gap  "CR25A "  to  "CR30A "  opening after  4.98500000E-02  sec.
C      New:  Gap  "CR25A "  to  "CR30A "  opening after  4.99000000E-02  sec.
PRINTED NUMBER WIDTH, 12, 2,   { Cancel wider output of preceding subcase
$BLANK DATA       { Request old treatment of true blanks (not 
 .000050    .050     60.  
       1       9       1       3       1      -1
       5       5      20      20     100     100
C
MODELS HYBRID

INPUT -- declare names and function of all MODELS inputs from EMTP

        vc3a {v(CR30A)}   --   vc3a = voltage at node CR30A
        vc2a {v(CR20A)}   --   vc2a = voltage at node CR20A
        vc3b {v(CR30B)}   --   vc3b = voltage at node CR30B
        vc2b {v(CR20B)}   --   vc2b = voltage at node CR20B
        vc3c {v(CR30C)}   --   vc3c = voltage at node CR30C
        vc2c {v(CR20C)}   --   vc2c = voltage at node CR20C
        iznoa {i(CRZ2A)}  --  iznoa = current through switch at node CRZ2A
        iznob {i(CRZ2B)}  --  iznob = current through switch at node CRZ2B
        iznoc {i(CRZ2C)}  --  iznoc = current through switch at node CRZ2C

OUTPUT -- declare names of MODELS outputs to EMTP

        gapa, gapb, gapc  -- switch control signals

MODEL flash

comment ****************************************************************
*                                                                      *
*    Function: set or cancel the gap firing control signal             *
*    Inputs  : voltage and current across ZnO resistor                 *
*    Output  : the firing signal to the electrical ZnO component       *
*                                                                      *
************************************************************* endcomment

  INPUT  vcap   --  voltage difference across series caps  [Volts]
         iczn   --  ZNO current                            [Amps]

  DATA   Pset   --  power setting                          [Megajoules/msec]
         Eset   --  energy setting                         [Megajoules]
         fdel   --  firing delay                           [msec]
         fdur   --  firing duration                        [msec]
 
  VAR  power    --  power into ZnO resistor                [Watts]
       trip     --  gap firing control signal              [0 or 1]
       energy   --  energy into ZnO resistor               [Joules]
       tfire    --  time at which the gap was last fired   [sec]
 
  OUTPUT  trip

  HISTORY  INTEGRAL(power) {DFLT:0}

  INIT trip:=0
       tfire:=0
  ENDINIT

  EXEC
    ------------------------------------------------------------------
    power:=vcap*iczn
    energy:=INTEGRAL(power)
    ------------------------------------------------------------------
    IF  trip>0              -- is already firing
    AND t-tfire>fdur*1.e-3  -- has exceeded firing duration
    THEN
      trip:=0     -- cancel the firing signal
      tfire:=0    -- null the firing time
    ENDIF
    ------------------------------------------------------------------
    IF trip=0                        -- is not signaling to fire
    AND tfire=0                      -- firing condition not yet detected
    AND (    power  >= Pset * 1.e9   -- power setting exceeded
          OR energy >= Eset * 1.e6 ) -- energy setting exceeded
    THEN
      tfire:=t    -- set the firing detection time
    ENDIF
    ------------------------------------------------------------------
    IF trip=0                 -- is not signaling to fire
    AND tfire>0               -- firing condition has been detected
    AND t-tfire>=fdel*1.e-3   -- firing delay exceeded
    THEN
      trip:=1     -- set the firing signal
    ENDIF
    ------------------------------------------------------------------
  ENDEXEC

ENDMODEL  

USE flash AS flasha   ----------------------------------- PHASE A ------

  TIMESTEP MIN: 0.0005

  INPUT  vcap:=vc3a-vc2a   -- voltage across caps
         iczn:=iznoa       -- current through ZnO resistor

  DATA   Pset:= 1   -- power setting    [MJ/mS]
         Eset:= 9   -- energy setting   [MJ]
         fdel:= 4   -- firing delay     [msec]
         fdur:=20   -- firing duration  [msec]

  OUTPUT  gapa:=trip  -- passes value of 'trip' to MODELS variable 'GAPA'

ENDUSE -----------------------------------------------------------------

USE flash AS flashb   ----------------------------------- PHASE B ------

  TIMESTEP MIN: 0.0005

  INPUT  vcap:=vc3b-vc2b   -- voltage across caps
         iczn:=iznob       -- current through ZnO resistor

  DATA   Pset:= 1   -- power setting    [MJ/mS]
         Eset:= 9   -- energy setting   [MJ]
         fdel:= 4   -- firing delay     [msec]
         fdur:=20   -- firing duration  [msec]

  OUTPUT  gapb:=trip  -- passes value of 'trip' to MODELS variable 'GAPB'

ENDUSE -----------------------------------------------------------------

USE flash AS flashc   ----------------------------------- PHASE C ------

  TIMESTEP MIN: 0.0005

  INPUT  vcap:=vc3c-vc2c   -- voltage across caps
         iczn:=iznoc       -- current through ZnO resistor

  DATA   Pset:= 1   -- power setting    [MJ/mS]
         Eset:= 9   -- energy setting   [MJ]
         fdel:= 4   -- firing delay     [msec]
         fdur:=20   -- firing duration  [msec]

  OUTPUT  gapc:=trip  -- passes value of 'trip' to MODELS variable 'GAPC'

ENDUSE -----------------------------------------------------------------

RECORD  flasha.vcap   AS VCAPA       
        flasha.iczn   AS IZNA
        flasha.power  AS PZNA       
        flasha.energy AS EZNA
        flasha.trip   AS GAPA

ENDMODELS      
C  ---------------------   End  MODELS  data,  begin electric network:  
C   **** Series Capacitors, Xc = 28.4 ohms
C
  CR20A CR30A                         93.40                                     
  CR20B CR30B                         93.40                                     
  CR20C CR30C                         93.40                                     
C
  CR20D CR30D                         93.40                                    
  CR20E CR30E                         93.40                                    
  CR20F CR30F                         93.40                                    
C
C   **** Jumpers to help EMTP convergence with ZNO
C
  CRZ1A CRZ2A              0.01                                                 
  CRZ1B CRZ2B              0.01                                                 
  CRZ1C CRZ2C              0.01                                                 
C
  CRZ1D CRZ2D              0.01                                                 
  CRZ1E CRZ2E              0.01                                                 
  CRZ1F CRZ2F              0.01                                                 
C
C   **** ZNO protection for Caps, I=1.8KA, 2.2 PU, alpha = 40, Vref =147.5
C
92CRZ1A CR20A                          5555.                                   1
                  147500.                      -1.                             
                      1.0                      40.                      .80
                    9999.
92CRZ1B CR20B CRZ1A CR20A              5555.                                  
92CRZ1C CR20C CRZ1A CR20A              5555.                                  
C
92CRZ1D CR20D CRZ1A CR20A              5555.                                   
92CRZ1E CR20E CRZ1A CR20A              5555.                                   
92CRZ1F CR20F CRZ1A CR20A              5555.                                   
C
C   **** Source Impedances
C
51SRC1A RAVBA               0.1  20.0
52SRC1B RAVBB               0.1  18.5
53SRC1C RAVBC
C
51SRC1A RAVBA              300.0
52SRC1B RAVBB              150.0
53SRC1C RAVBC
C
C
51SRC2A GRCBA               0.1  7.00
52SRC2B GRCBB               0.1  10.7
53SRC2C GRCBC
C
51SRC2A GRCBA              350.0
52SRC2B GRCBB              150.0
53SRC2C GRCBC
C
C  *** LINE IMPEDANCES  Values have been 1/4 to get high fault MVA ****
C
51CR50A CR30A               0.5  14.0                                          
52CR50B CR30B                8.  52.0                                          
53CR50C CR30C                                                                  
C
51CR50D CR30D               0.5  14.0                                          
52CR50E CR30E                8.  52.0
53CR50F CR30F
C
51CR01A CR20A               0.4   9.0
52CR01B CR20B                5.  36.0
53CR01C CR20C
C
51CR01D CR20D               0.4   9.0                                       
52CR01E CR20E                5.   36.
53CR01F CR20F
C
C  *****  BYPASS IMPEDANCE ****
  CR20A CR25A              5.    0.23                                           
  CR20B CR25B              5.    0.23                                           
  CR20C CR25C              5.    0.23                                           
  CR20D CR25D              5.    0.23                                          
  CR20E CR25E              5.    0.23                                         
  CR20F CR25F              5.    0.23                                         
C   ****  PARALLEL RESISTOR
  CR20A CR25A              200.                                                
  CR20B CR25B              200.                                                 
  CR20C CR25C              200.                                                 
  CR20D CR25D              200.                                                
  CR20E CR25E              200.                                               
  CR20F CR25F              200.                                               
BLANK card ending program branch cards.
C  **** TACS CONTROLLED CAP BYPASS SWITCHES *******
C  ** SWITCH 1
12CR25A CR30A                                                   GAPA          11
12CR25B CR30B                                                   GAPB          1 
12CR25C CR30C                                                   GAPC          1 
C 12CR25D CR30D                                                   GAPD           1
C 12CR25E CR30E                                                   GAPE            
C 12CR25F CR30F                                                   GAPF            
  CR25D CR30D  -1.0        10.
  CR25E CR30E  -1.         10.
  CR25F CR30F  -1.         10.      
C
C   **** Main Breaker Switches
C
  RAVBA CR50A  -1.0      10.0                                                   
  RAVBB CR50B  -1.0      10.0                                                   
  RAVBC CR50C  -1.0      10.0                                                 
C
  RAVBA CR50D  -1.0      10.0                                                   
  RAVBB CR50E  -1.0      10.0
  RAVBC CR50F  -1.0      10.0
C
  GRCBA CR01A  -.006     10.0                                                 
  GRCBB CR01B  -.006     10.0                                                 
  GRCBC CR01C  -.006     10.0                                                 
C
  GRCBA CR01D  -1.0      10.0
  GRCBB CR01E  -1.0      10.0
  GRCBC CR01F  -1.0      10.0
C
C  *** ZNO current measuring switches  ***
C
  CR30A CRZ2A        -1.        1.
  CR30B CRZ2B        -1.        1.
  CR30C CRZ2C        -1.        1.
  CR30D CRZ2D        -1.        1.
  CR30E CRZ2E        -1.        1.
  CR30F CRZ2F        -1.        1.
  CR20A          0.01998       10.   { Fault switch, phase "a" to ground }     1
BLANK card terminating program switch cards
14SRC1A    440000.   60.       -20.0                          -1.      10.
14SRC1B    440000.   60.       100.0                          -1.      10.
14SRC1C    440000.   60.       220.0                          -1.      10.
14SRC2A    440000.   60.         0.0                          -1.      10.
14SRC2B    440000.   60.       120.0                          -1.      10.
14SRC2C    440000.   60.       240.0                          -1.      10.
C --------------+------------------------------
C From bus name |  Names of all adjacent busses.
C --------------+------------------------------
C        CR20A  |TERRA *CR30A *CRZ1A *CR01A *CR25A *CR25A *
C        CR30A  |CR20A *CRZ2A *CR50A *CR25A *
C        CR20B  |CR30B *CRZ1B *CR01B *CR25B *CR25B *
C        CR30B  |CR20B *CRZ2B *CR50B *CR25B *
C        CR20C  |CR30C *CRZ1C *CR01C *CR25C *CR25C *
C        CR30C  |CR20C *CRZ2C *CR50C *CR25C *
C        CR20D  |CR30D *CRZ1D *CR01D *CR25D *CR25D *
C        CR30D  |CR20D *CRZ2D *CR50D *CR25D *
C       Total network loss  P-loss  by summing injections =   8.848099600964E+07
C Output for steady-state phasor switch currents.
C      Node-K    Node-M           I-real            I-imag            I-magn  
C      CR25A     CR30A            Open              Open              Open   
C      CR25B     CR30B            Open              Open              Open   
C      CR25C     CR30C            Open              Open              Open      
C      CR25D     CR30D        9.29172435E+02   -1.77481409E+02    9.45970963E+02
C      CR25E     CR30E       -3.10882808E+02    8.93427637E+02    9.45970963E+02
C      CR25F     CR30F       -6.18289626E+02   -7.15946228E+02    9.45970963E+02
BLANK card terminating program source cards.
C End inject:  SRC2C   -220000.       440000.   -1191.301327068  2345.3572211208
C End inject:   -381051.1776652  -120.0000000   -2020.272665456     -120.5267278
C Step   Time       CR25A      CR25A      TACS       TACS       TACS        TACS
C                   CR30A      CR30A      VCAPA      IZNA       PZNA        EZNA
C ***      Phasor I(0) =  9.2917243E+02   Switch  "CR25D "  to  "CR30D "  closed
C ***      Phasor I(0) = -3.1088281E+02   Switch  "CR25E "  to  "CR30E "  closed
C ***      Phasor I(0) = -6.1828963E+02   Switch  "CR25F "  to  "CR30F "  closed
C            < <  Etc.  (many switch currents)  > >
C   0     0.0  -143.44999        0.0 143.449993 .28697E-11 .411659E-9        0.0
C   1   .5E-4  602.099261        0.0 -602.09926 -.68002E-6 .411659E-9        0.0
C   2   .1E-3  1347.43459        0.0 -1347.4346 -.15198E-5 .411659E-9        0.0
C   3  .15E-3   2092.2912        0.0 -2092.2912 -.23592E-5 .411659E-9        0.0
BLANK card ending output variables requests (none here, since all column 80)
C   300  .015  -23363.319        0.0 23363.3193 .263161E-4 .614832009 .014309206
C ***          Close switch  "CR20A "  to  "      "  after  2.00000000E-02  sec.
C   400   .02  37570.0647        0.0 -37570.065 -.42322E-4 1.59004745 .016928315
C   500  .025  56444.8022        0.0 -56444.802 -.63583E-4 3.58891722 .022232887
C   600   .03  183226.734        0.0 -183226.73 -5786.1028 .106017E10 .3650637E7
C            Gap     "CR25A "  to  "CR30A "  closing after  3.15000000E-02  sec.
C   700  .035         0.0 -4433.8226 22742.1728 .256182E-4 .582612436 .4929451E7
C   800   .04         0.0 -5450.0316 25576.0355 .288104E-4 .736855415 .4929451E7
C   900  .045         0.0 4627.82167  -23165.41 -.26095E-4 .604499203 .4929451E7
C            Gap     "CR25A "  to  "CR30A "  opening after  4.99000000E-02  sec.
C  1000   .05   -816.5929        0.0 3143.46664 .354099E-5 .011130982 .4929451E7
C Variable max:184341.425 26931.1785 39549.6402 .445516E-4 .135464E10 .4929451E7
C Times of max:    .02805     .03175      .0125      .0125       .028        .05
C Variable min: -39549.64 -7466.2202 -184341.43 -7350.3614 .411659E-9        0.0
C Times of min:     .0125      .0378     .02805     .02805        0.0        0.0
  PRINTER PLOT                                                                  
 194 5. 0.0 40.         CR25A CR30A             { Axis limits: (-0.747,   2.220)
BLANK card ending plot cards
BEGIN NEW DATA CASE
PRINTED NUMBER WIDTH, 12, 2, { Keep dT loop precision the same, but 2 blank separators
C     10th of 16 subcases was added 10 April 2011  to illustrate very simple use
C     of  POINTLIST  within  MODELS.  Data comes from Orlando Hevia of Santa Fe,
C     Argentina.  It  establishes the standard of comparison for his alternative
C     SEEK  function (see 11th subcase,  which produces the same answer).  Here
C     the (X,Y) points of  POINTLIST  number only 7.  For much larger numbers of
C     points (e.g., 20200 for Orlando's TFGNMOD.DAT), Dube's  POINTLIST  is much
C     too slow.  For much larger numbers of points,  use Orlando's new  SEEK  as
C     will be illustrated by the following subcase.  WSM.
     1.0   10.00     { Take ten steps, from 1 second to 10 seconds.  This is "T"
       1      -1   { Print every step;  do not bother saving points for plotting
MODELS STAND ALONE
MODEL SAMPLE
VAR
   X,Y,Z
FUNCTION R POINTLIST
(0, 69)
(2, 50)
(4, 142)
(5, 142)
(6, 188)
(8, 73)
(10, 100)
INIT
        X:=0
ENDINIT
EXEC
X:=R(T)
Y:=SQRT(X)*0.3
Z:=R(SQRT(X)*0.3)
ENDEXEC
ENDMODEL
USE SAMPLE AS SAMPLE
ENDUSE
RECORD
  SAMPLE.X AS X
  SAMPLE.Y AS Y
  SAMPLE.Z AS Z
ENDRECORD
ENDMODELS
C Column headings for the  3   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  Next   3     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C   Step      Time       MODELS      MODELS      MODELS
C                        X           Y           Z
C      0       0.0          69.  2.49198716  72.6314093
C      1       1.0         59.5  2.31408729  64.4480155
C      2       2.0          50.  2.12132034  55.5807358
C      3       3.0          96.  2.93938769  93.2118338
C      4       4.0         142.  3.57491259  122.445979
C      5       5.0         142.  3.57491259  122.445979
C      6       6.0         188.  4.11339276  144.608033
C      7       7.0        130.5   3.4270979  115.646503
C      8       8.0          73.  2.56320112  75.9072517
C      9       9.0         86.5  2.79016129  86.3474191
C     10       10.         100.         3.0         96.
BLANK CARD ENDING PLOT CARDS
BEGIN NEW DATA CASE
$DEPOSIT, JSEEDR=0      { Reinitialize random # generator as if no preceding use
C        Note: the preceding $DEPOSIT  is required because random # use here is
C              not the 1st of data case.  Preceding subcases roll dice.  Memory
C              of this is erased by setting JSEEDR to zero.   WSM.  29 July 2011
C              This is for addition of GAUSS.  Variable  NORMAL  below should be
C              identical in value to what is produced by TACS within subcase 2
C              of BENCHMARK DC-18.  This works well for this subcase,  but later
C              subcases then would have different random numbers.  In order not
C              to change the random numbers of following subcases,  we want to
C              restore at the end the JSEEDR value that exists at the beginning.
C              New logic of DEPOSI makes this possible.  See 2nd $DEPOSIT below.
$PREFIX, []    { $INCLUDE files are located in same place as this main data file
C     Preceding line allows the SEEK file to be remote if this main data file is
PRINTED NUMBER WIDTH, 12, 2, { Keep dT loop precision the same, but 2 blank separators
C     11th of 16 subcases was added  10 April 2011  to illustrate that Orlando
C     Hevia's new  SEEK  function of  MODELS  will produce the same answer as
C     the preceding 10th subcase,  which uses  POINTLIST.  Data is too small for
C     there to be a reason to use  SEEK  instead of  POINTLIST,  of course.  For
C     a more realistic illustration of the advantage of  SEEK,  see explanation 
C     of the EEUG list server dated 11 June 2011.  This includes the following 
C     paragraph:
C   The inadequacy of  MODELS  POINTLIST  for large data (many points of the 
C   piecewise-linear X-Y function) was documented in detail in a paper that 
C   was presented at the year-2006 EEUG conference.  See  "Determination of 
C   an optimized energy storage size for a wind farm based on wind forecasts"  
C   by Steve Völler,  Christian Müller, and Johannes Verstege of the 
C   University of Wuppertal in Germany.  In this case,  POINTLIST  functioned 
C   correctly, but was too slow to be practical.  Those German users give new 
C   meaning to the term industrial-strength.  Would a reader believe 105K 
C   pairs of (X, Y)?  Yes, letter K as in one thousand.  Mr. Dube obviously 
C   did not imagine such use when he programmed his original logic.  But 
C   Orlando Hevia's new  SEEK  now provides an efficient alternative without 
C   the need for compilation and linking of user-supplied source code.
C     This data too comes from Orlando Hevia of Santa Fe, Argentina.  WSM.
     1.0   10.00     { Take ten steps, from 1 second to 10 seconds.  This is "T"
       1      -1   { Print every step;  do not bother saving points for plotting
MODELS STAND ALONE
C 3456789012345678901234567890123456789012345678 -- Ruler for following optional 
C                 < File Name>    IPRMDL  KOMPAC    miscellaneous data card :
MODELS MISC. DATA dc68seek.            0       0
C               Note about preceding.  Special value  KOMPAC = 8765  is required
C               to bypass the check on (X,Y) points of the function.  Logic will
C               ensure that X never decreases (a fatal error).  It also discards
C               all but the end points of each string of three or more identical
C               Y values.  That is the 3rd of 3 params shown.  As for the first,
C               see the substantial note below.  Finally,  the middle parameter,
C               IPRMDL,  is just  IPRSUP  for the new  MODELS  code.   I.e.,  it
C               provides diagnostic printout control.   All 3 values will remain
C               in force unless redefined.  So,  for stacked subcases,  there is
C               no need for repeated definition unless values are to be changed.
MODEL SAMPLE
VAR X,Y,Z, NORMAL
FUNCTION seek FOREIGN seek {ixarg: 2}
FUNCTION gauss FOREIGN gauss {ixarg: 2}
EXEC
C    Of the three following  EXEC  lines,  the first and the third use Orlando's
C    new  SEEK  function which is defined by the  (X, Y)  pairs of  DC68SE37.DAT
C    where the root name  "DC68SEEK."  was read from the MODELS misc. data card.
C    The  "37"  that overwrites the  "EK"  of DC68SEEK corresponds to the second
C    argument of  SEEK.  If another function were involved,  it would correspond
C    to a different subscript (any # < 100 million is legal).  An arbitrary # of
C    such functions is allowed,  and  each can involve an arbitrary # of  (X, Y)
C    as long as List 20 has adequate unused space.  As for  .DAT  being the file
C    type,  this follows from  DATTYP  of  STARTUP  because the user's file name
C    ended with a period.  If there is no MODELS miscellaneous data, the default
C    disk file name  SEEKDATA.DAT  will be assumed by ATP.  WSM.
X:=SEEK(T, 37 )
Y:=SQRT(X)*0.3
Z:=SEEK(SQRT(X)*0.3, 37 )
C    In order that all data be seen,  let's document the content of DC68SE37.DAT
C    which contains just the following 7 formatted lines:
C 0, 69
C 2, 50
C 4, 142
C 5, 142
C 6, 188
C 8, 73        
C 10, 100
C    This is just usual free-format data.  Each line begins in column 1, and
C    the final line is known as ATP encounters an end-of-file mark.  Neither
C    comment lines nor $-cards nor inline comments are allowed.  This is not
C    an oversight.  Rather,  it is a way to speed execution for large files.
C    If the user wants to terminate his function prematurely,  however,  he
C    is allowed to add a line such as  "-9999., 0"  where the  -9999  is
C    interpreted as a software EOF.  For an illustration,  see subcase # 15.
C       For some reason, Dube does not accept an in-line comment on the
C       following function.  Perhaps not on any function?  So, strip it off:
C NORMAL:=gauss( T )  { Argument "T" in fact is unused, so is arbitrary
NORMAL:=gauss( T )
ENDEXEC
ENDMODEL
USE SAMPLE AS SAMPLE
ENDUSE
RECORD
SAMPLE.X AS X
SAMPLE.Y AS Y
SAMPLE.Z AS Z
SAMPLE.NORMAL AS NORMAL
ENDRECORD
ENDMODELS
C   MODELS "SEEK" function # 37 involves 7 pairs of (X,Y).  These are stored in List-20 SCONST cells 3 through 16.
C Column headings for the  4   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  Next   4     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C   Step      Time       MODELS      MODELS      MODELS      MODELS
C                        X           Y           Z           NORMAL
C      0       0.0          69.  2.49198716  72.6314093  .890624441
C      1       1.0         59.5  2.31408729  64.4480155  -.26382016
C      2       2.0          50.  2.12132034  55.5807358  -.41782888
C      3       3.0          96.  2.93938769  93.2118338   .66418476
C      4       4.0         142.  3.57491259  122.445979  .242122394
C      5       5.0         142.  3.57491259  122.445979  -.38464212
C      6       6.0         188.  4.11339276        142.   -.9889959
C      7       7.0        130.5   3.4270979  115.646503  -.26589058
C      8       8.0          73.  2.56320112  75.9072517  -.73063403
C      9       9.0         86.5  2.79016129  86.3474191  -.75107385
C     10       10.         100.         3.0         96.  .288485151
C   Size 11-20:     0     3 -9999 -9999  1744     0     0     0    23    15
C   The preceding line from 80-column case-summary statistics confirms that
C   15 cells of List 20 storage are being used by  SEEK.  There are 6 points
C   for a total of 12 values.  But for each function,  there are two extra
C   cells.  Finally,  the last function must be bounded by a blank,  so the
C   total is  2 * 6 + 2 + 1 = 15  as shown.
C       That was before adding  GAUSS  use.  After the addition, this becomes:
C   Size 11-20:     0     4 -9999 -9999  1838     0     0     0    23    17
$DEPOSIT, JSEEDR= -98789  { Flag for DEPOSI to restore value before previous use
C   Note about preceding line:  As this use begins on  31 July 2011,  service is
C   limited to a single variable (in this case,  JSEEDR).  For use here,  JSEEDR
C   is a random variable.  The user wants to restore whatever value  JSEEDR  had
C   at the instant it was zeroed by the preceding $DEPOSIT.  The user can not do
C   this manually (old logic)  because he does not know what the value was.  But    
C   ATP remembers the value within  DEPOSI,  and can do it.   WSM.  31 July 2011
BLANK CARD ENDING PLOT CARDS
BEGIN NEW DATA CASE
C     12th of 16 subcases  was  added  13 April 2011  to illustrate that Orlando
C     Hevia's new  EGAUSS  function of MODELS will perform Gaussian elimination.
C     This was the 1st of 4 stacked subcases within Orlando's  file  TFGNMOD.DAT
C     Solve the following 3 linear equations in 3 unknowns:
C        X1 + 4*X2 +   X3 = 7
C        X1 + 6*X2 -   X3 = 13
C      2*X1 -   X2 + 2*X3 = 5
C     THE INITIAL SOLUTION IS X1= 5, X2= 1, X3= -2,  x4=18
C     THE ELEMENT N+1 OF SOLUTION IS THE DETERMINANT OF the coefficient matrix
     0.2     1.0   { Take 5 time steps, changing the matrix at each new step
       1      -1   { Suppress accumulation of plot points,  which are unused
MODELS STAND ALONE
C 3456789012345678901234567890123456789012345678 -- Ruler for following optional 
C                 < File Name>    IPRMDL  KOMPAC    miscellaneous data card :
MODELS MISC. DATA dc68seek.            0       0
C THE FUNCTION CAN BE USED AS FOLLOW:
C
C SOLUTION[1..N+1]:= EGAUSS(N,MATRIX[1..N*N+N])
C
C EGAUSS FROM GAUSSIAN ELIMINATION
C
MODEL solve
FUNCTION egauss FOREIGN egauss {ixarg:100}
 VAR arg[1..12]
 VAR x[1..4]
EXEC
if t=0 then
 arg[1]:=1.0
 arg[2]:=4.0
 arg[3]:=1.0
 arg[4]:=1.0
 arg[5]:=6.0
 arg[6]:=-1.0
 arg[7]:=2.0
 arg[8]:=-1.0
 arg[9]:=2.0
 arg[10]:=7.0
 arg[11]:=13.0
 arg[12]:=5.0
endif
x[1..4]:=egauss(3,arg[1..12])
write( 'Solution vector X =',  x[1], x[2], x[3],  '  Determinant =', x[4] )
arg[1]:= arg[1]*2.0
arg[12]:= arg[12]+1.0
C   Warning about the following use of RANDOM,  which is the  MODELS  connection
C   to ATP's random number generator.  For  KOMPAR = 0  or 1  (set in  STARTUP),
C   answers will be different  every time  the data is simulated,  and different
C   program versions can not be compared  (e.g., GNU vs. Salford).  For purposes
C   of a test case, we want repeatability and comparability, so set  KOMPAR = 3.
C   Also, position within the data case (here, the 12th subcase) may be critical
C   since future random numbers depend on past random numbers.
arg[2]:=random()
C arg[2]:= arg[1] * arg[12] --- a deterministic alternative to preceding line
write('arg(1,12,2) =', arg[1], arg[12], arg[2] )
ENDEXEC
ENDMODEL
USE solve AS solve
ENDUSE
ENDMODELS
C    Because there is no electric network,  there is no usual output of the
C    time-step loop.  Once the dT loop begins,  the only output is from MODELS:
C  Blank card ends electric sources.  KCONST = 1     |BLANK  card ending dummy source that ATP added to  TACS  or  MODELS STAND ALONE
C  X vector = 5.0         1.0         -2.          Determinant = 18.
C  arg(1,12,2) = 2.0         6.0         .2328306E-9
C  Blank card ending requests for output variables.  |BLANK  card ending output requests (none) that ATP added to this  STAND ALONE
C  Column headings for the  0   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  X vector = 3.111111111 1.777777778 .7777777784  Determinant = 8.999999999
C  arg(1,12,2) = 4.0         7.0         .1608161E-4
C  X vector = .7096617728 2.741941174 4.161308814  Determinant = 30.99993567
C  arg(1,12,2) = 8.0         8.0         .1107408979
C  X vector = .1640999973 3.031236365 5.351518185  Determinant = 74.55703641
C  arg(1,12,2) = 16.         9.0         .763080108
C  X vector = -.104457932 3.219802884 6.214359374  Determinant = 159.9476796
C  arg(1,12,2) = 32.         10.         .1799780347
C  X vector = -.007328865 3.275392315 6.645025023  Determinant = 338.2800879
C  arg(1,12,2) = 64.         11.         .902878134
BLANK card ending plot requests
BEGIN NEW DATA CASE
C     13th of 16 subcases  has the same solution as the preceding.  It was added
C     30 April 2011 to illustrate the  LINSOL  alternative to  EGAUSS.  The data
C     is compatible except that no value for the determinant  of the coefficient
C     matrix is available.  About history,  SUBROUTINE LINSOL  was provided by
C     MODELS author Laurent Dube,  paid by BPA,  as part of the original user-
C     supplied source code of MODELS.  This dates to 1995 or 1996 (not part of
C     the original MODELS).  Comment cards document the source and use:
C        Solves [a]*x=b using Crout's method with partial pivoting
C        Reference: Numerical Recipes chapter 2.
C        Modifies [a] and b, and returns b=x
C        Returns ierr=1 if [a] is singular, else returns ierr=0
C        Uses divzro as the smallest possible divisor value
     0.2     1.0   { Take 5 time steps, changing the matrix at each new step
       1      -1   { Suppress accumulation of plot points,  which are unused
MODELS STAND ALONE
MODEL solve
FUNCTION linsol FOREIGN linsol {ixarg:100}
 VAR arg[1..12]
 VAR x[1..4]
EXEC
if t=0 then
C     As written,  LINSOL  has been connected using the interface of Orlando's
C     EGAUSS.  There is a complication, however.  EGAUSS  stored the matrix by
C     rows whereas  LINSOL  assumes storage by columns.  Note the difference:
C        1    2    3     4    5     6     7    8    9   --  Index of  MODELS
C       1,1  1,2  1,3   2,1  2,2,  2,3   3,1  3,2  3,3  --  Content for EGAUSS
C       1,1  2,1  3,1   1,2  2,2,  3,2   1,3  2,3  3,3  --  Content for LINSOL
C     I.e.,  one is the transpose of the other.  Of the 9 elements that must
C     be assigned,  only the 3 diagonal elements are unchanged.  One at a time:
  arg[1]:=1.0 -- Element (1,1) of EGAUSS is a diagonal, so remains unchanged
C arg[2]:=4.0  { Element (1,2) of EGAUSS had index 2,  so index 4 of LINSOL:
  arg[4]:=4.0
C arg[3]:=1.0   { Element (1,3) of EGAUSS had index 3,  so index 7 of LINSOL
  arg[7]:=1.0
C arg[4]:=1.0   { Element (2,1) of EGAUSS had index 4,  so index 2 of LINSOL
  arg[2]:=1.0
  arg[5]:=6.0 -- Element (2,2) of EGAUSS is a diagonal, so remains unchanged
C arg[6]:=-1.0  { Element (2,3) of EGAUSS had index 6,  so index 8 of LINSOL
  arg[8]:=-1.0
C arg[7]:=2.0   { Element (3,1) of EGAUSS had index 7,  so index 3 of LINSOL
  arg[3]:=2.0
C arg[8]:=-1.0  { Element (2,3) of EGAUSS had index 8,  so index 6 of LINSOL
  arg[6]:=-1.0
  arg[9]:=2.0 -- Element (3,3) of EGAUSS is a diagonal, so remains unchanged
C   That was for the 3x3 square matrix. The right-hand-side vector is unchanged:
 arg[10]:=7.0
 arg[11]:=13.0
 arg[12]:=5.0
endif
x[1..4]:=linsol(3,arg[1..12])
write( 'Solution vector X =',  x[1], x[2], x[3] )
arg[1]:= arg[1]*2.0
arg[12]:= arg[12]+1.0
C arg[2]:=random()  { Element (1,2) of EGAUSS had index 2,  so index 4 of LINSOL
  arg[4]:=random()
C write('arg(1,12,2) =', arg[1], arg[12], arg[2] )
write(  'arg(1,12,4) =', arg[1], arg[12], arg[4] )
ENDEXEC
ENDMODEL
USE solve AS solve
ENDUSE
ENDMODELS
C    Because there is no electric network,  there is no usual output of the
C    time-step loop.  Once the dT loop begins,  the only output is from MODELS:
C  Blank card ends electric sources.  KCONST = 1     |BLANK  card ending dummy source that ATP added to  TACS  or  MODELS STAND ALONE
C  X vector = 5.0         1.0         -2.
C  arg(1,12,4) = 2.0         6.0         .2328306E-9
C  Blank card ending requests for output variables.  |BLANK  card ending output requests (none) that ATP added to this  STAND ALONE
C  Column headings for the  0   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  X vector = 3.111111111 1.777777778 .7777777784
C  arg(1,12,4) = 4.0         7.0         .1608161E-4
C  X vector = .7096617728 2.741941174 4.161308814
C  arg(1,12,4) = 8.0         8.0         .1107408979
C  X vector = .1640999973 3.031236365 5.351518185
C  arg(1,12,4) = 16.         9.0         .763080108
C  X vector = -.104457932 3.219802884 6.214359374
C  arg(1,12,4) = 32.         10.         .1799780347
C  X vector = -.007328865 3.275392315 6.645025023
C  arg(1,12,4) = 64.         11.         .902878134
BLANK card ending plot requests
BEGIN NEW DATA CASE
C     14th of 16 subcases  was  added  17 April 2011  to illustrate that Orlando
C     Hevia's new  INVERT  function of  MODELS  will indeed invert a matrix.
C     Data comes from Orlando Hevia of Santa Fe, Argentina.
     1.0     5.0   { Take 5 time steps, changing the matrix at each new step
       1      -1   { Suppress accumulation of plot points,  which are unused
MODELS STAND ALONE
C    The input matrix is on the left,  its inverse is on the right:
C   ORIGINAL matrix           INVERSE CALCULATED WITH QPRO
C  1        2        3        -0.33333  0.151515  0.121212
C  4        2        5        0.666667  -0.48485  0.212121
C  6        3        2          -2E-17  0.272727  -0.18182
C
C DEFINE ALL THE VARIABLES TO BE PASSED TO INVERT AS OUTPUT
C    Since the matrix to be inverted is 3x3,  9 scalar variables are used :
VAR a1, a2, a3, a4, a5, a6, a7, a8, a9
C    Next come 4 scalar dimensions of the call to inversion routine INVERT:
C IXDATA: NUMBER ELEMENTS IN AN ARRAY XDATA OR ARGUMENTS, FOR EXAMPLE, THE
C         NUMBER OF EQUATIONS, A KEY TO INVERT OR TO CALCULATE DETERMINANT
C
C IXIN:   NUMBER OF ELEMENTS IN THE INPUT ARRAY XIN. SQUARE OF NUMBER OF
C         EQUATIONS IN THIS CASE
C
C IXOUT:  NUMBER OF ELEMENTS IN THE OUTPUT ARRAY XOUT. SQUARE OF NUMBER OF
C         EQUATIONS IN THIS CASE, 1 FOR DETERMINANT (NOt IMPLEMENTED)
C
C IXVAR:  NUMBER OF ELEMENTS IN HISTORY ARRAY, NOt USED IN THIS CASE
C
C    The following is the call to invert the matrix.  IXIN  is # of elements (9)
C    of the input matrix,  IXOUT  is the # of elements (9) of the output matrix.
C    The 1st and the 4th parameters  (IXDATA and IXVAR, respectively)  do not
C    depend on the matrix order (3 in this case).
MODEL invert FOREIGN invert {ixdata:1, ixin:9, ixout:9, ixvar:0}
USE invert AS invert
 DATA xdata[1]:=3   -- Order of the matrices involved (here, 3)
C    The following defines the 9 cells of the input matrix XIN.  Storage is by
C    rows. Note the trailing "t" which adds simulation time to each element of
C    the matrix at each time step :
 INPUT xin[1..9] := [ 1.0, 2.0, 3.0,
                      4.0, 2.0, 5.0,
                      6.0, 3.0, 2.0 ] + t
C    The following connects the 9 elements of the output matrix (the inverse)
C    with variable names  A1  through  A9.  This is done to interface with the
C    ATP output vector (the preceding subcase had no such connection, note) :
 OUTPUT a1:=xout[1],a2:=xout[2],a3:=xout[3],
        a4:=xout[4],a5:=xout[5],a6:=xout[6],
        a7:=xout[7],a8:=xout[8],a9:=xout[9]
ENDUSE
RECORD
       a1 AS a11
       a2 AS a12
       a3 AS a13
       a4 AS a21
       a5 AS a22
       a6 AS a23
       a7 AS a31
       a8 AS a32
       a9 as a33
ENDMODELS
C Column headings for the  9   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  Next   9     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C   Step    Time       MODELS      MODELS      MODELS      MODELS      MODELS      MODELS      MODELS      MODELS      MODELS
C                      A11         A12         A13         A21         A22         A23         A31         A32         A33
C      0     0.0   -.33333333  .151515152  .121212121  .666666667  -.48484848  .212121212         0.0  .272727273  -.18181818
C      1     1.0   -.31914894   .14893617  .127659574  .574468085  -.46808511  .170212766   -.0212766  .276595745  -.19148936
C      2     2.0   -.31147541  .147540984  .131147541  .524590164  -.45901639  .147540984  -.03278689  .278688525  -.19672131
C      3     3.0   -.30666667  .146666667  .133333333  .493333333  -.45333333  .133333333        -.04         .28         -.2
C      4     4.0   -.30337079  .146067416  .134831461  .471910112   -.4494382  .123595506  -.04494382  .280898876  -.20224719
C      5     5.0   -.30097087  .145631068   .13592233   .45631068  -.44660194  .116504854  -.04854369  .281553398   -.2038835
BLANK card ending plot requests
BEGIN NEW DATA CASE
$PREFIX, []    { $INCLUDE files are located in same place as this main data file
C     Preceding line allows the SEEK file to be remote if this main data file is
PRINTED NUMBER WIDTH, 12, 2, { Keep dT loop precision the same, but 2 blank separators
C     15th of 16 subcases  was added  21 April 2011  to illustrate extensions to
C     Orlando Hevia's  SEEK  function that were not tested by the preceding 11th
C     subcase.  Data is similar.  That is, function # 82 began as a copy of # 32  
C     but is more general in that two redundant (X, Y) points have been added
C     (to be omitted by ATP) and the software EOF is used, permitting arbitrary
C     comment information at the bottom of  DC68SE82.DAT  Erasure of a no longer
C     used  SEEK  function also is illustrated.
     1.0   10.00     { Take ten steps, from 1 second to 10 seconds.  This is "T"
       1      -1   { Print every step;  do not bother saving points for plotting
MODELS STAND ALONE
C 3456789012345678901234567890123456789012345678 -- Ruler for following optional 
C                 < File Name>    IPRMDL  KOMPAC    miscellaneous data card :
C MODELS MISC. DATA dc68seek.            9       0
MODEL SAMPLE
C   Of the 6 following variables,  the last is a dummy variable that is used for
C   SEEK function erasure only.  The returned funct value must be put somewhere.
VAR X,Y,Z, X37, X49, XDUM
FUNCTION seek FOREIGN seek {ixarg: 2}
EXEC
C    Before using function # 82,  let's use a second function, the old # 37:
  IF t<4.0 THEN   -- Execute only at the first 3 time steps
C          About the following  WRITE(  statement,  note lack of indentation and
C          use of apostrophes (not quotation marks)  to delimit the text.  These
C          two details are required to maintain any lower case within the text :
write( '    Use SEEK model # 37 for T < 4 only.' )
  X37:=SEEK(T, 37 )
  ENDIF           -- Terminate 4-line IF block  
  IF t=6.0 THEN   -- execute only on the 6th time step
write( '    Erase SEEK function # 37 at T =6.' )
  XDUM:=SEEK(T, -37 )
  ENDIF           -- Terminate 4-line IF block 
C    The following 3 lines use SEEK function # 82,  which is the same as # 37:
X  :=SEEK(T, 82 )
Y:=SQRT(X)*0.3
Z:=SEEK(SQRT(X)*0.3, 82 )
C    The following nested IF loop results in action only for steps 8 and 9:
  IF t>7.0 THEN   -- Execute only for time steps 8 or later
  IF t<10.0 THEN   -- Execute only for time steps 9 or earlier
write( '    Use SEEK model # 49 for T = 7 or 8 only.' )
  X49:=SEEK(T, 49 )
  ENDIF           -- Terminate inner 4-line IF block
  ENDIF           -- Terminate outer 6-line IF block
C    The following IF loop disconnects SEEK function # 49 on step 10:
  IF t=10.0 THEN   -- execute only on the 10th time step
write( '    Erase SEEK function # 49 at T =10.' )
  XDUM:=SEEK(T, -49 )
  ENDIF           -- Terminate 4-line IF block 
C    In order that all data be seen,  let's document the content of DC68SE82.DAT
C    which contains the following formatted lines.  Note the 4 consecutive
C    Y values 142 (points 3 through 6 inclusive).  Thus points 4 and 5 will be
C    omitted.  ATP will store only 7 (X,Y) pairs.  Note also the software EOF
C    (value X = -9999.) which allows arbitrary comment information to follow.
C 0, 69
C 2, 50
C 4, 142
C 4.3, 142,
C 4.7, 142,
C 5, 142
C 6, 188
C 8, 73        
C 10, 100
C -9999., 0,
ENDEXEC
ENDMODEL
USE SAMPLE AS SAMPLE
ENDUSE
RECORD
SAMPLE.X AS X
SAMPLE.Y AS Y
SAMPLE.Z AS Z
SAMPLE.X37 AS X37
SAMPLE.X49 AS X49
ENDRECORD
ENDMODELS
C    Use SEEK model # 37 for T < 4 only.
C    MODELS "SEEK" function # 37 involves 7 pairs of (X,Y).  These are stored in List-20 SCONST cells 3 through 16.
C    3 or more consecutive, identical Y has allowed the omission of     2  incoming (X,Y) points.  The next line reflects this fact.
C    MODELS "SEEK" function # 82 involves 7 pairs of (X,Y).  These are stored in List-20 SCONST cells 19 through 32.
C Blank card ending requests for output variables.  |BLANK  card ending output requests (none) that ATP added to this  STAND ALONE
C Column headings for the  5   EMTP output variables follow.  These are divided among the 5 possible classes as follows ....
C  Next   5     output variables belong to MODELS (with "MODELS" an internally-added upper name of pair).
C  Step      Time       MODELS      MODELS      MODELS      MODELS      MODELS
C                       X           Y           Z           X37         X49
C     0       0.0          69.  2.49198716  72.6314093         69.  88888.8889
C   Use SEEK model # 37 for T < 4 only.
C     1       1.0         59.5  2.31408729  64.4480155        59.5  88888.8889
C   Use SEEK model # 37 for T < 4 only.
C     2       2.0          50.  2.12132034  55.5807358         50.  88888.8889
C   Use SEEK model # 37 for T < 4 only.
C     3       3.0          96.  2.93938769  93.2118338         96.  88888.8889
C     4       4.0         142.  3.57491259  122.445979         96.  88888.8889
C     5       5.0         142.  3.57491259  122.445979         96.  88888.8889
C   Erase SEEK function # 37 at T =6.
C MODELS erasure of SEEK function # 37  at time  T = 6.00000E+00 sec has gained   16  cells of List 20.  IFSEM =   16
C     6       6.0         188.  4.11339276        142.         96.  88888.8889
C     7       7.0        130.5   3.4270979  115.646503         96.  88888.8889
C   Use SEEK model # 49 for T = 7 or 8 only.
C   MODELS "SEEK" function # 49 involves 4 pairs of (X,Y).  These are stored in List-20 SCONST cells 19 through 26.
C     8       8.0          73.  2.56320112  75.9072517         96.         73.
C   Use SEEK model # 49 for T = 7 or 8 only.
C     9       9.0         86.5  2.79016129  86.3474191         96.        86.5
C   Erase SEEK function # 49 at T =10.
C MODELS erasure of SEEK function # 49  at time  T = 1.00000E+01 sec has gained   10  cells of List 20.  IFSEM =   16
C    10       10.         100.         3.0         96.         96.        86.5
BLANK CARD ENDING PLOT CARDS
BEGIN NEW DATA CASE
C     16th of 16 subcases  was  added  29 April 2011  to  illustrate  high-order
C     comparison between  EGAUSS  (see subcase 12)  and LINSOL (see subcase 13).
C     The testing is ordered by using EGAUSS for order two (two equations in two
C     unknows).  Normally the data would consist of 4 cells for the 2x2 matrix 
C     [A] followed by 2 cells for the right hand side.  Instead,  cell 1 of the
C     6 must be the special reserved key 23456.  This is to be followed by five
C     parameters that dedicated code uses to generate the equation set.  Matrix
C     order N  is in cell 2.  The matrix will be symmetrical,  so transposition
C     (a difference between data of EGAUSS and data of LINSOL)  is not an issue.
C     An offidiagonal in cell (I,J) will be given value  -COEFF / SEPAR ** POWER
C     where  SEPAR = IABS ( I - J )  is the distance from the diagonal.  Each of
C     the diagonals  A(I,I)  will be the negative of the sum of all offidagonals
C     plus constant  DIAG.  The right hand side  B(I)  will equal  RHS.  So this
C     is the matrix of a set of resistors connecting each node with every other
C     node.  Ground is a node,  with  DIAG  being the admittance to ground.  RHS
C     is the injected current at each node.
C NEW LIST SIZES
C        0       0       0       0       0       0       0       0       0       0
C        0       0       0       0       0       0       0       0       0  100000
C        0       0       0       0       0   10000       0  126000       0       0
C   240000     742
C   About the preceding program dimensions, two are critical for high order use.
C   These are List 20 and  List 28.  Both are set to the limits of  LISTSIZE.BPA
C   A third,  List 26, must be twice the order of the equation set being solved.
C   In order to simplify things,  the limiting value for this too has been used.
C   If order 31 below is changed to order 310  (this is  arg[2]),  uncomment the
C   preceding five  NEW LIST SIZES  data cards by removing  "C "  from the left.
C   Expect execution to take about 1000 times as long, of course.  Note that 310
C   is close to the limit for List 20 equal to 100K since 3.1**2 is close to 10.
     0.2     1.0   { Take 5 time steps, changing the matrix at each new step
       1      -1   { Suppress accumulation of plot points,  which are unused
MODELS STAND ALONE
MODEL solve
FUNCTION egauss FOREIGN egauss {ixarg:100}
 VAR arg[1..6]
 VAR x[1..3]
EXEC
if t=0 then
 arg[1]:=23456.0 -- Reserved key to request test of LINSOL and EGAUSS
C arg[2]:=310.0  -- Matrix order of the equations being solved.  This is N
  arg[2]:= 31.0  -- Matrix order of the equations being solved.  This is N
 arg[3]:=2.0     -- Exponent to apply to separation.  This is POWER
 arg[4]:=1.0     -- Coefficient of matrix elements.   This is COEFF
 arg[5]:=1.0     -- Right hand side of each equation. This is RHS
 arg[6]:=31.0    -- Constant term of each diagonal.   This is DIAG
endif
x[1..2]:=egauss(2,arg[1..6]) -- The 1st arg, 2, is required for this test
ENDEXEC
ENDMODEL
USE solve AS solve
ENDUSE
ENDMODELS
BLANK card ending plot requests
BEGIN NEW DATA CASE
BLANK