Index Circular Stirling Stirling program Source code Links

Optimized Stirling engine design quickbasic program
The engine of the future project
The following is a Stirling engine design program, based on a translation to Quick Basic of the isothermal second order calculation program,
written by William R. Martini in Fortran.
I found the Martini program on:
http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19830022057_1983022057.pdf
and translated it into quickbasic.
This is a version that makes more use of the input and output features that modern programming languages offer.
There is an optimizer in this program.
This means that you can chose the program to find a value for an input that predicts the highes efficiency or the highest net power, or a combination of both.
You can even optimize a lot of inputs at the same time.
I even used the optimizer to change the program so that the outputs are very close to the General Motors test results (you can find them on my site).
In the future I wil also compare the outputs with the GPU-3 test results.

How do you download the program? In fact it is already on your computer.
You can select the part between the ===== lines and copy it in a simple editor and save it with a name with the extention .bas
Download a version of Quick Basic or QBasic and load the program. On the end of the page you find the help file.
Copy the text between the ############# lines and save it with the name STRLD.HLP in the same directory.
This gives you the possibility to see what the inputs, outputs and options mean.
You select the item where you want to know more about by moving the -> arrow and then type ?
When you think something does not work like it should, please e-mail me (hvd@haghen.nl).
Currently I'm working on a better version.
You can find the old version of the program here:
stirlingmartiniold.html
For the html-javascript version of the Isothermal second order calculation Program (Written by William R. Martini in Fortran) Click here
The program is in the source code.
When you put the right inputs in the html page, and you click on:
Save current inputs in cookie
you see the inputs in the message box.
You can put those inputs in this basic program.

Problems


For an Alpha type Stirling make the input DD Pist.D.ROD DIAM zero.
In a newer version I will make an extra page to chose for different types of Stirling engines.


Program STRLG33.BAS
'====================================================================

DECLARE SUB makestr (printv!, printv$)
DECLARE SUB makestr (printv!, printv$)
DECLARE SUB existfile (filename$, filecount%, mes$)
DECLARE SUB dosfilename (filename$, mes$)
DECLARE SUB givefilenr (filecount%, filenr%)
DECLARE SUB closefilenr (filecount%, filenr%)

true = 1
false = -1

'Isothermal second order calculation

'Program 150 -10 oct 1979-
'Written by William R. Martini
CLS



'init calculator *********************************************************
'C --- PAG. 334

'sets up arrays (dimensions)
DIM H(13), C(13), P(13), FH(13), FC(13), V(14)
'sets up integers
'INTeger A1%, OG%, ZA%,ZB%,ZZ%,CRT%,TRM%
'Sets up real numbers
'real IC,ID,IH,IP,KA,KB,KG,KM,K3,L1,LB,LD,LE,LI,LR,LX,LY,M,
'ME,MF,MR,MU,MW,MX,M1,M2,M3,NP,NU,LC,LH,L(14),NT,ND
DIM L(14)
'Sets up logical unit numbers. 'CRT' is the logical unit number for
'the input file, and 'LP' is the logical unit for the output
'file
PI = 4 * ATN(1)'3.141593'more exact

P4 = PI / 4
CM = 2.54 'cm/inch   not used

'init optimizer
CONST nivar = 42' number of input variables
CONST novar = 23' number of output vaiables

DIM vi(nivar)
DIM ni$(nivar)

DIM vo(novar)
DIM no$(novar)
DIM eo$(novar)

inifactor# = .1#
tuning# = .9#


skipit = false
do$ = ""

GOSUB 2006 'Init load inputs abbreviationsni$()
GOSUB 2035 'init load output vars names no$()
'op$ = "                       "
'ip$ = "                                          "

'start ---------------
GOSUB 1335 'make or read STIRLSYS.TXT file
CLS
GOSUB 1124 'execute start settings

10 'loop calculate
CLS
IF do$ = "calculate" THEN
  GOSUB 1212 'start optimizer
ELSEIF do$ = "outpscreen" THEN
  GOSUB 2034 'load output values to vo()
  GOSUB 2011 'outputs screen
ELSEIF do$ = "optimize" THEN
  GOSUB 1241 'optimizer screen
ELSE
  GOSUB 1100 'start screen
END IF
GOTO 10 'loop

END

main:

'Program reads in engine dimensions, operating conditions, and
'conversion constants from the input file. Also this is the return
'point after a case has been completed. If there are no more cases to
'run (i.e. an end of file occurs), the program cals exit.


mes$ = ""
IF (FF >= 1) OR (FF <= 0) THEN
  mes$ = "< EROR > Wrong generator fill factor (FF must be between 0 and 1)."
  RETURN
END IF

'#@@@ converting inputs to value per unit
ZH = ZH / N% 'static heat loss
'BF = BF * 4
'BF = 1.448
'BF=3.7
BF = 1.5

'The degree increment is set at 30 degrees
ND = 30
'A correction factor is calculated which increases the accuracy in
'calculating the work integrals with 30 degree increments.
XX = 1 + 5.321E-05 * ND ^ 1.9797
'Temperature change for cold space (DU) and temperature change for hot
'space (DV) are set.
DU = 64
'C --- PAG. 335
DV = 64
'The first thing the program does is to compute a list of engine
'volumes.

'Conversion to Kelvin degrees from input fahrenheit degrees.
TM = (TF + 460) / 1.8
TW = (TY + 460) / 1.8
'conversion to hertz and to npa.
NU = SP / 60
PG = .006894 * PS

IF skipit <> true THEN
PRINT "Heater tube inside temp.="; TF; "fahrenheit="; TM; "Kelvin="; TM - 273.15; "Celcius"
PRINT "Cooling water temperature="; TY; "fahrenheit="; TW; "Kelvin="; TW - 273.15; "Celcius"
PRINT "engine speed="; SP; " RPM="; NU; "RPS"
PRINT "average pressure="; PS; " psia="; PG; "MPa"
END IF

'Determines gas property values from 'OG%' (
'if 'OG%'= 1, the property values for hydrogen are used.
'If 'OG%' = 2, the property values for oxygen are used. (MUST BE HELIUM)
'If 'OG%' = 3. the property values for air are used.)
'Property values for additional gases may be added if desired'
'IF OG% = 1 THEN GOTO 20
'IF OG% = 2 THEN GOTO 21
IF OG% = 3 THEN 'air propperty values
  IF skipit <> true THEN
  PRINT "WORKING GAS= AIR"
  END IF
  KA = -12.6824
  KB = .782
  CP = 1.0752
  CV = .7883
  M1 = 1.8194E-04 'MU=M1+M2*(TR-293)+M3*(PAVG) ,g mass/cm sec, PAVG in MPa
  M2 = 5.36E-07
  M3 = 1.22E-06
  MW = 29
  PR = .9071
ELSEIF OG% = 1 THEN 'hydrogen propperty values
  IF skipit <> true THEN
  PRINT "WORKING GAS= HYDROGEN"
  END IF
  KA = -11.0004
  KB = .813
  CP = 14.62
  CV = 10.49
  M1 = 8.873E-05
  M2 = .0000002
  M3 = 1.18E-07
  MW = 2.02
  PR = .8408
ELSEIF OG% = 2 THEN 'HELIUM propperty values (NOT OXIGEN !!!!)
  IF skipit <> true THEN
  PRINT "WORKING GAS= HELIUM" 'OXYGEN"
  END IF
  KA = -10.1309
  KB = .6335
  CP = 5.2
'C --- PAG. 336
  CV = 3.12
'  M1 = 1.6614E-04' 196.14E-6, see pag. 301
  M1 = 1.9614E-04' , see pag. 301
  M2 = 4.63E-07' 4.64E-7 on page 301
  M3 = -9.3E-08' was -9.3E8 -0.093E-6 pag. 301
  MW = 4
  PR = .8018
ELSE
  PRINT "  _______ #### ERROR, INPUT FAILURE #### ___________"
  PRINT "OG% = wrong, must be 1, 2 or 3 for hydrogen, oxygen, or air."
  INPUT "ENTER for end program"; i$
END IF
'Conversion of cooling water flow to grams/second. Initially cooler
'tube metal temperatue is made the same as the inlet cooling water
'temperature. The total heay transfer areas for all the engines
'coolers and all the engine heaters are calculated.

IF (MS% AND 4) = 0 THEN
  FW = 63.12 * FX
' FX= GPM @ 2000 rpm per cilinder
ELSE
' Has to be:
  FW = 63.12 * FX * SP / 2000
' when the waterflow is proportional to engine speed.
END IF

IF skipit <> true THEN
PRINT "flow of cooling water="; FX; "GPM at 2000 rpm per cilinder="; FW; "g/sec"
END IF

TX = TW

'The multiplication with N% must be a mistake,
' because later in the program AC is multiplied by N% again.

AC = PI * IC * LD * NC%

IF skipit <> true THEN
PRINT "heat transfer area for cooler="; AC; "cm^2"
END IF

'The multiplication with N% must be a mistake,
' because later in the program AH is multiplied by N% again.

AH = PI * IH * LI * NH%

IF (MS% AND 16) = 16 THEN
 'Cylinder head area is added to heat transfer area.
  AH = AH + P4 * DC
END IF

IF skipit <> true THEN
PRINT "heat transfer area of heater="; AH; "cm^2"
GOSUB 1188 'Hit key
END IF

'calculates engine dead volumes and initializes pressures and volumes.
'Initializes for determination of average pressure and maximum and
'minimum volumes.

HD = P4 * IH * IH * LH * NH% + EC * DC ^ 2 * P4

IF skipit <> true THEN
PRINT
PRINT P4 * DC ^ 2 * EC; "= piston end clearance, dead volume, cm^3"
PRINT P4 * IH ^ 2 * LH * NH%; "= heater tube dead volume, cm^3"
PRINT "+ -----------"
PRINT HD; "= total hot dead volume, cm^3"
PRINT
END IF

'CX = P4 * ID * LE * NE%
'          ID   inside diameter of connecting duct, cm
'                LE   length of connecting duct, cm
'                      NE%   number of connecting ducts per cylinder
IF (MS% AND 2) = 0 THEN
  CX = P4 * ID * LE * NE% '
ELSE
  CX = P4 * ID ^ 2 * LE * NE%'  ^2 was forgotten in original
END IF

RD = (1 - FF) * P4 * DR ^ 2 * LR * NR% + PI * DC * G * LB

IF skipit <> true THEN
PRINT (1 - FF) * P4 * DR ^ 2 * LR * NR%; "= regenerator inside dead volume, cm^3"

PRINT PI * DC * G * LB; "= gap in hot cap dead volume, cm^3"
'pi* DiameterEngineCylinder*GapInHotCap*LengthHotCap

PRINT "+ -----------"
PRINT RD; "= total regenerator dead volume, cm^3"
PRINT
END IF

CTDV = P4 * IC ^ 2 * LC * NC%
ECCDV = EC * P4 * (DC ^ 2 - DD ^ 2)
CD = CX + CTDV + ECCDV

IF skipit <> true THEN
PRINT CX; "= cold dead volume outside cooler tubes, in connecting duct cm^3"
PRINT CTDV; "= cooler tube dead volume cm^3"
PRINT ECCDV; "= end clearance cold dead volume cm^3"
'EC  piston end clearance, cm
'DC  diameter engine cylinder, cm
'DD  diameter of piston drive rod, cm
PRINT "+ -----------"
PRINT CD; "= total cold dead volume, cm^3 (incl. connecting duct)"
PRINT
GOSUB 1188 'Hit key
END IF

PM = 0
VX = 0
VN = 1E+30
'initially sets the effective hot space temperature to the hot metal
'temperature and the effective cold space temperature to the cooling
'mater temperature for the first time around. Calculates the log mean
'temperature for the regenerator. Calculates the leakage coefficient
'for 30 dergree increments.
TH = TM
TC = TW
TR = (TM - TW) / LOG(TM / TW)
LX = L1 * ND / (360 * NU)

'Since the thermoconductivity enters the calculation only at the
'regenerator temperature it can be calculated before the main
'iteration loop
KG = EXP(KA + KB * LOG(TR))
'Start of do loop 23 to calculate engine volumes.

'C --- PAG. 337
FOR i = 1 TO 13
'Calculates the hot volume and cold volume for each angle increment for
'crank operated pistons. Since a double acting machine has a piston
'drive rod (DD) and a single acting machine does not, 'DD' is used as
'an indicator of wether the cold volume of the engine is above the
'piston or below it.
  X = 30 * (i - 1) * RA
  J = i
  IF DD <> 0 THEN
    Y = (30 * (i - 1) + AL) * RA
  ELSE
    Y = (30 * (i - 1) - AL) * RA
  END IF
  H(J) = P4 * DC ^ 2 * (RC - SQR(CR ^ 2 - (RC * SIN(X)) ^ 2) + RC * COS(X) + CR) + HD

  IF DD <> 0 THEN
    IF (CR ^ 2 - (RC * SIN(Y)) ^ 2) < 0 THEN
      mes$ = "< ERROR > Crank radius is to large."
      RETURN
    END IF
    C(J) = P4 * (DC ^ 2 - DD ^ 2) * (SQR(CR ^ 2 - (RC * SIN(Y)) ^ 2) - RC * COS(Y) - CR + RC) + CD
  ELSE
    C(J) = P4 * DC ^ 2 * (RC - SQR(CR ^ 2 - (RC * SIN(Y)) ^ 2) + RC * COS(Y) + CR) + CD
  END IF
  'calculates the total gas volume and finds the maximum volume.
  V(J) = H(J) + RD + C(J)
  IF V(J) > VX THEN VX = V(J)
  'finds the minimum volume.
  IF V(J) < VN THEN VN = V(J)
  'calculates the initial gas inventory.
  IF J = 3 THEN L(1) = PG * (H(J) / TH + RD / TR + C(J) / TC)
  'end of loop to calculate engine volumes
NEXT i'continue
' 'ZA%' is set at zero so that the fastest way of arriving at the proper
'effective hot space and cold space temperature will be tried first.
'also a counter, 'ZB%', is set at zero.
ZA% = 0
ZB% = 0
GOSUB 200 'initialization
IF mes$ <> "" AND ZZ% = 0 THEN
  ZH = ZH * N%
 'Because ZH can be input and output, it has to be reset in case of an error.
