Chapter 6 Computational Aspects

CHAPTER 6

Computational Aspects

6.1. Introduction and Calculation of Fk Tables This chapter comprises a detailed description of the dynamic-programming computational algorithm solutions for obtaining optimal xenon shutdown programs. The digital computer code for the solution of problems (a) and (b), called DYNPROG, is written in ALTAC 111 (FORTRAN 11) algebraic programming language. A detailed listing of the computer instructions in the above language, a glossary of symbols used, and a sample printout are given in the ensuing sections. This computer program (code) contains approximately 300 ALTAC I11 instructions. For a twenty-stage (N=20) problem, with the 20 Fk tables, each consisting of two-dimensional 30 x 30 matrices, and u* =O or M , the average computer running time per problem is about 5.5 minutes on the Philco TRANSAC 2000, Model 210, or about 2 minutes on the Model 212. The fast access memory of these computers contains 16,384 48-bit words. Because of the relative paucity of fast memory storage with respect to this class of problems, using this class of computers, the Fk tables are successively stored on magnetic tape after computation. To compute the Fk table, the Fk-l table must be searched to find the optimal F k - l and uk* values, as will be explained later. This necessitates recalling the Fk- table from the magnetic tape back into the fast (core) memory, which operation consumes a relatively large amount of computer time. The point is that most of the computer running time is consumed by reading from and writing on the Fk table magnetic tape. 81

88

16

Computational Aspects

If the core memory is large enough to accommodate all Fk tables, k=O, 1, 2, ..., N , simultaneously, the problem running time can be reduced by at least an order of magnitude. However, if the core memo-

ry is truly large, e.g., of the order of 250,000 words or more, higher core-memory priority should be given to expanding the Fk tables for greater accuracy. Double interpolation, used presently in DYNPROG, is time-consuming and less accurate than a finer Fk table grid mesh would be [18]. Using (5.46), (5.47), and (5.48), which comprise the dynamic-programming computational algorithm, the computation proceeds “backward” in time through each stage of the multistage control decision process by which the xenon optimal shutdown problem is construed. Assume for illustrative purposes that the xenon override constraint is initially not invoked. It will be incorporated easily later in the discussion. The computation begins at the possible termination set of states ( x T ,y T ) = ( x O ,yo). With zero stages remaining in the control phase, a set of ( x o , y o ) values are chosen, where O < X ~ < Xand ~ ~ ~ ~ O
or

(6.1)

The Fo table is shown in Table 6.1. Using problem (a) to illustrate,

=

min @(xo

Obu,
+ ko(ul)d,

y o +)io(ul)d),

(6.2)

so that the Fl table is manufactured by first choosing a particular pair of (xo,y o ) values. For arbitrary u l , the corresponding transformed state is xl=xo+ko(ul)d and y1=yo+)io(ul)d, where f o ( u l ) and Jjo(ul) are obtained from the xenon-iodine state equations (3.12) and (3.13), respectively. For given valuest of u l , such that 0 < u1 < M , and t Since the optimal control is bang-bang, as shown earlier, u1 = 0, M are the only allowed u values from which to choose in the minimum Fk search computation. The DYNPROG code has the incremental-u-choice option available for future work, when the optimal control may not be bang-bang. Exercising this option even for bang-bang control serves as a good code check to see if the discontinuous nature of the control is actually evidenced.

6.11

89

Introduction and Cahdation Of Fk Tables

TABLE 6.1. Fo(x,,yT)

YT

XI0

X20

x30

x40

x50

X6O

Y 10

Fil

Ff

FA3

Fi4

Fi5

F’,B

~

u;, = 0

YZO

u;, = 0

Fil

F2,Z u;, = 0

u;, = 0

Y 30

F;1 u;, = 0

FY

uj,

=0

Ui3 =

0

FF l4i3 = 0

F3,3 Uj, = 0

=0

Fi4 Ui4 = 0

FS,4

uj,

=0

u;, = 0

F:’ u;, = 0

F3,5 Uj, = 0

u;, = 0

Fie u;, = 0

Fie u;, = 0

corresponding ( x l ,y,), a set F, ( x l ,yl) is computed, using double interpolation from the Fo table. Then the minimum Fo is chosen to give

F l ( X 1 , Y l ) = min FO(X0 OBu,BM

+ ~ o ( u , ) A , Y o+ 3 0 ( u , ) A ) = FO(U,*).

(6.3) The (xl*,yl*) are the locations in the F, table where F l ( x l * , yl*) and u l * obtained above are written. The Fl table is shown in Table 6.2. The Fk tables k = 2, 3, ..., N are successively computed in the same way and stored contiguously on magnetic tape, so that the Fk-l table can be read into the core memory to be searched to compute the Fk table. After all the Fo, F,, ..., FN tables have been computed and stored on tape, the extremal (xk, yk), k = 0, 1, ..., N , in the phase space, and also the optimal flux shutdown sequence {uN*}is obtained by numerically integrating the equations of state (3.12) and (3.13). This is done in the forward direction, and initiated by using the equilibrium conditions x = y = u = l . The “code obtains” the uN*, since N control decisions remain at the outset, corresponding to the x N = Y N = 1 .entry in the FN table. Using uN*,the state equations are integrated to a time A = T / N to obtain ~ ~ - ~ If ~A isy small ~ - enough, ~ . xN-,=xN+iN(uN*)A, y N - , = y , + j N ( u N * ) A . The F N - , table is now entered at XN-l,YN-I, again using double interpolation to obtain the corresponding uN- l*.

90

16

ComputationalAspects

TABLE 6.2.

Fl(x1, y l ) X

Y

Xll

x21

x31

x41

x51

X6'

Yl'

F;'

Fi2

Fi3

Fi4

Fi5

Fie

u11

u12

'13

u14

u15

u1.3

F2,'

F2,2

F2,3

F2,4

u2 1

u22

u23

u24

'25

u2.3

F;'

F?

F:3

F:4

F:'

F:.3

'31

u32

u33

u34

u35

;6'

Y2'

Y31

This procedure is repeated until all the F, tables are exhausted, thus generating the extremal locus, the sequence {xN,yN}.

6.2. The Xenon Override Constraint As mentioned earlier, the preceding explanation of the DYNPROG code did not include the xenon override constraint x < x,. One method by which this constraint can be aecomrnodated into the code is to simply constrain all x, < x, using an appropriate computer instruction. x=x, is a horizontal straight line in the xenon-iodine phase space cutting through the unconstrained extremals and denying them the region of phase space beyond. Evidently this reduces the search space when computing the Fk tables, which reduces the computing-time and fast-memory-storage requirements. The memory storage thereby freed is available for other purposes, such as reducing the mesh size of the (xk, Yk) grid* The second way in which the xenon override constraint can be introduced is through the addition of an artificial cost to the criterion functional, as mentioned earlier. In the DYNPROG code the artifice used is the addition of 10(x/x,)20 to @ or Y. This cost becomes tremendous when x is very close to x,, and rapidly dwindles to a negligible amount for x
6.31

DYNPROG and COAST Input-Data Format

91

extremal never crosses the x = x c line, but on the contrary is diverted from x = x c in a series of sawtooth-shaped segments in the phase space, as discussed later. These segments will be seen to correspond to a flux pulse train which “burns out” sufficient xenon so that the optimal extremal can be rejoined at a later point in the phase space. As mentioned in Section 5.6, another mode of suboptimal shutdown is to make a portion of the extremal coincident with the line x = x c . This, however, requires a relaxation of the upper bound M of the flux constraint, as will be discussed later as well.

6.3. DYNPROG and COAST Input-Data Format The input-data parameters which are required for a given computerproblem run of DYNPROG are punched on a single hollerith card (punched card) according to the format to be described. The card-totape converter enters this data card together with the DYNPROG code stack of punched hollerith cards, to be read onto the computer input tape for processing prior to the problem run. The input data consists of the following parameters. INDEX (1 or 2) determines whether (1) problem (a) or (2) problem (b) is to be solved for the particular computer run number of stages, N = T/A, or uN* decision intervals NMAX length of the control duration (units of iodine mean T life 9.58 hours) ro =ayl,/Al, equilibrium power, or flux, parameter RO post-shutdown time, at which the xenon concenTO tration is minimized; for problem (b) only flux-constraint upper bound, M. M = l is rated UMAX power flux increments used in search for uk*. DU = UMAX DU corresponds to bang-bang control Ax and/or Ay; size of Fk-table grid mesh DXY determines whether or not Fk tables are printed out SKIP (0, or 1) on off-line printer following the computerproblem run MAX I J dimension of Fk tables, < 100 x 100

92

Computational Aspects

CON

16