END IF
RETURN

'initialization
200 A = 0
WHILE (X > .0001) OR (Z > .01)

 PM = 0
LY = 0
'start of loop (to calculate pressures for one engine cycle).
FOR i = 1 TO 13
'C --- PAG. 338
  'calculate pressure
  P(i) = L(i) / (H(i) / TH + RD / TR + C(i) / TC)
  'Calculate gas inventory for next increment due to leakage
  L(i + 1) = L(i) * (1 - LX * (P(i) - PG))
  'Accumulate values, mean pressure and mean gas inventory.
  IF i <> 1 THEN
    PM = PM + P(i)
    LY = LY + L(i)
  END IF
  'end of loop (to calculate pressures for one engine cycle)
NEXT i
'indexes cycle counter, calculates mean pressure, readjusts gas
'inventory taking into account gas leakage.
A = A + 1
PM = PM / 12
IF A >= 3 THEN
  L(1) = L(13)
ELSE
  L(1) = L(13) * PG / PM
END IF
'convergence criteria: pressure from beginning to the end of cycle
'must not change by more than one hundredh of a percent and the mean
'pressure must be within one percent of the desired gas pressure.
'Usually one or two cycles are required to meet this criteria.
X = ABS(P(1) - P(13))
Z = ABS(PM - PG)

IF skipit <> true THEN
PRINT X; "Pressure difference after loop"; A
PRINT Z; "Volume difference after loop"; A
END IF
WEND

IF skipit <> true THEN
PRINT "Pressure and Volume are O.K.";
GOSUB 1188 'Hit key

END IF

'initializing
W1 = 0
PX = P(1)' was 0, this is better
PN = P(1)' was 10000, this is better
MR = LY * ND / 360
'Start of loop (finds the maximum and minimum pressure).
FOR i = 1 TO 13
  IF P(i) > PX THEN PX = P(i)
  IF P(i) < PN THEN PN = P(i)
NEXT i'continue

IF skipit <> true THEN
PRINT PX; "= Maximum pressure"
PRINT PN; "= Minimum pressure"
GOSUB 1188 'Hit key
END IF

'C --- PAG. 339
'Start of loop (finds the work per cycle by integrating the
'pressure volume loop).
FOR i = 1 TO 12
  W1 = W1 + (P(i) + P(i + 1)) * (V(i + 1) - V(i)) / 2
NEXT i'continue
'Basic power for the whole engine is calculated from the integrated
'power using the correction factor XX which compensates for the
'truncation error of using only a small number of points to integrate.
BP = NU * XX * W1

'Initializing
HX = 0
CY = 0
HN = 1
CN = 1
'Calculates an array giving the fraction of the total gas inventory in
'the hot space and in the cold space for each point during the cycle.

FOR i = 1 TO 13
  FH(i) = P(i) * H(i) / (MR * TH)
  IF FH(i) > HX THEN HX = FH(i)
  IF FH(i) < HN THEN HN = FH(i)
  FC(i) = P(i) * C(i) / (MR * TC)
  IF FC(i) > CY THEN CY = FC(i)
  IF FC(i) < CN THEN CN = FC(i)
NEXT i' continue
'If FH(i) and FC(i) are graphed as a function of the angle, it is seen
'that a good approximation of the graph is to have two periods per
'cycle of constant mass flow interspersed with periods of no flow at
'all. F1 to F4 are the fractions of the total cycle time when
'different flows are assumed to occur (see nomenclature).
'When 'FH1' and 'FC1' are calculated, the average cycle time, when flow
'is assumed to occur either into or out of the hot space and either
'into or out of the cold space, is calculated.
F1 = (HX - HN) / (6 * (FH(1) - FH(3)))
F2 = (HX - HN) / (6 * (FH(10) - FH(8)))
F3 = (CY - CN) / (6 * (FC(8) - FC(10)))
F4 = (CY - CN) / (6 * (FC(3) - FC(1)))
FH1 = (F1 + F2) / 2
FC1 = (F3 + F4) / 2
'C --- PAG. 340
'Effective mass flow into or out of the hot space is calculated.
M = MR / r
WH = (HX - HN) * M * MW * NU / FH1
'Effective mass flow into or out of the cold space is calculated.
WC = (CY - CN) * M * MW * NU / FC1
'Fraction of the time the flow is assumed to pass through the
'regenerator and the flow rate of the regenerator is calculated as the
'average between the hot and cold flows.
FR = (FH1 + FC1) / 2
WR = (WH + WC) / 2
'Regenerator gas density.
RM = .1202 * MW * PG / TR
'Calculates regenerator windage loss.
MU = M1 + M2 * (TR - 293) + M3 * PG
GR = WR / (P4 * DR ^ 2 * NR%)
RR = DW * GR / MU
CW = 2.7312 * (1 + 10.397 / RR)
DP = CW * GR ^ 2 * LR / (2E+07 * DW * RM)
A = 1 / RM
RW = DP * WR * 2 * FR * A
'Calculates heater windage loss. In this calculation the viscosity for
'the input temperature and subroutine 'rest' returns the friction
'factor for the input reynolds number. The calculation takes into
'account friction losses, as well as 4.4 velocity heads for an
'entrance and an exit loss, one 180 degree bend, and two 90 degree
'bends.
MU = M1 + M2 * (TM - 293) + M3 * PG
RM = .1202 * MW * PG / TM
A = 1 / RM
GH = WH / (P4 * IH ^ 2 * NH%)
RE = IH * GH / MU
RT = RE

IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF

AF = P4 * IH ^ 2 * NH%
'C --- PAG. 341
VH = WH / (RM * AF)
DP = 2 * CW * GH ^ 2 * LH / (1E+07 * IH * RM) + VH ^ 2 * 4.4 * RM / 2E+07
HW = DP * WH * 2 * FH1 * A
'This calculates the windage loss through the gas cooler and the
'connecting tube. The same comments for the gas heating windage loss
'apply here as well. The velocity heads charge to the gas cooler is
'1.5 for a smple entrance and exit loss. In the connecting head line,
'three velocity heads are changed to account for entrance and exit loss
'plus two 90 degree bends.
'NU = M1 + M2 * (TX - 293) + M3 * PG ??? NU is input
MU = M1 + M2 * (TX - 293) + M3 * PG
'MU   gas viscosity, g/cm sec
RM = .1202 * MW * PG / TX
A = 1 / RM
GC = WC / (P4 * IC ^ 2 * NC%)

RE = IC * GC / MU
RZ = RE
IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF
AF = P4 * IC ^ 2 * NC%
VC = WC / (RM * AF)
DP = 2 * CW * GC ^ 2 * LC / (1E+07 * IC * RM) + VC ^ 2 * 1.5 * RM / 2E+07
GD = WC / (P4 * ID ^ 2 * NE%)
RE = ID * GD / MU

IF RE >= 2000 THEN
  X = LOG(RE)
  X = -3.09 - .2 * X
  CW = EXP(X)
ELSE
  CW = 16 / RE
END IF

AF = P4 * ID ^ 2 * NE%
VC = WC / (RM * AF)

IF (MS% AND 8) = 0 THEN
  DP = DP + 2 * CW * GD ^ 2 * LE / (1E+07 * ID * RM) + VC ^ 2 * 3 * RM / 2E+07
ELSE
'Italian
  DP = DP + 2 * CW * GC ^ 2 * LE / (1E+07 * ID * RM) + VC ^ 2 * 3 * RM / 2E+07
END IF

IF DP * WC * 2 * FC1 * A > 1E+20 THEN
  mes$ = "< EROR > Insufficient cool pipes number or diameter."
  RETURN 'GOTO 300
END IF
CF = DP * WC * 2 * FC1 * A
'Calculates indicated power.
IP = BP - HW - RW - CF
'C --- PAG. 342
'Calculates mechanical friction loss.
MF = (1 - ME / 100) * IP
'Calculates net power.
NP = IP - MF
'Calculates basic heat input.
BH = BP / (1 - TC / TH)
'Calculates reheat loss for net net .05-.20 which is used in the 4l23
'machine. This section is specific for this type of regenerator
'material.

IF RR < 42 THEN
  IF RR <= 0 THEN
    mes$ = "< ERROR > Regenerator Reynolds number <0 (RR)."
    RETURN
  END IF
'  IF (MS% AND 64) = 0 THEN
    X = EXP(-.1826 - .0583 * LOG(RR))
'  ELSE
'   'Italian
'    X = EXP(-.1826 - .5044 * LOG(RR))
   'I think this is a mistake because the factor .5044
   'is also in an other line.
'  END IF
ELSEIF RR < 140 THEN
  X = EXP(.5078 - .243 * LOG(RR))
ELSE
  X = EXP(1.78 - .5044 * LOG(RR))
END IF

NT = X * LR / DW
X = WR * CP * (TH - TW)
Y = RD * CV * (PX - PN) * NU * MW / (r * FR)
K3 = FR * (X - Y)
RH = K3 / (NT + 2) * 2
'Calculates temperature swing loss
MX = NR% * P4 * DR ^ 2 * LR * FF * 7.5
TS = K3 / (NU * MX * 1.05)
SL = K3 * TS / (2 * (TH - TX))
'Calculates pumping or appendix loss.
X = (PI * DC / KG) ^ .6
'Italian C following line does not match origina code
'But is the same as in the NASA document
Y = ((PX - PN) * MW * NU * CP * 2 / ((TH + TX) * r)) ^ 1.6
Z = G ^ 2.6
QP = X * 2 * LB * (TH - TX) * Y * Z / 1.5
'Calculates shuttle heat loss.
QS = 2 * P4 * RC * RC * KG * (TH - TC) * DC / (G * LB)
'Calculates static heat loss. This can be either specified or
'calculated from the basic dimensions.
IF ZZ% = 1 THEN
  ZHMETAL = (TH - TC) * (KM * ((DR ^ 2 * P4 * FF + PI * DR * SR) / LR + PI * DC * (SC + SE) / LB))
  ZHGAS = (TH - TC) * (KG * (DR ^ 2 * P4 * (1 - FF) / LR + DC ^ 2 * P4 / LB))
  ZH = 3.2608033# * (ZHMETAL + ZHGAS)
  'ZH=static heat conduction loss (the specified loss was 13 times lower)
ELSEIF (MS% AND 32) = 32 THEN
 'Correction for different temperatures.
  ZH = (TH - TC) * ZH / 453.79452#
END IF

'Sums all losses to calculate net heat demand.
QN = BH + ZH + SL + RH - HW - RW / 2 + QS + QP

'C --- PAG. 343
'Calculates cooler heat load.
QC = QN - NP
'Temperature rise in cooling water.
DT = QC / (FW * 4.185)
'Effective cold metal temperature.
TX = TW + DT / 2
'Calculates heat transfer coefficient in the cold heat exchanger.
RE = RZ
J = 1
'goto subroutine rest
GOSUB 100
IF mes$ <> "" THEN RETURN

HC = ST * CP * GC / PR

 'Two different methods of arriving at the proper effective hot space
 'and cold space temperature are interspersed. The fastest way,
 'which is usually tried first, involves calculating what the
 'temperature difference has to be between the metal temperature and
 'the effective gas temperature considering the heat transfer
 'capability of the heat exchanger and the correction factor.
 'However, if the heat exchanger is too small, the first iteration
 'method goes unstable and a second, more cautious, method must be
 'employed. The 'ZA%' is the flag which shows that the second
 'method is caled in.
  IF ZA% = 1 THEN
    GOSUB 46
    RETURN
  END IF
 ' 'X' is used as a temporary variable for the previous cold
 'temperature. The cold temperature is calculated, assuming there is
 'no error between the heat that can be transferred and the heat that
 'should be transfered. Counter 'ZB%' is indexed. A test is now made
 'of the 'TC' value just calculated. If the effective cold gas
 'temperature is greater than the effective hot gas temperature or
 'less than the cooling water temperature this iteration method has
 'gone unstable and the second, more cautious, method is brought in.
 'Also if the first iteration method has not come to an answer within
 '10 iterations, ('ZB%' greater than 10), the second iteration method
 'is brought in. The initial change in the hot gas temperature, 'DV',
 'and in the cold gas temperature. 'DU', are both set at 64 degrees.
 'The flag 'ZA%' is set at 1 and 'TC' and 'TH' are set at the initial
'C --- PAG. 344
 'values. Control pass to 46 where the second approach begins. If
 'the value of 'TC' does not indicate the second approach is needed
 'control passes to 48 to start calculation of the effective
 'temperature in the hot space.
X = TC
YY = HC * FC1 * AC * BF

TC = QC / YY + TX
E2 = QC - YY * (TC - TX)
ZB% = ZB% + 1
IF (TC > TH) OR (TC < TX) OR (ZB% > 10) THEN
  GOSUB 47
  RETURN
END IF
GOSUB 48
RETURN


47 DV = 64
DU = 64
ZA% = 1
TC = TW
TH = TM
GOSUB 46
RETURN


'On the first time through 'TC' = 'TW' and the error in the cold space,
'E2, is made equal to the required heat transfer through the gas
'coolers, 'QC'. Then the next estimate for 'TC' is made by adding
' 'DU', 64 degrees, to 'TX', the average temperature of the gas
'cooler metal. The program then goes to 48, skipping over the rest of
'the adfusiment program for the cold space.
46 IF TC = TW THEN
  'Location of control if 'TC' equals 'TW'.
49 E2 = QC
  TC = TX + DU
  GOSUB 48
  RETURN
END IF

'if 'TC' is not equal to 'TW', as it will be for anything except
'for the first time through, the previous error is saved as 'E1'.
'than 'E2' is calculated as the difference between the heat that
'should be transferred and the heat that can be transferred by the
'capabilities of the heat exchanger.
E1 = E2

E2 = QC - HC * FC1 * AC * (TC - TX) * BF

'if this error is positive, then the correction number, 'DU', is
'added to the cold temperature, 'TC', and the program goes on to the
'hot space analysis.
IF E2 > 0 THEN
' Location of control if 'E2' is greater than 0.
  TC = TC + DU
  GOSUB 48
  RETURN
END IF