conversion factor for xenon concentration level to poison reactivity, in dollars xenon override constraint; units of normalized xenon concentration; e.g., XMAX= 10 means an override reactivity constraint corresponding to 10 times the equilibrium poison reactivity

XMAX

The above parameters are tabulated in Table 6.3 in terms of their location on a hollerith card. TABLE6.3. DYNPROG INPUT Card column

Number type

1 2-4 5-9 10-14

I1 I3 F5.0

F5.0

15-19 20-24 25-29 30-34 35-39

F5.0 F5.0

40-43 44-49 50-57

I4 F6.3 F8.1

F5.0 F5.0 F5.0

Parameter

Description

1 for PHI, 2 for PSI problems number of stages (usually 20) end of control phase (units of 11-1) ro parameter corresponding to operating power prior to shutdown TO xmin time, in PSI problems only UMAX upper bound on u* (usually 1) DU u* increment (usually UMAX) DXY size of Fk-table grid mesh SKIP F k table dump bypass (0 prints F k tables, 1 does not) MAX I J size of F k table array poison reactivity conversion factor to dollars CON XMAX xenon override constraint INDEX NMAX T RO

The COAST program, also written in ALTAC I11 language, computes families of solutions of the xenon and iodine differential equations (3.16), which describe the conditions following shutdown where u* ~ 0COAST . is used, with input initial conditions from DYNPROG at shutdown, to continue the extremal in its post-shutdown “coasting” phase. It also generates the data from which Fig. 3.2 is drawn. That is, COAST can provide the data for the xenon and iodine following immediate shutdown as well, where x(O)=y(O)=l and u*=O. The input data to COAST consists of the following parameters.

Appendix to Chapter 6

6.A]

xo

93

xenon-concentration shutdown level, or xenon initial condition for immediate shutdown iodine-concentration shutdown level, or iodine initial condition for immediate shutdown time increment ro = a p o / l , power or flux parameter number of time increments per run xenon poison reactivity conversion factor to dollars

YO DEL RO IMAX CON

The above parameters are tabulated in Table 6.4 as they would appear on a hollerith (punched) data card. TABLE 6.4. COAST INPUT Card column

Number type

1-9 10-18 19-27 28-36 37-41 4247

F9.4 F9.4 F9.4 F9.4

I5

F6.3

Parameter XO YO DEL RO IMAX CON

Description xenon concentration initial condition iodine concentration initial condition time increment (usually 0.1) equilibrium flux parameter number of time increments reactivity conversion factor

6.A. Appendix to Chapter 6

6.A.I. DYNPROG Glossary A B C CON COY DEL=d =T/NMAX DRU DU=du DXY 8 DYNPROG, ~ 2 3 , ~ DY1

=

wW/(l-w)

yl(w

= (1 - W ) / Y l ( W

+ ro)

+ ro)

= 1/1 - w factor for converting x to reactivity in dollars = --(y~/(l - w)[(exp - A ) - (exp - ~ 4 1 ) time increment used throughout code

= rodu u increment in u* search loop x and y interval which fixes Fk matrix grid interval size

code name; ALTAC 3; compiler TAC; binary deck and code check ; execute with conditional post-mortem dumps = du(1 - exp - d)

Computational Aspects

94

EDEL EDEO EX EU FKJ) FNJ) 4700-3 1000F, 0-51OOc, 30000-35000H FU

= exp -

= exp -

A

wA

+

exp - A(w rodu) = exp - (du)(rod) generic symbol for Fk(x,y) interpolated F k corresponding to jth column floating-point dump, command format, =

hollerith (input-output) dump interpolated F k with next value of u for comparison with FR(J) in U* search loop F k table entry on solution print out format = yl(w ro)

FWRDF GAM 1 GAM 2 = ya(w ro) index of X(rows) of F matrix, 1 < I < MAXIJ - 1 I INDEX 1,2 1 is @(x,y), problem (a), 2 is Y ( x , y ) , problem (b) tape-read check for correct I IPRIME IXNMTROTO INDEX, NMAX, T, RO, TO, UMAX, DU, DXY, SKIP, UMDUDXY MAX, IJ, CON; SKP index of y(columns) of F k matrix, 1 < J < MAXIJ - 1 J format for punching data cards MIJCON max number of rows(1) and columns(J) in F k table matrix S 1MJ MAXIJ instruction to mount extra tape used to store F tables MOUNT TAPE maximum number of stages of F k NMAX tape-positioning parameter NMA =NMAX 1 tape-read check for correct N NPRIME tape backspace-positioning parameter NRBACK RAYDMP special dump =w+ro RC RO = ro = o ( p o / L l ; equilibrium power parameter =w rodu RU computer time-sharing card /" = x,y/DXY; linear interpolation interval THX,THY = 1 -THX, THY THXC,THYC u1 = yl(w ro)(y - u ) / ( w rodu - 1) = (w ro)u/(w roAu) u2 UMAX urnax M = - y1/l - w uo1 = u*(Jdu),1 < J < MAXIJ 1 US(J)

+ +

+

+

+

us0

USTAR

u* U*

E O

+

+

+

+

6.A]