'If this error is negative and the previous error was positive,
'then de degree increment 'DU', is just divided by 4, for future
'corrections.
IF (E2 < 0) AND (E1 > 0) THEN DU = DU / 4
'The degree increment is subtracted from 'TC'. If 'TC' becomes
'greater than 'TM', the hot metal temperature, obviously there is
'insufficient cooler heat transfer area and the programm stops for
'this case. This can occur for small cooler areas and specified heat
'C --- PAG. 345
'leaks.
TC = TC - DU

IF skipit <> true THEN
PRINT "TC="; TC
PRINT "DU="; DU
PRINT "TM="; TM
PRINT "tc-tm="; TC - TM

INPUT i$
END IF

IF TC > TM THEN
 'Because of insufficient cooler area the program is terminated for
 'this case.

  mes$ = "< ERROR > Insufficient cooler area."
  
  RETURN 'GOTO 300
END IF
GOSUB 48
RETURN

'Calculates heat transfer coefficient for gas heater. Flag 'ZA%'
'indicates whether the fast method of convergence at 59 or the slow
'method at 52 should be used.
48 RE = RT
J = 2
'goto subroutine rest
GOSUB 100
IF mes$ <> "" THEN RETURN

HH = ST * CP * GH / PR
'HH   heat transfer coefficient in heater, w/cm^2 K
'ST   Stanton number X(Pr)^(2/3)
'CP  heat capacity of hydrogen at constant p=14.62 j/g K @ 700 K (assumed
'    not to vary importantly with temperature) (What about O2 and He?)
'GH   mass velocity in heater, g/sec cm^2
'PR   Prandtl number to the 2/3 power = (Pr)^(2/3)

IF ZA% = 1 THEN

'This is analogous to 46 to 48, except this is for the hot space.
  IF TH = TM THEN
   'Location of control if 'TH' equals 'TH'.
    E4 = QN
    TH = TM - DV
  ELSE
    E3 = E4
    E4 = QN - HH * FH1 * AH * (TM - TH) * BF

    IF E4 <= 0 THEN
      IF (E4 < 0) AND (E3 > 0) THEN DV = DV / 4
        TH = TH + DV
        IF TH < TW THEN
         'Because of insufficient heater area the program is terminated for
         'this case
          mes$ = "< ERROR > Insufficient heater area."
          PRINT mes$
          INPUT i$
          RETURN 'GOTO 300
        END IF
      END IF
'C --- PAG. 347
    ' Location of control if 'TH' is not less than 'TM'.
      TH = TH - DV
    END IF

'   Convergence criteria for the slower method of iteration.
'   Convergence is complete when the air in the hot space and the air in
'   the cold space are both less than 1% of the heat transfered through
'C --- PAG. 346
'   the heat exchangers.
    X1 = ABS(E4)
    X2 = ABS(E2)
    x3 = QN / 100
    x4 = QC / 100
    IF (X1 > x3) OR (X2 > x4) THEN
      GOSUB 200 '?????
      RETURN
    END IF
 '  Location of control for the second iteration method.
    GOSUB 47
    RETURN
  END IF

'This is analogous to the coment made after 44 on the cold space,
'except this is for the hot space.
  Y = TH
  YY = HH * FH1 * AH * BF

  TH = TM - QN / YY
  E4 = QN - YY * (TM - TH)
  IF (TH > TM) OR (TH < TC) THEN
    GOSUB 47
    RETURN
  END IF
 'Convergence criteria for the first iteration method. The iteration
'is complete when change in the effective hot space and cold space
'temperature is less than one degree kelvin per iteration.
  X1 = ABS(TH - Y)
  X2 = ABS(TC - X)
  IF (X1 > 1) OR (X2 > 1) THEN
    GOSUB 200
    RETURN
  END IF

'completes preparation for output
A = -HW - RW / 2
B = 100 * IP / QN
C1 = QN * (100 / FE - 1)
D = FE * NP / QN
E = 100 * QN / FE

BP = BP * N%
HW = HW * N%
RW = RW * N%
CF = CF * N%
IP = IP * N%
MF = MF * N%
NP = NP * N%
BH = BH * N%
RH = RH * N%
QS = QS * N%
QP = QP * N%
SL = SL * N%
ZH = ZH * N%
A = A * N%
QN = QN * N%
C1 = C1 * N%
E = E * N%

IF skipit <> true THEN GOSUB 1000 ' print outputs

IF MID$(ss$, 10, 1) = "*" THEN GOSUB 1170'save outputs in file

'Reinitializing
i = i + 1
ZA% = 0
ZB% = 0

RETURN 'GOTO 300

100 'Subroutine REST
'Calculates stanton number from reynolds number
IF RE >= 10000 THEN
  ST = EXP(-3.57024 - .229496 * LOG(RE))
ELSEIF RE >= 7000 THEN
  ST = .0034
ELSEIF RE >= 4000 THEN
  ST = EXP(-13.3071 + .861016 * LOG(RE))
ELSEIF RE >= 3000 THEN
'C --- PAG. 348
  ST = .0021
ELSE
  'gives an error when optimizing regenerators
  PRINT "RE="; RE
  IF RE <= 0 THEN
    IF J = 1 THEN
      mes$ = "< ERROR > Insufficient regenerators, number or diameter.0>=(RE)"
    ELSE
      mes$ = "< ERROR > Insufficient regenerators, number or diameter 0>=(RE)"
    END IF
    RETURN
  END IF
  ST = EXP(.337046 - .812212 * LOG(RE))
END IF
RETURN


1000 ' print outputs

'This is where the printing of the output starts. To compress output
'the operating conditions and engine dimensions are identified only by
'their fortran symbol.

'Prints program heading
'Prints current operating conditions

PRINT "Current operating conditions are:"
PRINT "SP="; SP, "PS="; PS, "ND="; ND, "TF="; TF
PRINT "L1="; L1, "TY="; TY, "FX="; FX, "OG%="; OG%
PRINT
'Prints current dimensions
PRINT "current dimensions are:"
PRINT "DC="; DC, "DR="; DR, "IC="; IC, "OC="; OC
PRINT "DW="; DW, "DD="; DD, "IH="; IH, "OH="; OH
PRINT "G="; G, "LB="; LB, "LR="; LR, "CR="; CR
PRINT "RC="; RC, "LC="; LC, "LD="; LD, "LH="; LH

PRINT "LI="; LI, "NC%="; NC%, "NR%="; NR%, "N%="; N%
PRINT "NH%="; NH%, "FF="; FF, "AL="; AL, "CX="; CX
PRINT "ME="; ME, "FE="; FE, "EC="; EC, "SC="; SC
PRINT "SE="; SE, "SR="; SR, "ZZ%="; ZZ%, "ZH="; ZH
PRINT "KM="; KM, "ID="; ID, "LE="; LE, "NE%="; NE%
PRINT "BF="; BF
PRINT
GOSUB 1188 'Hit key

'Prints power outputs and heat inputs
PRINT "power, watts", , "heat requirement, watts"
PRINT "basic", BP, "basic", BH
PRINT "heater f.l.", HW, "reheat", RH
PRINT "regen. f.l", RW, "shuttle", QS
PRINT "cooler f.l.", CF, "pumping", QP
PRINT "net       ", IP, "temp,swing", SL
PRINT "mech.fric.", MF, "conduction", ZH
PRINT "brake", NP, "flow fric.credit"; A
PRINT "---------------------------",
PRINT "heat to engine"; QN
PRINT "indicated eff.%="; B, "furnace loss", C1
PRINT "overall   eff.%="; D, "fuel input", E
PRINT "----------------------------------------------------------------------"
PRINT "hot metal temp.   k="; TM; "  cooling water inlet; temp., k ="; TW
PRINT "effec.hot sp.temp.k="; TH; "  effec.cold sp.temp.k.="; TC
PRINT "----------------------------------------------------------------------"

GOSUB 1188 'Hit key


'Prints work diagram from data

PRINT "final work diagram:"
'PRINT "angle", "hot vol", "cold vol.", "tot. vol,", "pressure"; "  gas inv."
PRINT "ANGLE";
PRINT RIGHT$("             " + "HOT VOL.", 13); " ";
PRINT RIGHT$("             " + "COLD VOL.", 12);
PRINT RIGHT$("             " + "TOT. VOL.", 13);
PRINT RIGHT$("             " + "PRESSURE", 13);
PRINT RIGHT$("             " + "GAS INV.", 13)

FOR i = 1 TO 13
  F = ND * i - 30
  PRINT RIGHT$("     " + STR$(F), 5);
  CALL makestr(H(i), printv$)
  PRINT printv$;
  CALL makestr(C(i), printv$)
  PRINT printv$;
  CALL makestr(V(i), printv$)
  PRINT printv$;
  CALL makestr(P(i), printv$)
  PRINT printv$;
  CALL makestr(L(i) / r, printv$)
  PRINT printv$
NEXT i
GOSUB 1188 'Hit key

CLS
RETURN

1100 'start screen
IF ss$ = "" THEN ss$ = "  **      "
IF setfilename$ = "" THEN setfilename$ = "STIRLSET.TXT"
IF inifilename$ = "" THEN inifilename$ = "STIRLSET.TXT"
IF inpfilename$ = "" THEN inpfilename$ = "STIRLINP.TXT"
IF newfilename$ = "" THEN newfilename$ = "STIRLINP.TXT"
IF outfilename$ = "" THEN outfilename$ = "STIRLOUT.TXT"

CLS
PRINT "< START SCREEN >                                         'tab' or '->' for next";
PRINT
PRINT "                   =================================="
PRINT "                   = STIRLING ENGINE DESIGN PROGRAM ="
PRINT "                   =================================="
PRINT "              On the bases of the Wiliam R Martini program:"
PRINT "                   Isothermal second order calculation"
PRINT "        Converted to Quick Basic, changed and enhanced with an optimizer"
PRINT "             by Harrie van der Haghen. hvd@haghen.nl Amsterdam 2013"
IF mes$ = "" THEN
PRINT "  ############################################################################"
ELSE
  PRINT " MESSAGE: ===<<  "; mes$; "  >>==="
  mes$ = ""
END IF
PRINT
PRINT "DO:  START UP SETTINGS:                          FILENAMES:"
PRINT "     > START CALCULATIONS >>>"
PRINT "     Skip extra information on screen.", ,
PRINT "  switch is ";
IF skipit = true THEN
  PRINT "ON"
ELSE
  PRINT "OFF"
END IF
PRINT "     Start opstimizer in single step with space bar mode.";
PRINT " switch is ";
IF spacebarstep = true THEN
  PRINT "ON"
ELSE
  PRINT "OFF"
END IF
PRINT "     Optimize on all records in input file.",
PRINT "  switch is ";
IF optimizeall = true THEN
  PRINT "ON"
ELSE
  PRINT "OFF"
END IF

PRINT "     Load original Martini mode with exemple GM4L23 inputs."
PRINT "     Optimized mode and load changed GM4L23 inputs."
'PRINT
PRINT "     Load settings from file: settings file name= "; inifilename$
PRINT "     Save settings to file:   settings file name= "; setfilename$
PRINT "     Load inputs from file:      input file name= "; inpfilename$
PRINT "     Save inputs to file:    new input file name= "; newfilename$
PRINT "     Save outputs in file:      output file name= "; outfilename$

PRINT "---------------------------------------------------------------------"
PRINT "  Select with arrow keys, change with * or =, Enter= do, ?= help, q= quit.";
IF sselect = 0 THEN sselect = 1
FOR pline = 1 TO 11
' IF pline < 6 THEN
   LOCATE pline + 11, 1
' ELSE
 '  LOCATE pline + 12, 1
' END IF
 IF sselect = pline THEN
   PRINT "-> ";
 ELSE
   PRINT "   ";
 END IF
 IF MID$(ss$, pline, 1) = "*" THEN
   PRINT "*";
 ELSE
   PRINT " ";
 END IF
NEXT pline

1120 i$ = INKEY$
IF i$ = "" THEN GOTO 1120

IF i$ = "q" OR i$ = "Q" THEN
  END' quit program
ELSEIF i$ = "+" OR (i$ = CHR$(0) + CHR$(80)) THEN
  sselect = sselect + 1
  IF sselect > 11 THEN sselect = 0
ELSEIF i$ = "-" OR (i$ = CHR$(0) + CHR$(72)) THEN
  sselect = sselect - 1
  IF sselect < 1 THEN sselect = 11
ELSEIF i$ = "*" OR i$ = " " THEN
  IF MID$(ss$, sselect, 1) = " " THEN
    MID$(ss$, sselect, 1) = "*"
    IF sselect = 5 THEN
      MID$(ss$, 6, 1) = " "
    ELSEIF sselect = 6 THEN
      MID$(ss$, 5, 1) = " "
    END IF
  ELSE
    MID$(ss$, sselect, 1) = " "
  END IF
ELSEIF i$ = "=" THEN
1121 'new input
  LOCATE 3, 1
  FOR pline = 5 TO 10
   PRINT "                                                                              ";
  NEXT pline
  LOCATE 3, 1
  PRINT mes$
  PRINT ">>>>> INPUT NEW FILE NAME >> ";
  INPUT filename$
  CALL dosfilename(filename$, mes$)  'cheque filename
  IF mes$ <> "" THEN GOTO 1121'new input

  IF sselect = 7 THEN
    inifilename$ = filename$
    CALL givefilenr(filecount%, filenr%)
    filename$ = "STIRLSYS.TXT"
    OPEN filename$ FOR OUTPUT AS filenr%
    PRINT #filenr%, CHR$(34); inifilename$; CHR$(34)
    CALL closefilenr(filecount%, filenr%)
    CLOSE filenr%
  ELSEIF sselect = 8 THEN
    setfilename$ = filename$
  ELSEIF sselect = 9 THEN
    inpfilename$ = filename$
  ELSEIF sselect = 10 THEN
    newfilename$ = filename$
  ELSEIF sselect = 11 THEN
    outfilename$ = filename$
  END IF
ELSEIF i$ = "?" THEN
  CALL givefilenr(filecount%, varfilenr%)
  varfile$ = "STIRLD.HLP"
  PRINT "varfilenr%="; varfilenr%, "varfile$="; varfile$
  OPEN varfile$ FOR INPUT AS varfilenr%
  CLS
  PRINT
  line$ = ""
  WHILE line$ <> "start"
   INPUT #varfilenr%, line$
  WEND
  IF sselect < 11 THEN
    WHILE line$ <> ("START" + CHR$(48 + sselect))
     INPUT #varfilenr%, line$
    WEND
  ELSE
    WHILE line$ <> ("START1" + CHR$(38 + sselect))
     INPUT #varfilenr%, line$
    WEND
  END IF
  PRINT "INFORMATION ABOUT START SCREEN OPTION:"
  INPUT #varfilenr%, line$
  PRINT line$
  PRINT "----------------------------------------------------------------"
  INPUT #varfilenr%, line$
  PRINT line$
  CALL closefilenr(filecount%, varfilenr%)
  PRINT

  GOSUB 1188 'Hit key
  
ELSEIF i$ = CHR$(0) + CHR$(75) THEN
  spacebarstep = true
  GOTO 1241 'optimizer screen
ELSEIF (i$ = CHR$(0) + CHR$(77)) OR (ASC(i$) = 9) THEN
  GOTO 2008'inputs screen
ELSEIF i$ = CHR$(13) THEN ' Enter
  donow = sselect
  GOSUB 1125 'execute start screen selected option
  IF do$ = "calculate" THEN RETURN
END IF
GOTO 1100

1124 'execute start settings
FOR donow = 2 TO 11
' IF donow = 6 THEN donow = 9
 IF MID$(ss$, donow, 1) = "*" THEN
   GOSUB 1125
 END IF
NEXT donow
IF MID$(ss$, 1, 1) = "*" THEN
  do$ = "calculate"' start calculations
END IF
RETURN

1125 'execute start screen selected option
  IF donow = 1 THEN
    do$ = "calculate"' start calculations
  ELSEIF donow = 2 THEN
    IF sselect = 2 THEN
      IF skipit = true THEN
        skipit = false
      ELSE
        skipit = true
      END IF
    ELSE
      skipit = true
    END IF
  ELSEIF donow = 3 THEN
    IF sselect = 3 THEN
      IF spacebarstep = true THEN
        spacebarstep = false
      ELSE
        spacebarstep = true
      END IF
    ELSE
      spacebarstep = true
    END IF
  ELSEIF donow = 4 THEN
    IF sselect = 4 THEN
      IF optimizeall = true THEN
        optimizeall = false
      ELSE
        optimizeall = true
      END IF
    ELSE
      optimizeall = true
    END IF
  ELSEIF donow = 5 THEN
    GOSUB 2050 'read exemple values
    GOSUB 1188 'Hit key
    MS% = 0
    GOSUB 2022 'load input values in vi()
  ELSEIF donow = 6 THEN
    GOSUB 2050 'read exemple values
    LB = 10.03
    SC = .4305
    SE = 1.25
    SR = .41
    LR = 2.79
    MS% = 255
    GOSUB 2022 'load input values in vi()
    mes$ = "Changed GM4L23 inputs loaded."
  ELSEIF donow = 7 THEN
    filename$ = inifilename$
    CALL existfile(filename$, filecount%, mes$)
    IF mes$ = "" THEN
      GOSUB 1340 'read settings file
    END IF
  ELSEIF donow = 8 THEN
    GOSUB 1330'save settings file
  ELSEIF donow = 9 THEN
    CALL existfile(filename$, filecount%, mes$)
    IF mes$ = "" THEN
      GOSUB 1160 'load inputs from file
    END IF
  ELSEIF donow = 10 THEN
    GOSUB 1150 'save inputs in file
  ELSEIF donow = 11 THEN
    GOSUB 1170 'save outputs in file
  END IF
RETURN

1150 'save inputs in file
CALL givefilenr(filecount%, newfilenr%)
LOCATE 10, 1
PRINT "saving inputs in file newfilenr%="; newfilenr%;
PRINT " newfilename$="; newfilename$; " ";
OPEN newfilename$ FOR OUTPUT AS newfilenr%
PRINT #newfilenr%, CHR$(34); ".null."; CHR$(34)
PRINT #newfilenr%, DC; ","; LC; ","; LD; ","; IC; ","; OC; ","; NC%
PRINT #newfilenr%, DW; ","; FX; ","; ME; ","; FE; ","; OG%; ","; ZZ%
PRINT #newfilenr%, ZH; ","; LH; ","; LI; ","; IH; ","; OH; ","; NH%; ","; DD
PRINT #newfilenr%, RA; ","; G; ","; LB; ","; PS; ","; KM; ","; SC; ","; SE
PRINT #newfilenr%, SR; ","; LR; ","; DR; ","; NR%; ","; FF; ","; CR; ","; RC
PRINT #newfilenr%, N%; ","; AL; ","; TF; ","; TY; ","; SP; ","; AA; ","; BA
PRINT #newfilenr%, ID; ","; LE; ","; NE%; ","; BF; ","; PP; ","; CM; ","; FQ
PRINT #newfilenr%, r; ","; HP; ","; EC; ","; L1; ","; ASAS
PRINT #newfilenr%, MS%
PRINT #newfilenr%, "bottom"
CALL closefilenr(filecount%, newfilenr%)
mes$ = "Save inputs in file " + newfilename$ + " : Ready"

RETURN

1160 'open inputs file
CLS

IF inpfilenr% = 0 THEN
  CALL givefilenr(filecount%, inpfilenr%)
ELSE
  CLOSE inpfilenr%
END IF

LOCATE 10, 1
PRINT "loading inputs from file filenr%="; inpfilenr%;
PRINT " inpfilename$="; inpfilename$; " "

OPEN inpfilename$ FOR INPUT AS inpfilenr%
inpcols = 0
inprows = 0
start$ = ""
WHILE (NOT EOF(inpfilenr%)) AND ((start$ <> ".null.") AND (start$ <> "start"))
 INPUT #inpfilenr%, start$
WEND

IF start$ = ".null." THEN
  inprows = 1
ELSEIF start$ = "start" THEN
  inpcolvars$ = ""
  op$ = "                       "
  INPUT #inpfilenr%, inpcols
  inpcolvars$ = ""
  FOR inpcol = 1 TO inpcols
   INPUT #inpfilenr%, varname$
   PRINT "inpcol"; inpcol, varname$,
   FOR var = 1 TO nivar
    IF varname$ = ni$(var) THEN
      inpcolvars$ = inpcolvars$ + CHR$(var)
      PRINT "input var="; varname$
      EXIT FOR
    END IF
   NEXT var
   IF var > nivar THEN
     FOR var = 1 TO novar
      IF varname$ = no$(var) THEN
        inpcolvars$ = inpcolvars$ + CHR$(var + nivar)
        MID$(op$, var, 1) = "*"
        PRINT "output var="; varname$
        EXIT FOR
      END IF
     NEXT var
   END IF
  NEXT inpcol
  GOSUB 1188 'Hit key
  CLS
  inprows = 0

  inp$ = ""
  WHILE NOT EOF(inpfilenr%) AND inp$ <> "end"
   inp$ = ""
   WHILE NOT EOF(inpfilenr%) AND inp$ = ""
    INPUT #inpfilenr%, inp$
   WEND
   IF (NOT EOF(inpfilenr%)) AND (inp$ <> "") AND (inpcols > 1) AND inp$ <> "end" THEN
     inprows = inprows + 1
     PRINT : PRINT "inprows="; inprows, inp$;
     FOR inpcol = 2 TO inpcols
      INPUT #inpfilenr%, inp$
      PRINT " "; inp$;
     NEXT inpcol
   ELSE
     IF inpcols = 1 AND inp$ <> "end" THEN inprows = inprows + 1
   END IF
   IF inprows MOD 20 = 0 THEN
     GOSUB 1188 'Hit key
     CLS
   END IF
  WEND
END IF
PRINT : PRINT "ready"
'GOSUB 1188 'Hit key

IF inprows = 0 THEN

  CALL closefilenr(filecount%, inpfilenr%)
'  CLOSE inpfilenr%
  mes$ = "-- ERROR -- No '.null.' or 'start' in input file "
  mes$ = mes$ + inpfilename$
ELSE
  inprow = 1
  GOSUB 1162 'read next row entries from inputfile
END IF
RETURN

1162 'read next row entries from inputfile
IF inprow = 1 THEN
  CLOSE inpfilenr%
  OPEN inpfilename$ FOR INPUT AS inpfilenr%
END IF

IF start$ = ".null." THEN
  PRINT "START$="; start$
  inp$ = ""
  WHILE (NOT EOF(inpfilenr%)) AND (inp$ <> ".null.")
   INPUT #inpfilenr%, inp$
  WEND

  INPUT #inpfilenr%, DC, LC, LD, IC, OC, NC%
  INPUT #inpfilenr%, DW, FX, ME, FE, OG%, ZZ%
  INPUT #inpfilenr%, ZH, LH, LI, IH, OH, NH%, DD
  INPUT #inpfilenr%, RA, G, LB, PS, KM, SC, SE
  INPUT #inpfilenr%, SR, LR, DR, NR%, FF, CR, RC
  INPUT #inpfilenr%, N%, AL, TF, TY, SP, AA, BA
  INPUT #inpfilenr%, ID, LE, NE%, BF, PP, CM, FQ
  INPUT #inpfilenr%, r, HP, EC, L1, ASAS
'  INPUT #inpfilenr%, MS%
'  INPUT #inpfilenr%, last$
  GOSUB 2022 'load input values in vi()

  mes$ = "Load inputs from file " + inpfilename$ + " : Ready"
  CALL closefilenr(filecount%, inpfilenr%)

ELSE
  IF inprow = 1 THEN
    CLOSE inpfilenr%
    OPEN inpfilename$ FOR INPUT AS inpfilenr%
    inp$ = ""
    WHILE inp$ <> "start"
     INPUT #inpfilenr%, inp$
    WEND
    FOR inpcol = 0 TO inpcols
     INPUT #inpfilenr%, inp$
    NEXT inpcol
 END IF
 inp$ = ""
 WHILE inp$ = ""
  INPUT #inpfilenr%, inp$
 WEND

 IF ASC(inpcolvars$) <= nivar THEN
  vi(ASC(inpcolvars$)) = VAL(inp$)
 ELSE
  vo(ASC(inpcolvars$) - nivar) = VAL(inp$)
 END IF

 IF inpcols > 1 THEN
  FOR var = 2 TO inpcols
   IF ASC(MID$(inpcolvars$, var, 1)) <= nivar THEN
     INPUT #inpfilenr%, vi(ASC(MID$(inpcolvars$, var, 1)))
   ELSE
     INPUT #inpfilenr%, eo$(ASC(MID$(inpcolvars$, var, 1)) - nivar)
   END IF
  NEXT var
 END IF
END IF
RETURN

1170 'save outputs in file
CALL givefilenr(filecount%, outfilenr%)
LOCATE 10, 1
PRINT "saving inputs in file outfilenr%="; outfilenr%;
PRINT " outfilename$="; outfilename$; " ";
OPEN outfilename$ FOR OUTPUT AS outfilenr%

PRINT #outfilenr%, "current operating conditions are:"
PRINT #outfilenr%, "SP = "; SP, "PS = "; PS, "ND = "; ND, " TF = "; TF
PRINT #outfilenr%, "L1 = "; L1, "TY = "; TY, "FX = "; FX, " OG = "; OG%
PRINT #outfilenr%,
PRINT #outfilenr%, "current dimensions are:"
PRINT #outfilenr%, "DC="; DC, "DR="; DR, "IC="; IC, "OC="; OC
PRINT #outfilenr%, "DW="; DW, "DD="; DD, "IH="; IH, "OH="; OH
PRINT #outfilenr%, "G="; G, "LB="; LB, "LR="; LR, "CR="; CR
PRINT #outfilenr%, "RC="; RC, "LC="; LC, "LD="; LD, "LH="; LH

PRINT #outfilenr%, "LI="; LI, "NC%="; NC%, "NR%="; NR%, "N%="; N%
PRINT #outfilenr%, "NH%="; NH%, "FF="; FF, "AL="; AL, "CX="; CX
PRINT #outfilenr%, "ME="; ME, "FE="; FE, "EC="; EC, "SC="; SC
PRINT #outfilenr%, "SE="; SE, "SR="; SR, "ZZ%="; ZZ%, "ZH="; ZH
PRINT #outfilenr%, "KM="; KM, "ID="; ID, "LE="; LE, "NE%="; NE%
PRINT #outfilenr%, "BF="; BF
PRINT #outfilenr%, "------------------------------------------------"
PRINT #outfilenr%, "PI="; PI, "P4="; P4, "FX="; FX
PRINT #outfilenr%,
'Prints power outputs and heat inputs
PRINT #outfilenr%, "power, watts", , "heat requirement, watts"
CALL makestr(BP, printv$)
PRINT #outfilenr%, "basic", printv$,
CALL makestr(BH, printv$)
PRINT #outfilenr%, "basic", , printv$
CALL makestr(HW, printv$)
PRINT #outfilenr%, "heater f.l.", printv$,
CALL makestr(RH, printv$)
PRINT #outfilenr%, "reheat", , printv$
CALL makestr(RW, printv$)
PRINT #outfilenr%, "regen. f.l", printv$,
CALL makestr(QS, printv$)
PRINT #outfilenr%, "shuttle", , printv$
CALL makestr(CF, printv$)
PRINT #outfilenr%, "cooler f.l.", printv$,
CALL makestr(QP, printv$)
PRINT #outfilenr%, "pumping", , printv$
CALL makestr(IP, printv$)
PRINT #outfilenr%, "net       ", printv$,
CALL makestr(SL, printv$)
PRINT #outfilenr%, "temp,swing", , printv$
CALL makestr(MF, printv$)
PRINT #outfilenr%, "mech.fric.", printv$,
CALL makestr(ZH, printv$)
PRINT #outfilenr%, "conduction", , printv$
CALL makestr(NP, printv$)
PRINT #outfilenr%, "brake", printv$,
CALL makestr(A, printv$)
PRINT #outfilenr%, "flow fric.credit", printv$
PRINT #outfilenr%, "---------------------------",
CALL makestr(QN, printv$)
PRINT #outfilenr%, "heat to engine", printv$
CALL makestr(C1, printv$)
PRINT #outfilenr%, "indicated eff.%="; B, "furnace loss", , printv$
CALL makestr(E, printv$)
PRINT #outfilenr%, "overall   eff.%="; D, "fuel input", , printv$
PRINT #outfilenr%, "----------------------------------------------------------------------"
PRINT #outfilenr%, "hot metal temp.   k="; TM; "  cooling water inlet; temp., k ="; TW
PRINT #outfilenr%, "effec.hot sp.temp.k="; TH; "  effec.cold sp.temp.k.="; TC
PRINT #outfilenr%, "----------------------------------------------------------------------"