Appendix to Chtzpter 6

XMAX XY XY MAX

95

xe, xenon override constraint

x or y maximum value of x or y in

Fk

matrix

6.A.2. Sample Printout of DYNPROG NMAX = 20 UMAX = 2.00 D U = 2.00 XYMAX = 29.0 DEL = 0.05 B = 0.01 EDEL = 0.95 DYl = 0.1 CAM2 = 1.1 EU = 0.14

MAXIJ = 30 GAMl = 19.7 D R U = 40.00

Objective Function PSI

RO

= 20

T = 1.000

N

X

20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0

1.00000 1.90647 2.73477 3.48991 4.17663 4.79938 5.36238 5.86957 6.32469 6.73125 7.09256 7.41173 7.69168 7.93517 8.14478 1.33704 1.82637 2.27211 2.67710 0.62092 0.38246

DEL

= 0.05000

USTAR

Y 1.m 0.95123 0.90484 0.8607 1 0.81873 0.77880 0.74082 0.70469 0.67032 0.63763 0.60653 0.57695 0.54881 0.52205 0.49659 0.56991 0.54211 0.51567 0.49052 0.56414 0.63417

0.

0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.

0. 0. 2.000 0. 0. 0. 2.000 2.000 0.

TO

=

1.00

YMAX

=

FN

CON

0.3 13142 0.300607 0.264368 0.204034 0.1 19307 0.321586 0.193685 0.372604 0.199694 0.353619 0.131254 0.259472 0.381486 0.497165 0.175914 0.23 1 196 0.480335 0.224100 0.448630 0.386202 0.634170

8.1250 15.4901 22.2200 28.3555 33.9351 38.9950 43.5693 47.6903 51.3881 54.6914 57.6271 60.2203 62.4949 64.4732 66.1764 10.8635 14.8393 18.4609 21.7514 5.0450 3.1075

6.A.3. D YNPROG and COAST AL TAC III Instruction Listing 20

40 ABC I

8.15

HLT 14-087-8 ON MAGTAPE, LIBS ALTAC3 DYNPROGS IOWS IDENTIFYF,

4

96 C

ComputationalAspects

C FUNCTION AFMIF(1,DXY) COMPUTES VARIABLE FROM SUBSCRIPT C AFMIF (l,DXY)=(l-1.) *DXY C SUBROUTINE IFMA (XY,TH,IJ) COMPUTES SUBSCRIPT FROM X AND Y VARIABLE SUBROUTINE IFMA(XY,TH.IJ) COMMON MAXIJ,XYMAX,DXY IF(XY) 10,10,20 10 IJ=I IH=XY/DXY RETURN 20 IF(XY -XYMAX)40,30,30 30 IJ=MAXIJ-1 IH=(XY -XYMAX+DXY)/DXY RETURN 40 IJ=XY/DXY+l. T H = (XY -AFMIF(IJ.DXY))/DXY RETURN END C MAKES HOLLER CODE FOR PHI A N D OR PSI C SUBROUTINE TYPE (P) DIMENSION P (1) PRINT 902,P 902 FORMAT ( I H0,20x,l9HOBJECTIVE FUNCTION A3) RETURN END C C C MAIN PROGRAM BEGINS C COMMON MAXIJ,XYMAX,DXY DIMENSION F( 100,100),FR(l00). US(100) CONTINUE C READ 900 S GETS BY CASE CARD C COLS TYPE PARAM DESCRIPTION C 1 I1 INDEX 1 FOR PHI OR 2 FOR PSI PROBS C 2-4 13 NMAX NUMBER O F STAGES (USUALLY 20) C 5-9 FS.0 T END OF CONTROL PHASE (UNITS O F INVERSE LAMBDA ONE) C 10-14 FS.0 RO PHI NAUGHT PARAMETER C 15-19 FS.0 TO XMlN TIME (IN PSI PROBS ONLY) C 20-24 F5.0 UMAX UPPER BOUND ON U (USUALLY I.) C 25-29 FS.0 DU U INCREMENT (USUALLY I . ) C 30-34 FS.0 DXY COORDINATE INTERPOLATION INCREMENT (USUALLY 0.1) FN TABLE DUMP BYPASS(0. DOESNT AND 1. DOES C 35-39 FS.0 SKIP BYPASS) C 40-43 14 MAXlJ SIZE O F FN TABLE (NOT GREATER THAN 100 x 100) C 44-49 F6.3 CON K MULTIPLIER FOR X C 50-57 F8.1 XMAX XMAX XENON OVERRIDE CONSTRAINT 1 READ 900, INDEX,NMAX,T,RO,TO,UMAX,DU,DXY,SKIP,MAXIJ,CON,XMAX REWIND 6 900 FORMAT (II,I3,7FS.O. I4,F6.3,F8.1) IF(1NDEX) 2,Z.S