PRINT #outfilenr%,

'Prints work diagram from data

PRINT #outfilenr%, "final work diagram:"
'PRINT "angle", "hot vol", "cold vol.", "tot. vol,", "pressure"; "  gas inv."
PRINT #outfilenr%, "ANGLE";
PRINT #outfilenr%, RIGHT$("             " + "HOT VOL.", 13); "  ";
PRINT #outfilenr%, RIGHT$("             " + "COLD VOL.", 11);
PRINT #outfilenr%, RIGHT$("             " + "TOT. VOL.", 13); " ";
PRINT #outfilenr%, RIGHT$("             " + "PRESSURE", 13); " ";
PRINT #outfilenr%, RIGHT$("             " + "GAS INV.", 13)

FOR i = 1 TO 13
  F = ND * i - 30
  PRINT #outfilenr%, RIGHT$("     " + STR$(F), 5);
  CALL makestr(H(i), printv$)
  PRINT #outfilenr%, printv$;
  CALL makestr(C(i), printv$)
  PRINT #outfilenr%, printv$;
  CALL makestr(V(i), printv$)
  PRINT #outfilenr%, printv$;
  CALL makestr(P(i), printv$)
  PRINT #outfilenr%, printv$;
  CALL makestr(L(i) / r, printv$)
  PRINT #outfilenr%, printv$
NEXT i

CALL closefilenr(filecount%, outfilenr%)
CLOSE outfilenr%
mes$ = "Print outputs to file " + outfilename$ + " : Ready"

RETURN

1188 'Hit key
IF INKEY$ <> "" THEN skipit = false
IF skipit = false THEN
  PRINT
  PRINT "#### Hit any key to continue. ##### (spacebar to skip most outputs)"
1189 i$ = INKEY$
  IF i$ = "" THEN GOTO 1189
  IF i$ = " " THEN skipit = true
END IF
RETURN

'========== optimizer  ==============

1212 'start optimizer
CLS
spacebarstep = true
GOSUB main' calculate outputs

PRINT "IC="; IC
'INPUT i$


do$ = "outpscreen"
IF (INSTR(IP$, "*") > 0) AND (INSTR(op$, "*") > 0) THEN
  FOR var = 1 TO novar
   IF MID$(op$, var, 1) = "*" AND eo$(var) <> "" THEN
     do$ = "optimize"
   END IF
  NEXT var
END IF
IF do$ = "outpscreen" THEN
'PRINT "no optimize": INPUT i$
  RETURN
END IF
'PRINT "yes do optimize": INPUT i$

factor# = inifactor#

loop2 = 0
loop3 = 0
GOSUB 1230 'evaluate

coolnhtl = LC - LD
heatnhtl = LH - LI

evalwas# = eval#

GOSUB 1241 'optimizer screen
IF do$ <> "optimize" THEN RETURN

1215 'loop1
loop1 = loop1 + 1

GOSUB 1225' loop3
loop2 = loop2 + 1

1220 'loop2
loop2 = loop2 + 1

GOSUB 1241 'optimizer screen
IF do$ <> "optimize" THEN RETURN

evalwas# = eval#