6.A1

Appendix to Chapter 6

2 STOP C COMPUTATION O F CONSTANTS DEL = T/NMAX B = 8.*59./(56*(21. +29.*RO)) XYMAX=(MAXII- l.)*DXY RC=21./29. +RO GAMI=56./59.*RC GAMZ=3./59.*RC EDEL=EXPF(-DEL) DY 1=DU*(l. -EDEL) DRU=RO*DU EU=EDEL**DRU EDEO = EDEL**(21./29.) UOl =GAMl*(-29./8.) COY = UO 1*(EDEL - EDEO) GO TO (100,200)JNDEX C PHICOMPUTATION 100 A =(21./29.)r*(21./8.)~56./59.~(21./29. RO) C=29./8. 101 DO105 I = l , MAXlJ X =AFMIF(I.DXY) 102 DO105 J=l,MAXIJ Y =AFMIF(J.DXY) IF(Y) 108,108.109 108 F(I,J) = 10.E10 GO TO I05 109 F(I,J)=A*Y*(l +B*X/Y)**C 105 CONTINUE 8 I101,J102 GO TO 250 c PSI COMPUTATION 200 A T 0 = EXPF(21./29.*(T -TO)) BTO=(ATO-EXPF(T-TO))/B 20I DO205 I=l.MAXIJ X =AFMIF(1,DXY) ATOX =ATO*X 202 DO205 J = 1.MAXIJ Y =AFMIF(J,DXY) F(I,l)=ATOX+BTO*Y 205 CONTINUE 8 IZOl,JZ02 C STORES F AND USTAR TABLES FOR N=O 250 DO251 J = l , MAXIJ 251 US(J) =O.

+

N=O

DO260 I=l,MAXIJ IF(SKIP)255,255.260

C C C PRINT FO AND US0 C

C

255

905 906

PRINT 905,N,I,(F(I,J).J= 1,MAXIJ) FORMAT (3HON=I2,SX2HI=I3/(lH 10Fl2.5)) PRINT 907 PRINT 906,(US(J),J= 1,MAXII) FORMAT(1H 10F12.5)

97

98

Computational Aspects

I6

C COMPUTATION O F F(N,I,J) AND US(N,I,J) WRITE TAPE 6,N,I,(F(I,J),J= l.MAXIJ),(US(J),J= 1,MAXIJ) 260 300 DO 499 N=l,NMAX 301 DO399 I=l.MAXIJ X=AFMIF (1,DXY) 302 DO365 J=I,MAXIJ Y =AFMIF(J,DXY) C COMPUTATION O F F(1.J) FOR U=O FOLLOWS XI =EXEO*X+COY*Y Y1 =Y*EDEL CALLIFMA(XI,THX,Il) CALLIFMA(Y l,THY,Jl) THXC = I. -THX THYC= 1. -THY FR(J) =THXC*THYC*F(II ,J1) +THXC*THY*F(Il,Jl 1)+ lTHXoTHYCoF(I1 l,JI)+THXoTHY*F(Il+ l,J1+ I ) + IO.r(X/XMAX)roZO US(I)=O. C INITIALIZATION FOR U LOOP U=DU RU =21./29. +DRU EX= EDEO*EU c u LOOP 310 U1= GAM l*(Y - U)/(RU - 1 .) U2= RC*U/RU X1 =(X - U l - U2)oEX +Ul*EDEL+ UZ Y1 =Y 1 +DY I 350 CALLIFMA(X I ,THX,Il) CALLIFMA(Yl.THY,Jl) THXC= I. -THX THYC = 1. -THY

+

+

+

FU=THXCoTHYCoF(II,Jl)+THXC*THYoF(II.JI 1)+

lTHX*THYCoF(Il +l.JI)+THX~THYoF(II +1,J1 + I ) 704 IF(FU- FR(J))355,360,360 355 FR(J)=FU US(J)=U C INCREMENT U.RU, EX 360 U = U + D U RU =RU +DRU EX=EXoEU IF( U - U MAX) 3 10.310,365 C END U LOOP 365 CONTINUE S J FROM 302 WRITE TAPE 6,N,I,(FR(J),J= I,MAXIJ),(US(J),J= 1,MAXIJ) IF(SKIP)370,370,399 C C C PRINT FN AND USN C 370 907 399

PRINT 905,N,I(FR(J),J= 1,MAXIJ) PRINT 907 PRINT 906,(US(J).J= 1,MAXIJ) FORMAT(1H) CONTINUE S I FROM 301 DO401 M=l,MAXIJ

6.A]

Appendix to Chopter 6

99

BACKSPACE 6 DO415 I=l,MAXIJ READ TAPE 6,NPRIME,IPRIME,(F(I,J),J= 1,MAXIJ) IF(N -NPRIME) 410,403,410 403 IF( I-IPRIME) 410.415, 410 4 10 PRINT 9 15,N,NPRIME,I,IPRIME 915 FORMAT (18HO TAPE NOT SYNCHED,415) STOP 415 CONTINUE S I FROM 402 499 CONTINUE S N FROM 250 C FORMATTING USTAR AND FN TABLE PLUS COMPUTATION OF BOTH I ,GAMZ,EDEL, PRINT 921,NMAX,UMAX.DU,MAXIJ.XYMAX,DEL,B,RC,GAM DY 1,DRUI ,EU,EDEO,UOl.COY 921 FORMAT(6HlNMAX=l2,4X,5HUMAX=F5.2,4X,3HDU=F5.2.4X, 6HMAXIJ =I4,4X, I6HXYMAX = F5.1,4X.4HDEL = F4.1.4X,ZH B = F6.2,4X. 3HRC=75.1/ ~HOGAM~SF~.~I,~X,~HGAM~=F~.I.~X,~HEDEL=F~.~,~X, 401

402

4HDYl=F4.1,4X,4HDRU=F7.2,4X,33HEU=F7.2,4X,5HEDEO=F6.3,4X, 4HUOt=F5.2.4X,4HCOY=F6.3)