1221 ' vari inputs
FOR ivarnr = 1 TO nivar
 IF MID$(IP$, ivarnr, 1) = "*" THEN
   oldeval# = eval#
   remember = vi(ivarnr)

   newvalue = remember * (1 + factor#)
   GOSUB 1223 'change var value

   GOSUB 1225  ' loop3
   IF eval# > oldeval# OR (mes$ <> "") THEN
    
     newvalue = remember / (1 + factor#)
     GOSUB 1223 'change var value

     GOSUB 1225 ' loop3
     IF eval# > oldeval# OR (mes$ <> "") THEN
       eval# = oldeval#

       newvalue = remember
       GOSUB 1223 'change var value
     END IF
   END IF
 END IF
NEXT ivarnr
IF eval# >= evalwas# OR (mes$ = "") THEN

 'IF ABS(factor#) < .00000000000001# THEN
  IF ABS(factor#) < .00000001# THEN
 
    factor# = inifactor#
    IF evalwas# = evalbefor# THEN
      spacebarstep = true
      do$ = ""
mes$ = " ++++++++ ======##################  READY  ##################====== +++++++++"

      GOSUB 1241 'optimizer screen

      RETURN
    ELSE
      evalbefor# = evalwas#
    END IF
  ELSE
    factor# = factor# * tuning#
  END IF
  GOTO 1215'loop1
END IF
GOTO 1220'loop2

1223 'change var value
IF newvalue <> remember THEN
  IF RIGHT$(ni$(ivarnr), 1) = "%" THEN 'integers
    IF newvalue > remember THEN
      newvalue = remember + 1
    ELSE
      newvalue = remember - 1
    END IF
  END IF
END IF

IF newvalue < 0 THEN
  IF (ni$(ivarnr) <> "TF") AND (ni$(ivarnr) <> "TY") THEN
    newvalue = 0
  END IF
END IF

vi(ivarnr) = newvalue
var$ = ni$(ivarnr)
var = ivarnr
GOSUB 2015 'load vi(ivarnr) in input value

IF var$ = "LD" THEN
  vi(ivarnr - 1) = newvalue + coolnhtl
  var$ = ni$(ivarnr - 1)
  var = ivarnr - 1
  GOSUB 2015 'load vi(var) in input value
ELSEIF var$ = "LI" THEN
  vi(ivarnr - 1) = newvalue + heatnhtl
  var$ = ni$(ivarnr - 1)
  var = ivarnr - 1
  GOSUB 2015 'load vi(var) in input value
END IF
RETURN

1225 ' loop3
loop3 = loop3 + 1
eval# = 0

IF optimizeall = true THEN
  FOR inprow = 1 TO inprows
   GOSUB 1162 'read next row entries from input file
   GOSUB main' calculate outputs
   GOSUB 1230'evaluate
  NEXT inprow
ELSE
  GOSUB main' calculate outputs
  GOSUB 1230'evaluate
END IF

RETURN

1230 'evaluate

FOR ovarnr = 1 TO novar
 IF MID$(op$, ovarnr, 1) = "*" AND eo$(ovarnr) <> "" THEN
   IF mes$ <> "" THEN
     IF VAL(eo$(ovarnr)) <> 0 THEN
       eval# = eval# + 100
     ELSE
       eval# = eval# + VAL(eo$(ovarnr))
     END IF
     RETURN
   END IF

   var$ = no$(ovarnr)
   var = ovarnr
   GOSUB 2036 'load value output variable in vo(var)

'  criterium= delta from ideal
'  eval# = eval# + ABS(vo(ovarnr) - VAL(eo$(ovarnr)))

  'criterium= delta % of ideal
  IF vo(ovarnr) <> VAL(eo$(ovarnr)) THEN
    IF VAL(eo$(ovarnr)) <> 0 THEN
      eval# = eval# + ABS(100 - (vo(ovarnr) / (VAL(eo$(ovarnr)) / 100)))
    ELSE
      eval# = eval# + ABS(vo(ovarnr))
    END IF
  END IF

  'criterium= delta factor from ideal
 ' eval# = eval# + ABS(1 - vo(ovarnr) / VAL(eo$(ovarnr)))

  'criterium= factor (deviation from 1) of delta factor from ideal

' PRINT VAL(eo$(ovarnr))
' PRINT eo$(ovarnr)
' PRINT "ovarnr="; ovarnr
' PRINT (ABS(1 - vo(ovarnr)))
'   eval# = eval# + 1 / (ABS(1 - vo(ovarnr) / VAL(eo$(ovarnr))))
 END IF
NEXT ovarnr
RETURN

1241 'optimizer screen
optimum = false
CLS

LOCATE 1, 1
PRINT "< OPTIMIZER SCREEN > space bar = single step, 1=inputs screen, 2=outputs screen,";
LOCATE 2, 1
PRINT " 3=reset factor, 4=toggle skip display info,";
PRINT " Enter=restart, else=continue";
LOCATE 3, 1
PRINT "loop:"; loop1; loop2; loop3; "                 ;"
LOCATE 3, 24
PRINT " factor="; factor#, ' " evalwas#="; evalwas#
LOCATE 3, 54
PRINT " eval="; eval#;
IF mes$ <> "" THEN PRINT mes$
mes$ = ""
LOCATE 5, 1

PRINT "OPTIMIZABLE INPUTS ------------------------------------"
IF INSTR(IP$, "*") = 0 THEN
  PRINT "NO OPTIMIZABLE INPUTS SELECTED IN THE INPUTS SCREEN"
  pline = 6
ELSE
  pline = 5
END IF
FOR var = 1 TO nivar
 IF MID$(IP$, var, 1) = "*" THEN
   pline = pline + 1
  IF pline < 24 THEN
   LOCATE pline, 2
   IF var < 10 THEN PRINT " ";
   PRINT var;
   PRINT ni$(var);
   LOCATE pline, 11
   RESTORE 2002
   FOR abbrev = 1 TO var
     READ abbrev$
   NEXT abbrev
   PRINT abbrev$;
   LOCATE pline, 31
   PRINT vi(var);
  END IF
 END IF
NEXT var
PRINT
PRINT "SELECTED OUTPUTS ----------------------- VALUE ------ OPTIMUM --"
pline = pline + 1
IF INSTR(op$, "*") = 0 THEN
  PRINT "NO OUTPUTS SELECTED IN THE OUTPUTS SCREEN"
  pline = pline + 1
END IF
pline = pline + 1
FOR var = 1 TO novar
 IF MID$(op$, var, 1) = "*" THEN
  pline = pline + 1
  IF pline < 24 THEN
   LOCATE pline, 2
   IF var < 10 THEN PRINT " ";
   PRINT var;
   PRINT no$(var);
   LOCATE pline, 11
   RESTORE 2038
   FOR abbrev = 1 TO var
     READ abbrev$
   NEXT abbrev
   PRINT abbrev$;
   LOCATE pline, 41
  END IF
  var$ = no$(var)
  GOSUB 2036 'load value output variable in vo(var)

  IF pline < 24 THEN

   PRINT vo(var);
   LOCATE pline, 55
  END IF
  IF eo$(var) = "" THEN
    IF pline < 24 THEN
     PRINT "NO OPTIMUM GIVEN";
    END IF
  ELSE
    IF pline < 24 THEN
     PRINT eo$(var)
    END IF
    IF eo$(var) <> "" THEN
       optimum = true
    END IF
  END IF
 END IF
NEXT var
IF optimum = false THEN
  PRINT " NO OPTIMUM VALUE TO OPTIMIZE TO GIVEN IN THE OUTPUTS SCREEN!";
END IF


1242 i$ = INKEY$
IF i$ <> "" THEN
  IF i$ = " " THEN
    spacebarstep = true
    skipit = true
  ELSE
    IF spacebarstep = true THEN
      spacebarstep = false
    END IF
  END IF
ELSE
  IF spacebarstep = false THEN RETURN
END IF

IF i$ = "" THEN GOTO 1242


IF i$ = CHR$(27) THEN
  GOTO 2008'inputs screen
ELSEIF ASC(i$) = 9 OR (i$ = CHR$(0) + CHR$(77)) THEN
  GOTO 1100'start screen
ELSEIF (i$ = CHR$(0) + CHR$(75)) OR i$ = CHR$(15) THEN
  GOTO 2011 'outputs screen

ELSEIF i$ = "1" THEN
  GOTO 2008'inputs screen

ELSEIF i$ = "2" THEN
  GOTO 2011

ELSEIF i$ = "3" THEN
  factor# = inifactor#
  GOTO 1241 'optimizer screen
ELSEIF i$ = "4" THEN
  IF skipit = true THEN
    skipit = false
  ELSE
    skipit = true
  END IF
  GOTO 1241 'optimizer screen
ELSEIF i$ = CHR$(13) THEN
  CLS
  do$ = "calculate"
  RETURN
END IF

RETURN

1330 'save settings file
CALL givefilenr(filecount%, setfilenr%)
LOCATE 10, 1
PRINT "saving settings file setfilenr%="; setfilenr%;
PRINT " setfilename$="; setfilename$; " ";
OPEN setfilename$ FOR OUTPUT AS setfilenr%
PRINT #setfilenr%, CHR$(34); ss$; CHR$(34)
PRINT #setfilenr%, CHR$(34); IP$; CHR$(34)
PRINT #setfilenr%, CHR$(34); op$; CHR$(34)
PRINT #setfilenr%, CHR$(34); inifilename$; CHR$(34)
PRINT #setfilenr%, CHR$(34); setfilename$; CHR$(34)
PRINT #setfilenr%, CHR$(34); inpfilename$; CHR$(34)
PRINT #setfilenr%, CHR$(34); newfilename$; CHR$(34)
PRINT #setfilenr%, CHR$(34); outfilename$; CHR$(34)

FOR var = 1 TO novar
 PRINT #setfilenr%, CHR$(34); eo$(var); CHR$(34)
NEXT var

CALL closefilenr(filecount%, setfilenr%)
CLOSE setfilenr%
mes$ = "Save settings file " + setfilename$ + " : Ready"
RETURN

1335 ' read STIRLSYS.TXT file
filename$ = "STIRLSYS.TXT"
CALL existfile(filename$, filecount%, mes$)

CALL givefilenr(filecount%, filenr%)

IF mes$ <> "" THEN
  OPEN filename$ FOR OUTPUT AS filenr%
  setfilename$ = "STIRLSET.TXT"
  PRINT #filenr%, CHR$(34); setfilename$; CHR$(34)
ELSE
  OPEN filename$ FOR INPUT AS filenr%
  INPUT #filenr%, setfilename$
END IF
CLOSE filenr%

filename$ = setfilename$
CALL existfile(filename$, filecount%, mes$)

IF mes$ <> "" THEN
  ss$ = "  *  *     "
  IP$ = "  -    ----  -    -    -   -  - ---  --- -"
  op$ = "*       *              "
  inifilename$ = "STIRLSET.TXT"
  inpfilename$ = "STIRLINP.TXT"
  newfilename$ = "STIRLINP.TXT"
  outfilename$ = "STIRLOUT.TXT"
  eo$(1) = "100"
  eo$(9) = "100000"
  GOSUB 1330 'save settings file
ELSE
  GOSUB 1340 'read settings file
END IF
RETURN

1340 'read settings file
CALL givefilenr(filecount%, inifilenr%)
LOCATE 10, 1
PRINT "reading settings file inifilenr%="; inifilenr%;
PRINT " filename$="; filename$
OPEN filename$ FOR INPUT AS inifilenr%
INPUT #inifilenr%, ss$
INPUT #inifilenr%, IP$
INPUT #inifilenr%, op$
INPUT #inifilenr%, inifilename$
INPUT #inifilenr%, setfilename$
INPUT #inifilenr%, inpfilename$
INPUT #inifilenr%, newfilename$
INPUT #inifilenr%, outfilename$
FOR var = 1 TO novar
 INPUT #inifilenr%, eo$(var)
NEXT var

CALL closefilenr(filecount%, inifilenr%)
CLOSE inifilenr%
mes$ = "Reading settings file " + filename$ + " : Ready"
RETURN

2000 'input variable names
DATA "N%","DC","EC","RC","CR"
DATA "AL","DD","G","LB","SC","SE"
DATA "NH%","IH","LH","LI"
DATA "NR%","DR","LR","SR","DW","FF"
DATA "NC%","IC","LC","LD"
DATA "NE%","ID","LE"
DATA "SP","PS","OG%","FX","ME","FE"
DATA "L1","TF","TY","ZZ%","ZH"
DATA "KM","BF","MS%"
DATA "END"

2002 'inputs abbreviations
DATA "CYLINDERS #","CYLINDER DIAM.cm","PIST.END CL.cm"
DATA "CRANK RADIUS cm","CON.ROD LENGTH cm"
DATA "CRANK ANGLE degr","PIST.D.ROD DIAM.cm"
DATA "HOT CAP GAP cm","HOT CAP LENGTH cm","HOT CAP THICK.cm"
DATA "CYL.WALL THICK.cm"
DATA "HEATER TUBES #"
DATA "HEAT T.I.DIAM.cm","HEAT T. LENGTH cm","HEAT T. HTL.cm"
DATA "REGENERATORS #","REGENR.DIAM.cm"
DATA "REGENER.LENGTH cm","REG.WALL THICK.cm"
DATA "REG.WIRE DIAM.cm","REGENR.FILL FACTOR"
DATA "COOLER TUBES #"
DATA "COOL.T.I.DIAM.cm","COOL.T. LENGTH cm","COOL.T. HTL.cm"
DATA "CONNECTING DUCTS #","CON.DUCT I.DIAM.cm","CON.DUCT LENGTH cm"
DATA "SPEED rpm","AV.PRESSURE psia"
DATA "GAS H=1/HE=2/O=3 sw","WATER FLOW gpm","MECH.EFFICIENCY %"
DATA "FURN.EFFICIENCY %"
DATA "GAS LEAK./MPa.p.sec","HEAT.TUBE TEMP. F","COOL.WATER TEMP. F"
DATA "STAT.HEAT C.LOSS sw","STAT.HEAT C.LOSS wt"
DATA "METAL T.COND.w/cm k","BUGGER FACTOR","MODE SWITCH sw"

2006 'init load inputs abbreviationsni$()
RESTORE 2000
FOR var = 1 TO nivar
 READ ni$(var)
' PRINT var; " ="; ni$(var)
NEXT var
RETURN

2008 'list input
IF ipselect = 0 THEN ipselect = 1

CLS
PRINT "< INPUTS SCREEN > Type number or + or - or tab to select.";
PRINT " *=toggle optimize"
PRINT " = =change value. ?=info. Esc=outputs screen Enter=start program"
IF mes$ = "" THEN
PRINT "  ____________________________________________________________________________"
ELSE
  PRINT " MESSAGE: ===<<  "; mes$; "  >>==="
  mes$ = ""
END IF

FOR var = 1 TO nivar
' PRINT "var="; var; " var\2+1="; var \ 2
 IF var < nivar / 2 + 1 THEN
   LOCATE (3 + 1 + var), 1
   IF ipselect = var THEN
     PRINT "->";
   ELSE
     PRINT "  ";
   END IF
   PRINT MID$(IP$, var, 1);
   IF (i$ >= "0") AND (i$ <= "9") THEN
     IF var < 10 THEN PRINT " ";
     PRINT var;
   ELSE
     PRINT " "; ni$(var);
   END IF

   LOCATE (3 + 1 + var), 8
   RESTORE 2002
   FOR abbrev = 1 TO var
     READ abbrev$
   NEXT abbrev
   PRINT " "; abbrev$;
   LOCATE (3 + 1 + var), 27
   PRINT vi(var);
 ELSE
   LOCATE (4 + var - nivar \ 2), 40
   IF ipselect = var THEN
     PRINT "->";
   ELSE
     PRINT "  ";
   END IF
   PRINT MID$(IP$, var, 1);
   IF (i$ >= "0") AND (i$ <= "9") THEN
     PRINT var;
   ELSE
     PRINT " "; ni$(var);
   END IF
   LOCATE (4 + var - nivar \ 2), 47
   RESTORE 2002
   FOR abbrev = 1 TO var
     READ abbrev$
   NEXT abbrev
   PRINT " "; abbrev$;
   LOCATE (4 + var - nivar \ 2), 67
   PRINT vi(var);
 END IF
NEXT var

2009 i$ = INKEY$
IF i$ = "" THEN GOTO 2009
IF i$ = CHR$(27) THEN
  GOTO 2011'outputs screen
ELSEIF i$ >= "0" AND i$ <= "9" THEN
  IF VAL(i$) + 10 * ipselect <= nivar THEN
    ipselect = VAL(i$) + 10 * ipselect
    GOTO 2008
  ELSE
    ipselect = VAL(i$)
    GOTO 2008
  END IF
ELSEIF (i$ = "*") OR i$ = " " THEN
  IF MID$(IP$, ipselect, 1) = " " THEN
    MID$(IP$, ipselect, 1) = "*"
  ELSEIF MID$(IP$, ipselect, 1) = "*" THEN
    MID$(IP$, ipselect, 1) = " "
  ELSE
    mes$ = "## ERROR, NON OPTIMIZABLE INPUT ##"
  END IF
  GOTO 2008
ELSEIF i$ = "+" OR (i$ = CHR$(0) + CHR$(80)) THEN
  ipselect = ipselect + 1
  IF ipselect > nivar THEN ipselect = 1
  GOTO 2008
ELSEIF i$ = "-" OR (i$ = CHR$(0) + CHR$(72)) THEN
  ipselect = ipselect - 1
  IF ipselect < 1 THEN ipselect = nivar
  GOTO 2008
ELSEIF ASC(i$) = 9 OR (i$ = CHR$(0) + CHR$(75)) OR (i$ = CHR$(0) + CHR$(77)) THEN
  IF ipselect < nivar / 2 + 1 THEN
    IF i$ = CHR$(0) + CHR$(75) THEN
      GOTO 1100'start screen
    ELSE
      ipselect = ipselect + nivar \ 2
      IF ipselect > nivar THEN ipselect = nivar
    END IF
  ELSE
    IF (i$ = CHR$(0) + CHR$(77)) OR (ASC(i$) = 9) THEN
      GOTO 2011'outputs screen
    ELSE
      ipselect = ipselect - nivar \ 2
    END IF
  END IF
  GOTO 2008
ELSEIF i$ = CHR$(13) THEN ' Enter
  CLS
  do$ = "calculate"
  RETURN
ELSEIF i$ = "=" THEN
  LOCATE 3, 1
  PRINT ">>>>> INPUT NEW VALUE >> ";
  INPUT vi(ipselect)

  var = ipselect
  var$ = ni$(ipselect)
  GOSUB 2015 'load vi(var) in input value
  LOCATE 3, 1
  PRINT "                                                                              ";
  GOTO 2008
ELSEIF i$ = "?" THEN
  CALL givefilenr(filecount%, varfilenr%)
  varfile$ = "STIRLD.HLP"
  PRINT "varfilenr%="; varfilenr%, "varfile$="; varfile$
  OPEN varfile$ FOR INPUT AS varfilenr%
  CLS
  PRINT
  line$ = ""
  WHILE line$ <> "start"
   INPUT #varfilenr%, line$
  WEND
  WHILE line$ <> ni$(ipselect)
   INPUT #varfilenr%, line$
  WEND
  PRINT "INFORMATION ABOUT ";
  RESTORE 2002
  FOR abbrev = 1 TO ipselect
   READ abbrev$
  NEXT abbrev
  PRINT abbrev$
  PRINT "================================================================="
  PRINT "variable name in program ="; line$
  INPUT #varfilenr%, line$
  PRINT "value="; vi(ipselect); " "; line$
  IF ni$(ipselect) <> "MS%" THEN
    RESTORE 2045 'initial input variables
    var$ = ""
    varname = 0
    WHILE var$ <> ni$(ipselect)
     varname = varname + 1
     READ var$
    WEND
  
    RESTORE 2040 'initial input variable values
    FOR var = 1 TO varname
     READ var$
    NEXT var
    READ varval
    PRINT "was   "; varval; " as initial value in this program for GM4L23"
  END IF
  PRINT
  INPUT #varfilenr%, line$
  PRINT line$
  CALL closefilenr(filecount%, varfilenr%)
  PRINT
  PRINT "<<< HIT ANY KEY TO CONTINUE >>>"; i$
2010 i$ = INKEY$
  IF i$ = "" THEN GOTO 2010
  GOTO 2008
END IF
GOTO 2008

2011 'outputs screen
GOSUB 2034 'load output values to vo()
CLS
IF opselect = 0 THEN opselect = 1
IF mes$ <> "" THEN
  PRINT mes$
  PRINT
  PRINT "The outputs are not reliable."
  PRINT
  PRINT "Change the inputs."
  skipit = false
  GOSUB 1188 'Hit key
  CLS
END IF

PRINT "< OUTPUTS SCREEN > Type number or + or - or tab to select.";
PRINT " *=toggle optimize"
PRINT " = =change optimum. ?=info. Esc=optimizer screen, Enter=calculate"
FOR var = 1 TO novar
   LOCATE (2 + var), 1
   IF opselect = var THEN
     PRINT "->";
   ELSE
     PRINT "  ";
   END IF
   PRINT MID$(op$, var, 1);
   IF var < 10 THEN PRINT " ";
   PRINT var;
   PRINT no$(var);
   LOCATE (2 + var), 11
   RESTORE 2038 'abbreviations output variables
   FOR abbrev = 1 TO var
     READ abbrev$
   NEXT abbrev
   PRINT " "; abbrev$;
   LOCATE (2 + var), 41
   IF INSTR(STR$(vo(var)), ".") THEN
     PRINT LEFT$("             ", 10 - INSTR(STR$(vo(var)), "."));
     PRINT vo(var);
   ELSE
     PRINT RIGHT$("            " + STR$(vo(var)), 9);
   END IF
   LOCATE (2 + var), 64
   PRINT eo$(var);
NEXT var

2013 i$ = INKEY$
IF i$ = "" THEN GOTO 2013

IF i$ = CHR$(27) THEN
  GOTO 1241 'optimizer screen
ELSEIF ASC(i$) = 9 OR (i$ = CHR$(0) + CHR$(77)) THEN
  spacebarstep = true
  GOTO 1241 'optimizer screen
ELSEIF (i$ = CHR$(0) + CHR$(75)) THEN
  GOTO 2008'inputs screen
ELSEIF i$ >= "0" AND i$ <= "9" THEN
  IF VAL(i$) + 10 * opselect <= novar THEN
    opselect = VAL(i$) + 10 * opselect
    GOTO 2011'outputs screen
  ELSE
    opselect = VAL(i$)
    GOTO 2011'outputs screen
  END IF
ELSEIF i$ = "*" OR i$ = " " THEN
  IF MID$(op$, opselect, 1) = " " THEN
    MID$(op$, opselect, 1) = "*"
  ELSE
    MID$(op$, opselect, 1) = " "
  END IF
  GOTO 2011'outputs screen
ELSEIF i$ = "+" OR (i$ = CHR$(0) + CHR$(80)) THEN
  opselect = opselect + 1
  IF opselect > novar THEN opselect = 1
  GOTO 2011'outputs screen
ELSEIF i$ = "-" OR (i$ = CHR$(0) + CHR$(72)) THEN
  opselect = opselect - 1
  IF opselect < 1 THEN opselect = novar
  GOTO 2011
ELSEIF i$ = CHR$(13) THEN 'Enter
  CLS
  do$ = "calculate"
  RETURN
ELSEIF i$ = "=" THEN
  LOCATE 2, 1
  PRINT "                                                                              ";
  LOCATE 2, 1
  PRINT ">>>>> INPUT NEW OPTIMUM VALUE >>";
  INPUT eo$(opselect)
  GOTO 2011
ELSEIF i$ = "?" THEN
  CALL givefilenr(filecount%, varfilenr%)
  varfile$ = "STIRLD.HLP"
  PRINT "varfilenr%="; varfilenr%, "varfile$="; varfile$
  OPEN varfile$ FOR INPUT AS varfilenr%
  CLS
  PRINT
  line$ = ""
  WHILE line$ <> "start"
   INPUT #varfilenr%, line$
  WEND
  WHILE line$ <> no$(opselect)
   INPUT #varfilenr%, line$
  WEND
  PRINT "INFORMATION ABOUT ";
  RESTORE 2038 'abbreviations output variables
  FOR abbrev = 1 TO opselect
   READ abbrev$
  NEXT abbrev
  PRINT abbrev$
  PRINT "================================================================="
  PRINT "variable name in program ="; line$
  INPUT #varfilenr%, line$
  PRINT "value="; vo(opselect); " "; line$

  INPUT #varfilenr%, line$
  PRINT line$
  CALL closefilenr(filecount%, varfilenr%)
  PRINT
  PRINT "<<< HIT ANY KEY TO CONTINUE >>>"; i$
2014 i$ = INKEY$
  IF i$ = "" THEN GOTO 2014
  GOTO 2011
END IF
RETURN

2015 'load vi(var) in input value
IF var$ = "DC" THEN DC = vi(var): RETURN
IF var$ = "LC" THEN LC = vi(var): RETURN
IF var$ = "LD" THEN LD = vi(var): RETURN
IF var$ = "IC" THEN IC = vi(var): RETURN
IF var$ = "NC%" THEN NC% = vi(var): RETURN
IF var$ = "DW" THEN DW = vi(var): RETURN
IF var$ = "FX" THEN FX = vi(var): RETURN
IF var$ = "NE%" THEN NE% = vi(var): RETURN
IF var$ = "FE" THEN FE = vi(var): RETURN
IF var$ = "OG%" THEN OG% = vi(var): RETURN
IF var$ = "ZZ%" THEN ZZ% = vi(var): RETURN
IF var$ = "ZH" THEN ZH = vi(var): RETURN
IF var$ = "IH" THEN IH = vi(var): RETURN
IF var$ = "LH" THEN LH = vi(var): RETURN
IF var$ = "LI" THEN LI = vi(var): RETURN
IF var$ = "NH%" THEN NH% = vi(var): RETURN
IF var$ = "DD" THEN DD = vi(var): RETURN
IF var$ = "G" THEN G = vi(var): RETURN
IF var$ = "LB" THEN LB = vi(var): RETURN
IF var$ = "PS" THEN PS = vi(var): RETURN
IF var$ = "KM" THEN KM = vi(var): RETURN
IF var$ = "SC" THEN SC = vi(var): RETURN
IF var$ = "SE" THEN SE = vi(var): RETURN
IF var$ = "DR" THEN DR = vi(var): RETURN
IF var$ = "LR" THEN LR = vi(var): RETURN
IF var$ = "SR" THEN SR = vi(var): RETURN
IF var$ = "NR%" THEN NR% = vi(var): RETURN
IF var$ = "FF" THEN FF = vi(var): RETURN
IF var$ = "CR" THEN CR = vi(var): RETURN
IF var$ = "RC" THEN RC = vi(var): RETURN
IF var$ = "N%" THEN N% = vi(var): RETURN
IF var$ = "AL" THEN AL = vi(var): RETURN
IF var$ = "TF" THEN TF = vi(var): RETURN
IF var$ = "TY" THEN TY = vi(var): RETURN
IF var$ = "SP" THEN SP = vi(var): RETURN
IF var$ = "ID" THEN ID = vi(var): RETURN
IF var$ = "LE" THEN LE = vi(var): RETURN
IF var$ = "NE%" THEN NE% = vi(var): RETURN
IF var$ = "BF" THEN BF = vi(var): RETURN
IF var$ = "EC" THEN EC = vi(var): RETURN
IF var$ = "L1" THEN L1 = vi(var): RETURN
IF var$ = "MS%" THEN MS% = vi(var): RETURN'new to switch program changes on
RETURN

2022 'load input values in vi()
vi(1) = N%
vi(2) = DC
vi(3) = EC
vi(4) = RC
vi(5) = CR
vi(6) = AL
vi(7) = DD
vi(8) = G
vi(9) = LB
vi(10) = SC
vi(11) = SE
vi(12) = NH%
vi(13) = IH
vi(14) = LH
vi(15) = LI
vi(16) = NR%
vi(17) = DR
vi(18) = LR
vi(19) = SR
vi(20) = DW
vi(21) = FF
vi(22) = NC%
vi(23) = IC
vi(24) = LC
vi(25) = LD
vi(26) = NE%
vi(27) = ID
vi(28) = LE
vi(29) = SP
vi(30) = PS
vi(31) = OG%
vi(32) = FX
vi(33) = ME
vi(34) = FE
vi(35) = L1
vi(36) = TF
vi(37) = TY
vi(38) = ZZ%
vi(39) = ZH
vi(40) = KM
vi(41) = BF
vi(42) = MS%
RETURN


2032 'output var names
DATA "B","D","BP","HW","RW"
DATA "CF","IP","MF","NP","BH"
DATA "RH","QS","QP","SL","ZH"
DATA "A","QN","C1","E","PX"
DATA "PN","TH","TC"
DATA "END"

2034 'load output values to vo()
RESTORE 2032
FOR var = 1 TO novar
 READ var$
 GOSUB 2036 'load value output variable in vo(var)
NEXT var
RETURN

2035 'init load output vars names no$()
RESTORE 2032
FOR var = 1 TO novar
 READ no$(var)
 NEXT var
RETURN

2036 'load value output variable in vo(var)
IF var$ = "B" THEN vo(var) = B: RETURN
IF var$ = "D" THEN vo(var) = D: RETURN
IF var$ = "BP" THEN vo(var) = BP: RETURN
IF var$ = "HW" THEN vo(var) = HW: RETURN
IF var$ = "RW" THEN vo(var) = RW: RETURN
IF var$ = "CF" THEN vo(var) = CF: RETURN
IF var$ = "IP" THEN vo(var) = IP: RETURN
IF var$ = "MF" THEN vo(var) = MF: RETURN
IF var$ = "NP" THEN vo(var) = NP: RETURN
IF var$ = "BH" THEN vo(var) = BH: RETURN
IF var$ = "RH" THEN vo(var) = RH: RETURN
IF var$ = "QS" THEN vo(var) = QS: RETURN
IF var$ = "QP" THEN vo(var) = QP: RETURN
IF var$ = "SL" THEN vo(var) = SL: RETURN
IF var$ = "ZH" THEN vo(var) = ZH: RETURN
IF var$ = "A" THEN vo(var) = A: RETURN
IF var$ = "QN" THEN vo(var) = QN: RETURN
IF var$ = "C1" THEN vo(var) = C1: RETURN
IF var$ = "E" THEN vo(var) = E: RETURN
IF var$ = "PX" THEN vo(var) = PX: RETURN
IF var$ = "PN" THEN vo(var) = PN: RETURN
IF var$ = "TM" THEN vo(var) = TM: RETURN
IF var$ = "TW" THEN vo(var) = TW: RETURN
IF var$ = "TH" THEN vo(var) = TH: RETURN
IF var$ = "TC" THEN vo(var) = TC: RETURN
RETURN

2037 'load value in output variable vo()
vo(1) = B
vo(2) = D
vo(3) = BP
vo(4) = HW
vo(5) = RW
vo(6) = CF
vo(7) = IP
vo(8) = MF
vo(9) = NP
vo(10) = BH
vo(11) = RH
vo(12) = QS
vo(13) = QP
vo(14) = SL
vo(15) = ZH
vo(16) = A
vo(17) = QN
vo(18) = C1
vo(19) = E
vo(20) = PX
vo(21) = PN
vo(22) = TH
vo(23) = TC
RETURN

2038 'abbreviations output variables
DATA "INDICATED EFFICIENCY %","OVERALL EFFICIENCY %"
DATA "BASIC POWER watts","HEATER WINDAGE LOSS watts"
DATA "REGENERATOR FLOW LOSS watts","COOLER WINDAGE LOSS watts"
DATA "NET POWER watts","MECHANICAL FRICTTION watts"
DATA "BRAKE POWER watts","BASIC HEAT REQUIRMENT watts"
DATA "REHEAT LOSS watts","SHUTTLE LOSS watts"
DATA "PUMPING LOSS watts","TEMPERATURE SWING LOSS watts"
DATA "CONDUCTION LOSS watts","FLOW FRICTION CREDIT watts"
DATA "HEAT TO ENGINE watts","FURNACE LOSS watts"
DATA "FUEL INPUT watts","MAXIMUM PRESSURE MPa","MINIMUM PRESSURE MPa"
DATA "EFFECTIVE HOT SPACE TEMP.k","EFFECTIVE COLD SPACE TEMP.k"

2040 'input exemple data
DATA ".null."
DATA 10.16 ,12.9 ,12.02, .115, .167, 312
DATA  .00432, 25.0, 90, 80, 1, 0
DATA 9680, 41.8, 25.58, .472, .640, 36, 4.06
DATA .0174533, .0406, 6.40, 1400, .2, .0635, .1016
DATA .0510, 2.500, 3.500, 6, .2, 13.65, 2.325
DATA 4, 90, 1200, 135, 2000,. 435, .1532
DATA 0.76, 71, 6, .4, .006894, 2.54, 60
DATA 8.314, 1.341e-3, .0406, 0
DATA 0
DATA bottom
'0=MS%=Program mode switch     , 179 =AS (not used)

2045 'initial input variables
DATA DC, LC, LD, IC, OC, NC%
DATA DW, FX, ME, FE, OG%, ZZ%
DATA ZH, LH, LI, IH, OH, NH%, DD
DATA RA, G, LB, PS, KM, SC, SE
DATA SR, LR, DR, NR%, FF, CR, RC
DATA N%, AL, TF, TY, SP, AA, BA
DATA ID, LE, NE%, BF, PP, CM, FQ
DATA R, HP, EC, L1, ASAS
DATA MS%

2050 'read exemple values
RESTORE 2040
READ i$
READ DC, LC, LD, IC, OC, NC%
READ DW, FX, ME, FE, OG%, ZZ%
READ ZH, LH, LI, IH, OH, NH%, DD
READ RA, G, LB, PS, KM, SC, SE
READ SR, LR, DR, NR%, FF, CR, RC
READ N%, AL, TF, TY, SP, AA, BA
READ ID, LE, NE%, BF, PP, CM, FQ
READ r, HP, EC, L1, ASAS

PRINT
PRINT "  Original Martini exemple input values:"
PRINT "----------------------------------------------------"
PRINT DC; ","; LC; ","; LD; ","; IC; ","; OC; ","; NC%; ","; PI
PRINT P4; ","; DW; ","; FX; ","; ME; ","; FE; ","; OG%; ","; ZZ%
PRINT ZH; ","; LH; ","; LI; ","; IH; ","; OH; ","; NH%; ","; DD
PRINT RA; ","; G; ","; LB; ","; PS; ","; KM; ","; SC; ","; SE
PRINT SR; ","; LR; ","; DR; ","; NR%; ","; FF; ","; CR; ","; RC
PRINT N%; ","; AL; ","; TF; ","; TY; ","; SP; ","; AA; ","; BA
PRINT ID; ","; LE; ","; NE%; ","; BF; ","; PP; ","; CM; ","; FQ
PRINT r; ","; HP; ","; EC; ","; L1; ","; ASAS
PRINT
RETURN

'
SUB closefilenr (filecount%, filenr%)
'call closefilenr(filecount%,filenr%)

IF filecount% AND 2 ^ filenr% = 0 THEN
  PRINT "<< ERROR >> filenr"; filenr%; " is already closed."
  PRINT "Hit the 'Enter' key to continue.";
  INPUT i$
ELSE
  CLOSE #filenr%
  filecount% = filecount% - 2 ^ filenr%
END IF

END SUB

SUB dosfilename (filename$, mes$)
mes$ = ""
IF filename$ = "" THEN
  mes$ = "NO FILENAME GIVEN, FILE NAME NOW IS THE DEFAULD ONE"
  EXIT SUB
END IF
WHILE LEFT$(filename$, 1) = " "
 filename$ = RIGHT$(filename$, LEN(filename$) - 1)
WEND
WHILE RIGHT$(filename$, 1) = " "
  filename$ = LEFT$(filename$, LEN(filename$) - 1)
WEND
IF INSTR(filename$, ".") = 0 THEN
  IF LEN(filename$) > 8 THEN
    mes$ = "FILE NAME IS TO LONG"
    EXIT SUB
  END IF
  FOR ptr = 1 TO LEN(filename$)
   IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
     MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
   ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
     IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
       mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
       mes$ = mes$ + " IN FILE NAME"
       EXIT SUB
     END IF
   END IF
  NEXT ptr
ELSE
  IF INSTR(filename$, ".") > 9 THEN
    mes$ = "FILE NAME IS TO LONG"
    EXIT SUB
  ELSEIF INSTR(filename$, ".") = 1 THEN
    mes$ = "WRONG FILE NAME, STARTING WITH CHARACTER '.'"
    EXIT SUB
  ELSEIF INSTR(filename$, ".") = LEN(filename$) THEN
    mes$ = "WRONG FILE NAME, ENDING WITH CHARACTER '.'"
    EXIT SUB
  ELSEIF LEN(filename$) - INSTR(filename$, ".") > 3 THEN
    mes$ = "FILE NAME EXTESION IS TO LONG."
    EXIT SUB
  END IF
  FOR ptr = 1 TO INSTR(filename$, ".") - 1
   IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
     MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
   ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
     IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
       mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
       mes$ = mes$ + " IN FILE NAME."
       EXIT SUB
     END IF
   END IF
  NEXT ptr
  FOR ptr = INSTR(filename$, ".") + 1 TO LEN(filename$)
   IF (MID$(filename$, ptr, 1) >= "a") AND (MID$(filename$, ptr, 1) <= "z") THEN
     MID$(filename$, ptr, 1) = CHR$(ASC(MID$(filename$, ptr, 1)) - 32)
   ELSEIF (MID$(filename$, ptr, 1) < "A") OR (MID$(filename$, ptr, 1) > "Z") THEN
     IF INSTR("1234567890(){}@#$%^&!-_'/", MID$(filename$, ptr, 1)) = 0 THEN
       mes$ = "WRONG CHARACTER " + CHR$(32) + MID$(filename$, ptr, 1) + CHR$(32)
       mes$ = mes$ + " IN FILENAME EXTENSION."
       EXIT SUB
     END IF
   END IF
  NEXT ptr
END IF

END SUB

'
SUB existfile (filename$, filecount%, mes$)
'call existfile(filename$,filecount%,mes$)
mes$ = ""
CALL givefilenr(filecount%, filenr%)
OPEN filename$ FOR APPEND AS filenr%
fileend = SEEK(filenr%)
CLOSE filenr%
IF fileend = 1 THEN
  KILL filename$
  mes$ = "file " + filename$ + " : does not exist."
END IF
END SUB

SUB givefilenr (filecount%, filenr%)
'call givefilenr(filecount%,filenr%)
IF filecount% = 0 THEN filecount% = 1
filenr% = 0
WHILE (filecount% OR 2 ^ filenr%) = filecount%
 filenr% = filenr% + 1
WEND
filecount% = filecount% + 2 ^ filenr%

END SUB

SUB makestr (printv, printv$)
'call makestr (printv, printv$)
IF INSTR(STR$(printv), ".") THEN
  printv$ = LEFT$("             ", 9 - INSTR(STR$(printv), "."))
  printv$ = printv$ + STR$(printv)
ELSE
  printv$ = RIGHT$("            " + STR$(printv), 8)
END IF
printv$ = LEFT$(printv$ + "       ", 13)
END SUB

'==================================================================

DO NOT COPY THIS !!!
Here comes the help file.
Copy the text between the ############# lines and save
it with the name STRLD.HLP in the same directory.
'##################################################################
"start"
"START1"
"START CALCULATIONS"
"In case you hit the 'Enter' key when this option is selected,
the program starts calculating with the inputs you have chosen.
When this option is checked with '*' and the settings are
saved in the settings file, the program skips this start screen
and begins calculating with the checked inputs option on 
tart up. In case the program hangs because of wrong inputs you
can delete the STIRL.SYS file, and restart the program."
"START2"
"skip extra information on screen"
"This option skips the screen output during the calculations,
so you don't have to hit the 'Enter' key to resume the
program execution every time the information is displayed."
"START3"
"start optimizer in single step mode (step with space bar)"
"This option makes the optimizer stop after every time
new outputs are generated. Every time you hit the space bar,
the optimizer makes one step. Hitting an other key makes 
the optimizer run until you hit the space bar or until the
optimizer has found the optimum values for inputs you have
checked with '*' in the , optimized to the
optimum values of the checked outputs in the ."
"START4"
"Optimize on all records of the loaded input file switch."
"In case there are more than one records in the input file, 
the optimizer optimizes for all the cases that are specified
in the records of the input file.
When this switch is off, or when there is only one record,
only the curent loaded inputs are considerdby the optimizer.
The records must start with a ',null.' field, or when there is
a 'start' field, the next field specifies the number of inputs
per record. The next fields are the variable names. The number
of records counted till the end of the file, or an 'end' field."
"START5"
"load original Martini exemple GM4L23 inputs"
"These inputs are listed in the Stirling Design Manual,
written by William R. Martini. Those are the values for the
Siemens type Stirling engine GM4L23 of General Motors."
"START6"
"load changed GM4L23 inputs"
"Changed inputs are:
    LB = 10.03 (hot cap length in cm.)
    SC = .4305 (hot cap thickness in cm.)
    SE = 1.25 (cylinder wall thickness in cm.)
    SR = .41 (regenerator wall thickness in cm.)
    LR = 2.79 (regenerator length in cm.)
    MS = 255 (mode switch to change calculations)
"
"START7"
"load settings from file"
"Here you can load settings you have previously saved with the
next option. When checked and saved as settings file,the settings
are automaticaly loaded on start up. In case you hit 'Enter' when
this option is selected, the settings are loaded from the named
settings file. With '=' you can change the settings file. When 
checked with '*' these settings are automaticaly loaded on start up."
"START8"
"save settings to file"
"In case you hit the 'Enter' key when this option is selected,
you save the current settings to the named settings file.
When this option is checked the settings are automaticaly saved when
you start the calculations, so your last settings are preserved.
In case this name is the same as the file name in the option
'load settings from file' and that option is checked with '*',
these settings are automaticaly loaded on start up."
"START9"
"load inputs from file"
"In case you have saved the inputs with the next option, you can load
them again from the named file by hitting the 'Enter' key when this
option is selected. When checked and the settings are saved with the
'save settings to file' option, and the 'load settings from file'
option is checked, whith this input file name, the inputs are automaticaly
loaded on start up."
"START10"
"save inputs to file"
"With the 'Enter' key you can save the inputs to the named file.
You can change the 'new inputs file name' by hitting the '=' key
when this option is selected.
In case this option is checked with '*' the last inputs are automaticaly
saved before the calculations start in the named new inputs file name."
"START11"
"save outputs in file"
"The outputs are saved in the named output file. You can change that name
by pressing the '=' key.
Pressing the 'Enter' key when this item is selected saves the current in and 
outputs to the file. Before the start of the calculations most outputs are zero.
When this item is checked with a '*' the outputs are automaticaly saved during 
the execution of the program.
You can read the output file with a text editor (like notepad or wordpad) in the
directory were the basic file is."
"N%"
"number of similar engine units per engine."
"Must be a whole number, 1 or more. It is the actual number of cylinders in a siemens type. With more units it is asif more engines are working together. It makes the engine run smoother."
"DC"
"cm, diameter cylinder."
"Assumed is that all cylinders have the same diameter."
"EC"
"cm, piston end clearance."
"Piston clearance at the hot side. In a gamma it is also the end clearance of the displacer at the cold side. When the end clearance in the second cylinder is specified as zero, this value is also used there."
"RC"
"cm, crank radius cylinder 1."
"Used to calculate volume changes. Cranck is exspected to be in line with the cylinder."
"CR"
"cm, length of connecting rod 1."
"Used to calculate volume changes."
"DC"
"cm, diameter cylinder 2."
"When zero the program assumes that this cylinder has the same diameter as cylinder 1."
"EC2"
"cm, piston end clearance at cylinder 2 side."
"Piston clearance at expansion cylinder side."
"RC2"
"cm, crank radius 2."
"Used to calculate volume changes. Cranck is exspected to be in line with the cylinder."
"CR2"
"cm, length of connecting rod 2."
"Used to calculate volume changes."
"AL"
"crank angle between the crank on cylinder 1 and cyl. 2."
"This is the angle between the piston rod and the displacer rod on the crank shaft. It determins how long after the start of the heating up of the working gas the expansion starts."
"CYLAL"
"Angle between the two cylinders. Often 0 or 90 degrees."
"This is the angle between cylinder 1 (hot) and cylinder 2 (cold or expansion cylinder). Zero in a beta type. Determines together with the crank angle the phase angle between the two piston movements."
"DD"
"cm, diameter of piston/displacer drive rod (zero for alpha). Negative means the rod is on the hot side."
"A double acting machine (like a beta, a gamma or a siemens) has a piston drive rod. A single acting machine as an alpha has not. A low temperature gamma often has the displacer drive rod on the hot side. Indicate that with a negative value."
"G"
"cm, gap in hot cap."
"This is the annular gap between the hot cylinder side wall and the piston. The heat input area of a stirling engine is generally called the hot cap or hot end of the engine."
"LB"
"cm, length of hot cap."
"Used for gap in hot cap dead volume, pumping or appendix loss, shuttle heat loss and static heat loss. The heat input area of a stirling engine is generally called the hot cap or hot end of the engine."
"SC"
"cm, wall thickness of hot cap."
"Used for static heat loss. The hot cap conducts heat from the hot side to the cold side in case of a beta, a gamma or a siemens type engine. The heat input area of a stirling engine is generally called the hot cap or hot end of the engine."
"SE"
"cm, wall thickness of expansion cylinder wall."
"Used for shuttle heat loss. The expansion cylinder loses heat, due to conduction through the wall."
"NH%"
"number of heater tubes per cylinder."
"Must be more than zero."
"LH"
"cm, heater tube length."
"Used for heater tube dead volume and the gas heating windage loss."
"LI"
"cm, heater tube heat transfer length."
"The length of each tube that actualy is part of the heat transfer area of heater."
"IH"
"cm, inside diameter of heater tubes."
"Used for heat transfer area of heater, heater tube dead volume and the  heater windage loss."
"OH"
"cm, heater tube outside diameter (not used)."
"Not used"
"NR%"
"number of regenerators per cylinder."
"More than zero. Must be a whole number."
"SR"
"cm, wall thickness of regenerator housing."
"Used to calculate the static heat loss. Not used when static heat loss is specified."
"LR"
"cm, length of regenerator."
"Used to calculate the volume an capacity."
"DR"
"cm, diameter of regenerator."
"This is the inside diameter, used to calculate the volume and capacity."
"DW"
"cm, diameter of 'wire' in regenerator."
"Wire used as heat capacitor, .0017 inch *(2.54)=0.00432 cm"
"FF"
"filler factor, fraction of regenerator volume filled with solid."
"Filler factor .2 means 20% of the regenerator is filled with wire, 80% of the regenerator is filled with working gas."
"NC%"
"number of cooler tubes per cylinder."
"Cooler tubes are used to increase the area where the water cools the working gas."
"LC"
"cm, length of cooler tube."
"This is the length that defines the dead space inside the tube."
"LD"
"cm, heat transfer length of cooler tube."
"This is the length of the tube that defines the area for heat transfere from the working gas to the water."
"IC"
"cm, inside diameter of cooler tube."
"Defines cold dead space and heat transfer area for cooler."
"OC"
"cm, outside diameter of cooler tubes. (not used)."
"Not used. Can be used for heat conduction for cooler."
"NE%"
"number of connecting ducts per cylinder."
""
"ID"
"cm, inside diameter of connecting duct."
""
"LE"
"cm, length of connecting duct."
""
"SP"
"RPM, engine speed."
""
"PS"
"psia, average pressure."
""
"OG%"
"operating gas, hydrogen, helium or air."
""
"FX"
"cooling water flow @ 2000 rpm per cylinder"
"63.12 * FX= gram per second. This value is not corrected for other speeds."
"ME"
"mechanical efficiency, %."
""
"FE"
"furnace efficiency, %"
""
"L1"
"fraction of total gas charge leaking per MPa P per second."
""
"TEMP"
"temperature specified in: celcius, kelvin or fahrenheit."
""
"TF"
"Inside heater tube wall temperature, specified in degree fahrenheit"
""
"TY"
"inlet cooling water temperature, specified in degree fahrenheit"
""
"ZZ%"
"static heat conduction loss"
"Switch, Specified (below) OR calculated"
"ZH"
"watts, specified static heat conduction loss"
"You can specify the heat conduction loss here, but it is ignored when the switch ZZ is zero."
"KM"
"w/cm K, metal thermal conductivity."
"Needed for calculated heat conduction loss."
"BF"
"Bugger factor to convert power outputs to nearly what General Motors says they should be"
"This factor should idealy be one. In the program this factor is multiplied with 4 when the correction in the calculation for the heat and cool exchange area is switcht on with the mode switch (SW)"
"MS%"
"Program mode switch."
"When this variable is zero, the program works the same as the original program. With the exemple inputs for the GM4L23 it should give the outputs as given in   the Stirling engine design manual.
In case bit 1 =1 the heat transfer area for the cooler (AC) and the heater (AH) are not multiplied by the number of cylinders (N%). Where the bugger factor (BF) is used, it is multiplied by 4.
In case bit 2 =1 the volume of the connecting ducts is calculated with diameter (ID) to the power of 2. 
In case bit 3 =1 the cooling water flow is divided by 2000 and multiplied by the   engine speed (SP).
"
"BP"
"watts, basic power."
""
"BH"
"watts, sic heat requirement."
""
"HW"
"watts, heater windage loss."
""
"RH"
"watts, reheat loss"
"Reheat loss. One way that extra heat is required at the heat source is due to the inefficiency of the regenerator. The regenerator reheats the gas as it returns to the hot space. The reheat not supplied by the regenerator must be supplied by the heater as extra heat input."
"RW"
"watts, regenerator f.l."
"The regenerator flow loss is the windage loss due to the friction by the obstruction of the gas flow in the regenerator."
"QS"
"watts, shuttle loss"
"This is the heat loss due to the motion of the
displacer piston. In the stationary state, the
temperature of the piston will be approximate
ly equal to that of the adjacent cylinder
wall. When the piston moves, each section of
its surface moves to confront a new part of
the cylinder wall, at a different temperature.
Heat is transferred between the two adjacent
surfaces at different temperatures."
"CF"
"watts, cooler windage loss"
""
"QP"
"watts, pumping loss"
"pumping losses are aresult of the fact that there is a gap between the cylinder wall and the lower end of the piston where the sealing (piston ring) is, while the gap is open at hot side at the top. As the pressure in the cylinder varies, gas flows into and out of this volume. Sincethe lower end of the gap is colder, extra heat must be added to this gas as it leaves the space"
"IP"
"watts, net power"
"Net is the amount of work generated by a motor without taking into consideration the losses due to mechanical friction and any of the various auxiliary components that may slow down the actual speed of the motor."
"SL"
"watts, temperature swing"
""
"MF"
"watts, mechanical friction"
"Mechanical friction loss is mostly due to friction in seals and the bearings."
"ZH"
"watts, conduction loss"
""
"NP"
"watts, brake power"
"Brake is the amount of work generated by a motor without taking into consideration any of the various auxiliary components that may slow down the actual speed of the motor. Sometimes referred to as pure power, brake power is measured within the engine's output shaft."
"A"
"watts, flow friction credit"
"Flow friction losses are due to fluid friction when pressure drops let the gas mass move. Flow friction losses in heater and cooler are usualy small in comparison with the regenerator. The mass flow at the cold end is much more than the mass flow at the hot end, mostley due to the gas density change. Most of the flow friction is converted into heat. When this happens in the hot parts of the engine, the energy is not lost, but added again to the energy that is brought into the engine. (therefore a negative value as heat requirement)"
"QN"
"watts, heat to engine"
""
"B"
"%, indicated efficiency."
""
"C1"
"watts, furnace loss"
""
"D"
"%, overall efficiency"
""
"E"
"watts, fuel input"
""
"TM"
"hot metal temperature"
""
"TW"
"K. cooling water inlet temperature"
""
"TH"
"K. effective hot space temperature"
"The temperature in the hot space is lower than the hot metal temperature, because the gas flows from the cold space to the hot space. The heating of the gas in the hot space takes time (depended on the heat transfer area of the heater)."
"TC"
"K. effec.cold sp.temp."
"effective cold space temperature in degree Kelvin. This temperature is higher than the cooling water temerature, because the gas flows from the hot space to the cold space. The cooling of the gas in the cold space takes time (depended on the heat transfer area for the cooler)."
"end"

'######################################################################