PRINT 920,SYSDATE FORMAT (IHO,ZOX.ZOH OPTIMAL USTAR POLICY 2XA8) GO TO (501,502)INDEX CALL TYPE (3HPHI) 50 I PRINT 925,RO.T.DEL.XMAX FORMAT (6HORO = ,FS.O,5X4HT = F5.3,5X6HDEL=F6.5.5X7HXMAX=F5.2) 925 GO TO 503 502 CALL TYPE (3HPSI) PRINT 930,RO,T,DEL,TO,XMAX FORMAT (6HORO =,F5,0,5X,4HT = F5.3,5X6HDEL= F7.5.5XSHTO=F5.2, 930 15X7HXMAX = F5.2) GO TO 503 503 PRINT 935 935 FORMAT(/lH,7X,I HN, 15X.lHX.15X. IHY,I2X,SHUSTAR, 14X,2HFN//) C FIRST LINE O F USTAR AND FNMAX TABLE X=l. Y=l. I=OS TAPE 6 NOW AT END O F F(NMAX)=ROW 1 O F F(NMAX4-I) NMA=NMAX+ 1 510 D O 525 LINE= I.NMA N=NMAX -LINE+ 1 NRBACK =MAXIJ+ I + I CALL IFMA(X.THX.1) CALL IFMA(Y,THY,J) THXC = 1. -THX THYC = 1. -THY IF(THX)544,545,545 544 IOTA = 1 GO TO 560 545 IF (THX-S) 549,549,546 546 I F (THY 1.)547,547.548 547 IOTA= I+ 1 GO TO 560 548 IOTA =MAXIJ GO TO 560 549 IOTA = 1 920

-

GO TO 560

100

ComputationalAspects

IF(THY)561,562,562 JOTA= I G O TO 588 562 IF(THY - .5)566.566,563 563 IF(THY - 1.)564.564,565 564 J O T A = J + l G O TO 588 565 JOTA= MAXIJ G O TO 588 566 JOTA=J G O T O 588 588 NRBACK=NRBACK-I DO520 M = L.NRBACK 520 BACKSPACE 6 C READS ROW I O F FN TABLE READ TAPE 6, NPRIME,IPRIME.(FR(L),L = I,MAXIJ),(US(L),L= I,MAXIJ) C TAPE CHECK SUBROUTINE 600 IF(NPRIME-N)601.605.603 C NPRIME LESS THAN N 601 KR=MAXIJ (N-NPR1ME)-I DO602 K=L.KR 602 READ TAPE 6 G O T O 600 C NPRIME GREATER THAN N 603 KR=MAXIJ (NPRIME-N)+l DO604 K = I,KR 604 BACKSPACE 6 G O TO 600 C IPRIME LESS THAN I 605 IF(IPR1ME-I) 606,515,608 606 KR=I-IPRIME-I DO607 K = I,KR 607 READ TAPE 6 G O T O 605 C IPRIME GREATER THAN I 608 KR=IPRIME-1+1 D O 609 K=I,KR 609 BACKSPACE 6 G O TO 605 C TAPE 6 NOW AT RECORD I + l O F USTAR A N D F N TABLE 515 FOO=FR(J) FIO=FR(J+I) IF(1OTA - 1)517,516,517 516 U=US(JOTA) 517 READ TAPE 6,NPRIME,IPRIME,(FR(L),L= I,MAXIJ),(US(L),L= 1,MAXIJ) BACKSPACE 6 C TAPE 6 RESTORED TO RECORD I + l FOI = FR(J) FII = F R ( J + I ) IF(1OTA -1)518,519,518 51 8 U = US(J0TA) 519 FWRDF=THXCoTHYCoFOO+THXC.THYIFOI+THXoTHYCoFIO+THXoTHY.FL1 XC=X*CON PRINT 940,N,X,Y,U,FWRDF,XC 940 FORMAT( IH,7X,I2, IOX,F9.5,8X,F9.5,9X,F5.3,7X,EI 3.6,2X,Fl1.4)

560 561

6.A]

Appendix to Chapter 6

C OPTIMAL TRAJECTORY EXTRAPOLATION RU = 21 ./29. ROIU UI = GAM Ir(Y - U)/(RU - 1 .) U2=RC*U/RU EX = EDEL**RU X=(X - U1 -U2)*EX+UIrEDEL+U2 702 Y=U+(Y-U)rEDEL 525 CONTINUE 8 LINE FROM 510 GO TO I END S SCOAST rF2 3C O35OOOF 035OOOC O35OOOH * N I COAST IDENTIFYFS DIMENSION X( lOOO),Y(lOOO) READ 9008 CASE CARD C FORMAT CARD COLUMN PARAMETER C F9.4 1-9 XOJNITIAL X C F9.4 10-18 Y0,INITIAL Y C F9.4 19-27 DEL , TIME INCREMENT c F9.4 28-36 R0,PHI NAUGHT C I5 37-41 IMAXJNDEX UPPER LIMIT C F6.3 42-47 C0N.K MULTIPLIER FOR X(I) C 1 READ 900,XO,YO,DEL,RO,IMAX,CON 900 FORMAT (4F9.4,15,F6.3) 3 PRINT 905,XO,YO,DEL,RO,IMAX 905 FORMAT (4HIXO=F9.4,5X,4DEL=F9.4,5X,3HRO=F9.4,5X. 1 SHIMAX = 15//IH,3X,IHI,SX,4HX(I),5X,4HY(l)//) A=56.*(21. +29.1R0)/(8.059.) DOllO I=I,IMAX 100 X(I)=(XO+ArYO)rEXPF( -21./29.r(1- I)*DEL) -A*YOrEXPF( -(I- 1)rDEL) 101 Y(I) = YO*EXPF( - (I - I)rDEL) XOC=X(I)*CON GO TO 105 105 PRINT 500,I,X(I),Y(I),XOC

+

500

FORMAT(IH,I5,2F9.4,2X.F11.4) I10 CONTINUES1 GO TO 1

ENDS (DATA END DATA

XO YO DEL RO IMAX CON

101