'====================================================================
DECLARE FUNCTION xlarge! (X!(), NDIV%)
DECLARE FUNCTION small! (X!(), NDIV%)
DECLARE SUB makestr (printv!, length, decimals, printv$)
DECLARE SUB grafiek (opdracht$, x1!, y1!, x2!, y2!, kleur!)
DECLARE SUB closefilenr (filecount%, FILENUMMER%)
DECLARE SUB givefilenr (filecount%, FILENUMMER%)
DECLARE SUB pdint (X!, DMW!(), DMC!(), RVT!, DVC!(), NDIV%, DM!, PR!(), XINT!, DPR!(), XI1!, XI2!, XNHT!)
DECLARE SUB VOLW (W!(), WI!(), DW!(), DWI!(), CFI!, SFI!, ZZW!, NDIV%, SIFI!(), COFI!(), SALF!(), CALF!(), DALF!)
DECLARE SUB VOLC (DALF!, NF%, C!(), CI!(), DC!(), DCI!(), ZZC!, NDIV%, SIFI!(), COFI!(), SALF!(), CALF!(), PDR!, ALF!)
DECLARE FUNCTION plot! (X(), H!)

'www.scribd.com/doc/20191898/Stirling-Engine-Design-Manual
'        Rios 41.23 job (),kosdruk-notify=wrmjc*region=120k
'exec fortgclg
'fort.sysin dd *



'Expects number of ducts= number of regerators!
MS% = 1' mode switch for pumping loss Martini (off in Richard's exemple)

CLS


'AFC = Cooler free flow area, cm^2
'AFH = Heater free flow area, cm^2
'AFR = Regenerator free flow area, cm^2
'ALF = 4.7123889 (270 degrees)
'ARG = sin(PV angle)
'BDR = (input) Regenerator diameter, cm
'BEC = (input) Piston end clearence, cm
'BPD = (input) Piston diameter, cm
'BPL = (input) Hot cap length, cm
'BRC = (input) Piston gap, cm (of brg)
'BRL = (input) Regenerator length, cm
'BRO = (input) Regenerator density factor
'BST = (input) Piston stroke, cm
'BTC = (input) Effective cold temperature, K
'BTC1 = Cold metal temperature, K
'BTR = Regenerator temperature, K
'BTW = (input) Hot effective gas temperature, K
'BTW1 = Hot metal temperature, K
'BWD = (input) Effective regenerator wire diameter, cm
'C() = Cold space as fraction of the stroke amplitude at mid-increment
'      C() varies from 0 to 2 and back.
'CALF() = Sin change per radian increment
'CALFP = CALF()
'CFI = Cos of phase angle

'CI() = Same as C() for beginning of increment
'CMMAX = Largest cold dimensionless mass
'CMU = Cold hydrogen viscosity
'CNTU = Number of heat transfer units in cold space
'COFI() = Cos values for cold space
'CON = Conduction loss, watts
'COR = PMAX ^ (XNHT - 2) * DALF ^ (XNHT - 1) REGEL 340, 345, 346
'CP1 = Hydrogen heat capacity
'CRC = SQR(ZZC^2-CALF()^2)
'CRW = CRC in hot space
'CTD = (input) Cooler tube inside diameter, cm
'CTLL = (input) Total cooler tube length, cm
'CTLS = (input) Cooler cool tube length, cm
'CV1 = Hydrogen heat capacity
'DALF = 2phi/ndiv
'DC() = Angle derivative of C()
'DCI() = Angle derivative of CI()
'DDD = (input) Cooler duct diameter, cm
'DLL = (input) Cooler duct length, cm
'DM = Sum of changes in mass (DMRE)
'DMC = Cold dimension mass change XDMC()
'DMRE = Sum of changes in mass (DM)
'DMW = Hot dimensionless mass change XDMW()
'DMX = Dimensionless change in mass relating to X, the fraction from
'      the cold end
'DP = Change in pressure
'DPR() = DP array
'DTC = Cooler metal temperature - effective temperature
'DTW = HOT metal temperature - effective temperature
'DTH = Delta Th
'DV = Dead volume, cm^3
'DVC = DC()
'DVCI = DCI()
'DVW = DW()
'DVWI = DWI()
'DW()= Angle derivative of W()
'DWI() = Angle derivative of WI()

'DX = 1/xnds
'EX1 = 1-XNHT
'EX2 = 2-XNHT
'FC = Cold friction factor
'FFF = Friction flow credit, watts
'FH = Hot friction factor
'FR() = Regenerator friction factors (3 pts.)
'FI = Phase angle, rad.
'FI1 = Phase angle in deg.
'FIPV = pv angle (output) arcsin (ARG)
'FR() = Regenerator friction factor
'G1 = Y value subplot
'G2 = Y value subplot
'GDMS() = Calculated mass flow values
'GGV = Dead volume at side of hot cap, cm^3
'GINT() = Flow loss variable
'GI2() = Pressure drop value
'GI3() = Pressure drop value
'GLH = Heater pressure drop integral
'GLR = Regenerator pressure drop integral
'GLS = Cooler pressure drop integral
'H(1) = Fraction of total reduced dead volume from cold end to midway
'       in cooler.
'H(2) = Fraction of total reduced dead volume from the cold end through
'       the cooler
'H(3) = Fraction of total reduced dead volume from the cold end through
'       half the regenerator
'H(4) = Fraction of total reduced dead volume through de regenerator
'H(5) = Fraction of total reduced dead volume through the middle of the gas
'       heater (1-H(5) includes the rest of the heater and clearence on the
'       end and sides of the hot cap)
'HAC = Cold active volume amplitude, cm^3
'HAV = Hot active volume amplitude, cm^3
'HCV = Reduced cooler and cold ducting dead volume, dimensionless
'HEC = Reduced cold end clearence dead volume, dimensionless
'HGV = reduced hot cap gap dead volume, dimensionless
'HHC = reduced hot clearance dead volume, dimensionless
'HHV = reduced heater dead volume, dimensionless

'HMU= Hot hydrogen viscosity
'HNTU = .1044 * HTLS / (HTD * REH ^ .2)' regel 419,420 hydrogen ??
'HRV = Reduced regenerator dead volume, dimensionless
'HT = Basic heat input, watts
'HTD = (input) Heater tube inside diameter, cm
'HTE = Heat to engine, watts
'HTLL = (input) Total heater tube length, cm
'HTLS = (input) Heated heater tube length, cm
'HTW = Hot end heat transfer integral, dimensionless
'IND() = Array that shows if mass change is positive, or negative in
'        warm and cold sides
'J% = V [Temporary angle variable, radians
'K% = 1 if warm mass change is positive, 2 if negative
'L% = 1 if cold mass change is positive, 2 if negative
'LUP% = Iterational counter
'M% = X value for plot calculation
'MBR% = (input) Number of regenerators
'MCT% = (input) Number of cooler tubes per cylinder
'MHT% = (input) Number of heater tubes for cylinder
'MW% = Dimensionless mass in hot space = (mass, grams)(R)(BTW)/(PMX1(HAV))
'N% = NDIV% or X value for plot subroutine
'NN% = 1 up to phase angle, 2 after
'NDIV% = (input) Number of divisions per crank rotation (must be a multiple
'       of 4 so that the phase angle at 90 degrees can be an even number of
'       divisions)
'       (Programm must be revised if NDIV% is not 360)
'NDIV1% = NDIV% + 1
'NDS% = (input) Number of divisions in dead space
'NE% = ndiv/4 + 1
'NET% = (input) Regenerator filler option (+5 = metnet, -5 = screen)
'NF% = ndiv/4
'NFF% = NF%+1
'NFI% = (input) (phase angle)(NDIV%)/360
'NFIN% = Main loop final counter, for first part = phase angle, for second
'       part = end of cycle
'NIN% = NDS%+1
'NITE% = Cycle counter (counts to 15)
'NL% = (ndiv/2)+1
'NLOP% = Option counter limits changes in option to 7 (removed in final
'       version)
'0; = IND(K%,L%)-1,2,3 or 4, starts as 1
'NOC% = (input) number of cylinders
'NS% = (ndiv/4)+2
'NST% = Main loop initial counter, for the first part =1, for the second
'      part = phase angle
'NT% = (ndiv/4)+2
'NWR% = (input) Governs printout, zero for overall results only, different
'      from zero added pv data
'P = pressure, dimensionless
'PALF = Thermal diffusitivity of piston
'PDR = (input) Piston rod diameter, cm
'PI4 = phi/4 = .78539816
'PIE = 3.1415926#
'PAVG = Dimensionless average pressure
'PMAX = Maximum pressure, dimensionless
'PMIN = Minimum pressure, dimensionless
'PMX = Maximum pressure (MPa)
'PMX1 = (input) Avg. pressure MPa
'PR() = Pressure, dimensionless, fraction of maximum pressure
'PO = Basic power, watts
'POF = Net power, watts
'PS = Dimensiomless pressure from end of previous cycle
'PW = Pressure at halfway point for increment
'QB = Beta for shuttle heat loss calculation
'QCP = Cooler windage, watts
'QDK = Reheat factor
'QFS = Pumping loss factor
'QHC = Shuttle loss, watts
'QHG = Pumping loss, watts
'QHP = Heater windage, watts
'QHR = Reheat loss, watts
'QLM = Reheat factor
'QL1 = Shuttle factor
'QNPH = Reheat pressurization effect
'QNTU = Regenerator transfer units, dimensionless
'QP = Windage factor
'QR() = Regenerator windage loss values, watts
'QR1, QR2, QR3 line 409,410,411,412
'QRP = Average regenerator windage, watts
'R = Gas constant, joules/(gm)(K%)
'R2 = Constant = R(gc)2
'RE() = Regenerator Reynolds number in cold, middle, and hot part
'REC = Cold Reynolds number
'REH = Hot Reynolds number
'RER = Regenerator Reynolds number
'RMU = Regenerator hydrogen viscosity, g/cm sec.
'RNTU = Regenerator heat transfer units
'RP = Maximum pressure/minimum pressure
'RVT = Displacer mass ratio
'S = Pressure at halfway point, dimensionless
'SALF() = Sin values for cold space
'SALFP = Average sin values for cold space
'SFI = Sin of phase angle
'SHR = (input) Specific heat ratio for working gas
'SIFI() = SALF()
'SIFIP = SALF(1)
'SMC = Cold mass + 1/2 change in mass
'SMW = Hot mass + 1/2 change in mass
'SPD = (input) Engine speed, rad/sec
'TEC() = Dimensionless cold gas temperature
'TEST = Ensures that difference in dimensionless mass <.001
'TEST1 = Ensures that difference in dimensionless pressure <.005
'TEW() = Dimensionless hot gas temperature
'TMPC = Average TEC()
'TMPW = Average TEW()
'TCDM = Dimensionless average cold temperature for entire cycle
'TWDM = Dimensionless average warm temperature for entire cycle
'UDM() = Critical mass flow values from subplot
'UIN() = Critical pressure drop integral values from subplot
'UI23, 24, 33, 34 = Critical pressure drop values
'UPA = Power piston area, cm^2
'UTR = Temperature ratio = (hot metal temp, K%)/(cold metal temp, K%)
'VC = C()
'VCC = Cold volume cm^3
'VCD = Cold dead volume, cm^3
'VCI = CI()
'VD = Reduced dead volume, dimensionless
'VH = Hot volume, cm^3
'VHD = Hot dead volume, cm^3
'VRD = Regenerator dead volume, cm^3
'VT = Total volume, cm^3
'VW = W()
'VWI = WI()
'W() = Hot space as fraction of the stroke amplitude, calculated at mid
'      increment
'WC = Dimensionless cold work
'WI() = Same as W() for beginning of increment
'WMMAX = Largest hot dimensionless mass
'WW = Hot work, dimensionless
'X = short term variable
'XDMC() = Change in cold mass, grams
'XDMW() = Change in hot mass, grams
'XI1 = Pressure drop integral - account for the relationship between the
'      shapes of mass and pressure fluctuations
'XI2 = Influence of mass flow time variation on the heat transfer
'XI3 = xi1/xi2
'XINT = Basic pressure drop integral - for windage
'XMC = Cold gas mass, relative to total inventory
'XMCX() = Cold gas mass, grams
'XMT() = Total mass, grams
'XMW = Hot gas mass, relative to total inventory
'XMWS = Hot dimensionless gas mass from previous cycle
'XMWX() = Hot gas mass, grams
'XND = NDIV%
'XNDS = NDS%
'XNHT = (input) Value for exponent in heat transfer relation of regenerator
'       matrix
'XX = Short term variable
'Y =|DMX|

'ZEF = Indicated efficiency, %
'ZZC = (input) Connecting rod length/.5 stroke for cold piston
'ZZW = (input) Connecting rod length/.5 stroke for hot piston


'###########  I N I T ###############
 'D.3 Fortran Listing

'If memory is a problem all dimensions to 720 can be changed to 361
DIM XDMC(720), XDMW(720), IND(2, 2), PR(720), XMCX(720), XMWX(720)
DIM TEC(720), TEW(720), XMT(720)
DIM W(720), C(720), WI(720), CI(720), DW(720), DC(720), DWI(720)
DIM DCI(720), DPR(720)
DIM SIFI(720), COFI(720), SALF(720), CALF(720), H(5), UIN(5)
DIM UDM(5), GDMS(20), GINT(20), GI2(20), GI3(20), RE(3), FR(3), QR(3)

'constants are now setup
PIE = 4 * ATN(1)'PIE = 3.1415926#
PI4 = ATN(1)'PI4 = .78539816#

GOSUB 9900 'read variable text strings
'#####  MAIN  ###########

'output basic file
 CALL givefilenr(filecount%, OUTFILE%)
 bestand$ = "STIRLOUT.TXT"
 OPEN bestand$ FOR OUTPUT AS #OUTFILE%
 PRINT #OUTFILE%, "Stirling program output: "; bestand$
 GOSUB inputs 'load defaulds
 IF (NWR%) <> 0 THEN GOSUB printout

CALL closefilenr(filecount%, OUTFILE%)

INPUT "End of programm"; I$
END



inputs:
'Specified engine dimensions and operating conditions are read
'in from data cards

RESTORE 9999


READ name$
IF name$ = "end" OR name$ = "END" OR name$ = "" THEN RETURN
READ ZZC
'stops if next data is end, END or blank

IF (ZZC) = 0 THEN
  END ' end of program
END IF

'read(5.5) ZZC,ZZW,XNHT,PMX1,SPD,NET%,NOC%
'read (5+1)  SHR,NFI%,NDIV%,NWR%,NDS%,DDD,DLL
'read (5+1710) BTC,BTW,BPD,BST,BPL,BRC,BEC
'read (5+1720) PDR,BRO,BDR,BWD,BRL,MBR%,MCT%
'read (5+1730) CTD,CTLL,CTLS,HTD,HTLL,HTLS,MHT%

READ ZZW, XNHT, PMX1, SPD, NET%, NOC%
READ SHR, NFI%, NDIV%, NWR%, NDS%, DDD, DLL
READ BTC, BTW, BPD, BST, BPL, BRC, BEC
READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT%
READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT%
READ CON

PRINT ZZC, ZZC$
PRINT ZZW, ZZW$
PRINT XNHT, XNHT$
PRINT PMX1, PMX1$
PRINT SPD, SPD$
PRINT NET%, NET$
PRINT NOC%, NOC$
PRINT
PRINT SHR, SHR$
PRINT NFI%, NFI$
PRINT NDIV%, NDIV$
PRINT NWR%, NWR$
PRINT NDS%, NDS$
PRINT DDD, DDD$
PRINT DLL, DLL$
INPUT I$
PRINT BTC, BTC$
PRINT BTW, BTW$
PRINT BPD, BPD$
PRINT BST, BST$
PRINT BPL, BPL$
PRINT BRC, BRC$
PRINT BEC, BEC$
PRINT
'READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT%
PRINT PDR, PDR$
PRINT BRO, BRO$
PRINT BDR, BDR$
PRINT BWD, BWD$
PRINT BRL, BRL$
PRINT MBR%, MBR$
PRINT MCT%, MCT$
INPUT I$
'READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT%
PRINT CTD, CTD$
PRINT CTLL, CTLL$
PRINT CTLS, CTLS$
PRINT HTD, HTD$
PRINT HTLL, HTLL$
PRINT HTLS, HTLS$
PRINT MHT%, MHT$
PRINT CON, CON$
INPUT I$

PRINT #OUTFILE%, "name= "; name$
PRINT #OUTFILE%, ZZC, ZZC$
PRINT #OUTFILE%, ZZW, ZZW$
PRINT #OUTFILE%, XNHT, XNHT$
PRINT #OUTFILE%, PMX1, PMX1$
PRINT #OUTFILE%, SPD, SPD$
PRINT #OUTFILE%, NET%, NET$
PRINT #OUTFILE%, NOC%, NOC$
PRINT #OUTFILE%,

PRINT #OUTFILE%, SHR, SHR$
PRINT #OUTFILE%, NFI%, NFI$
PRINT #OUTFILE%, NDIV%, NDIV$
PRINT #OUTFILE%, NWR%, NWR$
PRINT #OUTFILE%, NDS%, NDS$
PRINT #OUTFILE%, DDD, DDD$
PRINT #OUTFILE%, DLL, DLL$
PRINT #OUTFILE%,
PRINT #OUTFILE%, BTC, BTC$
PRINT #OUTFILE%, BTW, BTW$
PRINT #OUTFILE%, BPD, BPD$
PRINT #OUTFILE%, BST, BST$
PRINT #OUTFILE%, BPL, BPL$
PRINT #OUTFILE%, BRC, BRC$
PRINT #OUTFILE%, BEC, BEC$
     PRINT #OUTFILE%,

'READ PDR, BRO, BDR, BWD, BRL, MBR%, MCT%
PRINT #OUTFILE%, PDR, PDR$
PRINT #OUTFILE%, BRO, BRO$
PRINT #OUTFILE%, BDR, BDR$
PRINT #OUTFILE%, BWD, BWD$
PRINT #OUTFILE%, BRL, BRL$
PRINT #OUTFILE%, MBR%, MBR$
PRINT #OUTFILE%, MCT%, MCT$
PRINT #OUTFILE%,
'READ CTD, CTLL, CTLS, HTD, HTLL, HTLS, MHT%
PRINT #OUTFILE%, CTD, CTD$
PRINT #OUTFILE%, CTLL, CTLL$
PRINT #OUTFILE%, CTLS, CTLS$
PRINT #OUTFILE%, HTD, HTD$
PRINT #OUTFILE%, HTLL, HTLL$
PRINT #OUTFILE%, HTLS, HTLS$
PRINT #OUTFILE%, MHT%, MHT$
PRINT #OUTFILE%,


'Hot and cold metal temperature are set at given hot and cold
'gas temperature:
BTW1 = BTW
BTC1 = BTC

'Power piston area
UPA = BPD * BPD * PI4
PRINT "UPA = power piston area ="; UPA; " CM^2"

'temperarure ratio, inverse of normal
UTR = BTW / BTC

'hot cylinder volume amplitude is reference volume for all
'dimensionless volumes

HAV = UPA * BST / 2
PRINT "HAV = hot active volume amplitude ="; HAV; " CM^3"

GGV = PIE * BPD * BPL * BRC
PRINT "GGV = dead volume at side of hot cap ="; GGV; " CM^3"


'initialize loop counter set for 3 iterations
'LUP% = 1

'339 **** Return point after 1 cycle ****
FOR lup = 1 TO 3

'THE FOLLOWING ARE REDUCED DIMENSIONLESS DEAD VOLUMES RELATING TO
'MAXIMUM AND MINIMUM EFFECTIVE GAS TEMP, AND HAV
HGV = GGV * 2 * BTW / (HAV * (BTW + BTC))
HEC = (UPA - PDR * PDR * PI4) * BEC * UTR / HAV
AFC = PI4 * CTD * CTD * MCT%'cooler_free_flow_area
HCV = (AFC * CTLL + DDD * DDD * PI4 * DLL * MBR%) * UTR / HAV'cooler+cold_ducting vol.
AFR = PI4 * BDR * BDR * MBR% * BRO          'regenerator flow aeria
HRV = AFR * BRL * 2 * BTW / (HAV * (BTW + BTC))' reduced regenerator dead volume
AFH = PI4 * HTD * HTD * MHT%  'heater free flow earia
HHV = AFH * HTLL / HAV       'reduced heater dead volume
HHC = UPA * BEC / HAV        'reduced hot clearence volume
VD = HGV + HEC + HCV + HRV + HHV + HHC

'52 The following are hot, regenerator and cold dead volumes in cm^2
VHD = AFH * HTLL + UPA * BEC  'HOT_DEAD_VOLUME
VRD = GGV + AFR * BRL         'REGENERATOR_DEAD_VOLUME   VRD
VCD = (UPA - PDR * PDR * PI4) * BEC + HCV * HAV / UTR   'COLD_DEAD_VOLUME

IF lup = 1 THEN
PRINT "HGV = reduced hot gap dead volume ="; HGV
PRINT "HEC = cold end clearence dead volume dimensionless ="; HEC
PRINT "AFC = cooler free flow area ="; AFC; " CM^2"
PRINT "HCV = cooler + cold ducting volume ="; HCV
PRINT "AFR = regenerator flow aeria ="; AFR; " CM^2"
PRINT "HRV = reduced regenerator dead volume ="; HRV
PRINT "AFH = heater free flow earia ="; AFH; " CM^2"
PRINT "HHV = reduced heater dead volume ="; HHV
PRINT "HHC = reduced_hot_clearance_dead_volume ="; HHC
PRINT "VD  = reduced dead volume="; VD

INPUT I$


 HAC = HAV - PI4 * PDR * PDR * BST / 2

CALL makestr(UPA * BEC, 12, 4, printv$): PRINT printv$;
PRINT " = piston end clearance dead volume, cm^3"
CALL makestr(AFH * HTLL, 12, 4, printv$): PRINT printv$;
PRINT " = heater tube dead volume, cm^3"
PRINT "+ -----------"
CALL makestr(VHD, 12, 4, printv$): PRINT printv$;
PRINT " = hot dead volume, cm^3"
CALL makestr(2 * HAV, 12, 4, printv$): PRINT printv$;
PRINT " = hot active volume, cm^3"
PRINT "+ -----------"
CALL makestr(VHD + 2 * HAV, 12, 4, printv$): PRINT printv$;
PRINT " = total max hot volume, cm^3"
PRINT
CALL makestr((UPA - PDR * PDR * PI4) * BEC, 12, 4, printv$): PRINT printv$;
PRINT " = piston end clearance dead volume, cm^3"
CALL makestr(AFC * CTLL, 12, 4, printv$): PRINT printv$;
PRINT " = cooler tube dead volume, cm^3"
CALL makestr(HCV * HAV / UTR - AFC * CTLL, 12, 4, printv$): PRINT printv$;
PRINT " = connecting ducts dead volume, cm^3"
PRINT "+ -----------"
CALL makestr(VCD, 12, 4, printv$): PRINT printv$;
PRINT " = cold dead volume, cm^3"
CALL makestr(2 * HAC, 12, 4, printv$): PRINT printv$;
PRINT " = cold active volume, cm^3"
PRINT "+ -----------"
CALL makestr(VCD + 2 * HAC, 12, 4, printv$): PRINT printv$;
PRINT " = total max cold volume, cm^3"
PRINT
CALL makestr(AFR * BRL, 12, 4, printv$): PRINT printv$;
PRINT " = regenerator dead volume, cm^3"
CALL makestr(GGV, 12, 4, printv$): PRINT printv$;
PRINT " = hot gap dead volume, cm^3"
PRINT "+ -----------"
CALL makestr(VRD, 12, 4, printv$): PRINT printv$;
PRINT " = regenerator dead volume, cm^3"
PRINT
CALL makestr(VRD + VCD + VHD, 12, 4, printv$): PRINT printv$;
PRINT " = total dead volume, cm^3",

INPUT I$
END IF

'on the first iteration only the input values are written out
'58 IF(LUP%-1)343,343,346

IF lup = 1 THEN
'343
 HAC = HAV - PI4 * PDR * PDR * BST / 2
 PRINT 'write(6,1805)
 PRINT "Cold metal temp (k) ="; BTC, ,
 PRINT "Hot metal temp  (k) ="; BTW
 PRINT "Piston diameter(cm) ="; BPD,
 PRINT "Piston stroke  (cm) ="; BST

 PRINT "Hot cap length (cm)  ="; BPL, ,
 PRINT "Piston gap     (cm)  ="; BRC
 PRINT "Piston end clearence(cm)="; BEC,
 PRINT "Piston rod diameter (cm)="; PDR

 PRINT "Regenerator porosity    ="; BRO,
 PRINT "Regenerator diameter(cm)="; BDR
 PRINT "Regenerator wire dia(cm)="; BWD,
 PRINT "Regenerator length  (cm)="; BRL

 PRINT "Number of regenerators   ="; MBR%,
 PRINT "Number of cooler tubes   ="; MCT%
 PRINT "Cool tube diameter  (cm) ="; CTD,
 PRINT "Total cool tube length(cm)="; CTLL
' INPUT I$
 PRINT "Cool ct length       (cm)="; CTLS,
 PRINT "Heater tube diameter (cm)="; HTD
 PRINT "Total heater tube length(cm)="; HTLL,
 PRINT "Heat heater tube length (cm)="; HTLS
 PRINT "Number of heater tubes      ="; MHT%,
 PRINT "Heat transfer exponent      ="; XNHT

 PRINT "Average pressure MPA ="; PMX1,
 PRINT "Engine speed (rad/sec)="; SPD
 PRINT "NETNET option         ="; NET%, ,
 PRINT "Number of cylinders   ="; NOC%

 PRINT "Cold crank ratio ="; ZZC, ,
 PRINT "Hot crank ratio  ="; ZZW
 GOSUB printinput
 'INPUT I$
END IF


'68 The reduced dead volume is divided into fractions, re-evaluated
'each iteration because of change in tempertatures
'346
H(1) = (HEC + HCV / 2) / VD
H(2) = H(1) + HCV / (2 * VD)
H(3) = H(2) + HRV / (2 * VD)
H(4) = H(2) + HRV / VD
H(5) = H(4) + HHV / (2 * VD)
'75 Calculations from 349 to 295 are on first iteration only. Engine
'volumes and volume derivatives are calculated and decision matrix
'is defined.
'IF (LUP% - 1) 349,349,295

IF lup = 1 THEN
 '349
 XND = NDIV%
 NDIV1% = NDIV% + 1
 'DALF=0.0174 radians/incriment
 DALF = 2 * PIE / XND'  6.2831853# / XND
 NT% = NDIV% / 4 + 2
 NE% = NT% - 1  'NOT USED ?????????????????
 '85 Calculation starts at 270 deg=4.71 radians in radians if the
 'front side of the cold piston is used and at 90 deg=1.57 radians
 'in radians if the back side is used and room must be allowed for
 'a piston drive rod. This gives the proper curve shape.
 'Calculation always starts with zero cold live volume.
 'IF (PDR) 4060,4060,4070
 IF PDR > 0 THEN
   ALF = PIE / 2' 1.5707963#
 ELSE
   ALF = PIE + PI4' 4.7123889#
 END IF
 NF% = NDIV% / 4
 '95 Call subroutine to calculate duplacable space above or under cold
 'piston at the midpoint and at the beginning of each angle
 'increment as a fraction of the piston stroke amplitude.
 'Subroutine also calculates derivatives to be used later.
 CALL VOLC(DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(), COFI(), SALF(), CALF(), PDR, ALF)
 'subroutine VOLC(DALF,NF%,C(),CI(),DC(),DCI(),ZZC,NDIV%,SIFI(),COFI(),SALF(),CALF(),PDR,ALF)
 '101 The phase angle must be 90 deg. for a heat engine and 270 deg for
 'a heat pump
  FI = DALF * NFI%
 FI1 = FI * 180 / PIE
 'WRITE(6,11) SHR, NDIV%, FI1, NDS%, DDD, DLL
 'CALL makestr(SHR, printv$)
 PRINT "Specific heat ratio ="; SHR, ', printv$,
 'CALL makestr(NDIV%, printv$)
 PRINT , "Div. per cycle ="; NDIV%' , printv$
' CALL makestr(FI1, printv$)
 PRINT "Phase angle (deg.) ="; FI1, , ', printv$,
' CALL makestr(NDS%, printv$)
 PRINT "Incr. in dead space int.="; NDS%', printv$
' CALL makestr(DDD, printv$)
 PRINT "Duct diameter (cm) ="; DDD, ', printv$,
 'CALL makestr(DLL, printv$)
 PRINT , "Duct length (cm)="; DLL', printv$
 INPUT I$
 SFI = SIN(FI)
 CFI = COS(FI)
 'Now the hot space fractions axd, derivatives are calculated for
 'each increment
 CALL VOLW(W(), WI(), DW(), DWI(), CFI, SFI, ZZW, NDIV%, SIFI(), COFI(), SALF(), CALF(), DALF)
 'choice matrix is defined- see note #2
 IND(1, 1) = 1
 IND(1, 2) = 3
 IND(2, 1) = 4
 IND(2, 2) = 2
END IF
'First guess at initial quantities
'Counter NN% is 1 at start OF cycle when Give cold volume is dsad
'It becomes 2 after live hot volume becomes zero, phase angle is
'270 deg. for a heat pump and 90 deg. for a heat engine.
'295
NN% = 1
'None of the gas mass is assumed to be in the cold space to start,
'volume is zero
XMC = 0
'Relative pressure is assumed to be maximum, approximation.
P = 1
'At start all gas mass is assumed to be in the hot space, ignores
'dead volume
XMW = 1 - CFI
'mass from pevious cycle is assumed at start to be zero, to assure
'at least 2 cycles to convergence.
XMWS = 0
'previous cycle pressure
PS = 1
'135 initial assumption for value of IND(k,l) is correct for heat
'engine wrong for a heat pump.
NO = 1
'initialize dimensionless works
WW = 0
WC = 0
NITE% = 1
NST% = 1
NFIN% = NFI%
'displaced mass ratio

RVT = HAC * UTR / HAV
CI(NDIV1%) = CI(1)
WI(NDIV1%) = WI(1)
'148 ******************************************
'start of main do loop, return point after each increment

RELOOP = 1
WHILE RELOOP = 1
RELOOP = 0
'404 'do 102 I=NST%,NFIN%

FOR I = NST% TO NFIN%

'transfer volumes and derivtives from storage
VW = W(I)
VC = C(I)
VWI = WI(I)
VCI = CI(I)
DVW = DW(I)
DVC = DC(I)
DVWI = DWI(I)
DVCI = DCI(I)
'splits to 4 options
'go to (201,202,203,204),NO
'IF NO = 1 THEN gOTO 201
'IF NO = 2 THEN gOTO 202
'IF NO = 3 THEN gOTO 203
'IF NO = 4 THEN gOTO 204

IF NO = 1 THEN
'201 integration program for mass increasing in both hot and cold spaces.
'NO=1 (see note #3
'Compute pressure change based upon initial conditions
 DP = -SHR * P * (RVT * DVCI + DVWI) / (RVT * VCI + VWI + SHR * VD) * DALF

'change in pressure                     dispMr coldV  hotV heatR rdeadV
'finds pressure at mid increment
S = P + DP / 2

'calculates final pressure change basic upon mid point values
DP = -SHR * S * (RVT * DVC + DVW) / (RVT * VC + VW + SHR * VD) * DALF
'calculates mass changes
DMW = S * DVW * DALF + VW * DP / SHR
DMC = -(DMW + VD * DP) / RVT
'determins choice matrix

IF DMW >= 0 THEN
  K% = 1
ELSE
  K% = 2
END IF
IF DMC >= 0 THEN
  L% = 1
ELSE
  L% = 2
END IF
NO = IND(K%, L%)
'if choice is changed next iteration will be through a different
'option
'gOTO 400


ELSEIF NO = 2 THEN

'202 integration program for mass decreasing in both hot and cold spaces.
'NO=2 (see option 1 for detailed explanation)

IF XMC < 0 THEN
  XMC = 0!
END IF
IF XMW < 0 THEN
  XMW = 0
END IF
DP = -SHR * (XMC * RVT * DVCI / VCI + XMW * DVWI / VWI) / (XMC * RVT / P + XMW / P + SHR * VD) * DALF
DMC = XMC * (DVCI * DALF / VCI + DP / SHR / P)
DMW = -RVT * DMC - VD * DP
S = P + DP / 2
SMC = XMC + DMC / 2
SMW = XMW + DMW / 2
DP = -SHR * (SMC * RVT * DVC / VC + SMW * DVW / VW) / (SMC * RVT / S + SMW / S + SHR * VD) * DALF
DMC = SMC * (DVC * DALF / VC + DP / SHR / S)
DMW = -RVT * DMC - VD * DP

IF DMW <= 0 THEN
  K% = 2
ELSE
  K% = 1
END IF
IF DMC <= 0 THEN
  L% = 2
ELSE
  L% = 1
END IF
NO = IND(K%, L%)
'gOTO 400


ELSEIF NO = 3 THEN
'203 integration program for mass decreasing in cold space and increasing
'in hot space, NO=3 (see option 1 for detailed explanation)
' if(XMC) 704,703,703
IF XMC < 0 THEN
  XMC = 0
END IF
DP = -SHR * (P * DVWI + XMC * RVT * DVCI / VCI) / (VWI + XMC * RVT / P + SHR * VD) * DALF
DMC = XMC * (DVCI * DALF / VCI + DP / SHR / P)
DMW = -RVT * DMC - VD * DP
S = P + DP / 2
SMC = XMC + DMC / 2
SMW = XMW + DMW / 2
DP = -SHR * (S * DVW + SMC * RVT * DVC / VC) / (VW + SMC * RVT / S + SHR * VD) * DALF
DMC = SMC * (DVC * DALF / VC + DP / SHR / S)
DMW = -RVT * DMC - VD * DP
IF DMW >= 0 THEN
  K% = 1
ELSE
  K% = 2
END IF

IF DMC <= 0 THEN
  L% = 2
ELSE
  L% = 1
END IF
NO = IND(K%, L%)
'gOTO 400


ELSEIF NO = 4 THEN
'204 Integration program for mass decreasing in cold space and decreasing
'in hot space, NO=4 (see option 1 for detailed explanation)
IF XMW < 0 THEN
  XMW = 0
END IF
DP = -SHR * (P * RVT * DVCI + XMW * DVWI / VWI) / (RVT * VCI + XMW / P + SHR * VD) * DALF

DMW = XMW * (DVWI * DALF / VWI + DP / SHR / P)
DMC = -(DMW + VD * DP) / RVT
S = P + DP / 2
SMC = XMC + DMC / 2
SMW = XMW + DMW / 2
DP = -SHR * (S * RVT * DVC + SMW * DVW / VW) / (RVT * VC + SMW / S + SHR * VD) * DALF
DMW = SMW * (DVW * DALF / VW + DP / SHR / S)
DMC = -(DMW + VD * DP) / RVT
IF DMW <= 0 THEN
  K% = 2
ELSE
  K% = 1
END IF
IF DMC >= 0 THEN
  L% = 1
ELSE
  L% = 2
END IF
NO = IND(K%, L%)
'gOTO 400
ELSE
  INPUT "NO IS NOT 1,2,3 OR 4"; I$
END IF

'400 increment pressure and mass
P = P + DP
XMC = XMC + DMC
XMW = XMW + DMW
'calculates works
PW = P - DP / 2
WC = WC + PW * DVC * DALF
WW = WW + PW * DVW * DALF
'records results into arrays
PR(I) = P
DPR(I) = DP
XMCX(I) = XMC
XMWX(I) = XMW
XDMC(I) = DMC
XDMW(I) = DMW
' ******** end of main loop **********

NEXT I

IF NN% = 1 THEN
  'reset main loop for last part of cycle
  NST% = NFI% + 1
  NFIN% = NDIV%
  NN% = 2
  RELOOP = 1
'  gOTO 404
ELSE
'205 Tests for convergence at end of cycle. The change in the fraction
'of mass in the hot space from one cycle to the next must be less
'than 0.1%, and the change in the pressure from one cycle to the next
'must be less than 0.5%. However, NO more than 15 cycles are
'allowed
TEST = SQR((XMWS - XMW) ^ 2)
TEST1 = SQR((PS - P) ^ 2)
IF (NITE% <= 15) THEN
  IF (TEST > .001) OR (TEST1 > .005) THEN

'Reinitialise for next cycle
 NN% = 1
 XMC = 0
 PS = P
 XMWS = XMW
 WW = 0
 WC = 0
 NST% = 1
 NFIN% = NFI%
 NITE% = NITE% + 1
 NO = 4
 RELOOP = 1
' gOTO 404

  END IF
END IF
END IF
WEND'RELOOP


'307 The dimendionless pressures and works have been calcolated for one
'cycle. Now the additional heat and power losses will be calculated.
'Calculate average dimensionless pressure.
PAVG = 0
'311
FOR I = 1 TO NDIV%
  PAVG = PAVG + PR(I)
NEXT I
PAVG = PAVG / NDIV%
'determine maximum and minimum dimensionless pressure
PMAX = xlarge(PR(), NDIV%)
PMIN = small(PR(), NDIV%)
'Adjust dimensionless works to relate to newly determend maximum
'pressure
WC = WC / PMAX
WW = WW / PMAX
'pressure ratio

RP = PMAX / PMIN
'Find maximum masses and adjust them to maximum pessure.
CMMAX = xlarge(XMCX(), NDIV%)
WMMAX = xlarge(XMWX(), NDIV%)
CMMAX = CMMAX / PMAX
WMMAX = WMMAX / PMAX
'calc. max. pressure, mpa
PMX = PMAX * PMX1 / PAVG
'Calculates angle between pressure wave and volume wave for a heat
'engine.
ARG = 2 * RP / (RP - 1) * WW / PIE' 3.1416
'if (1-ARG**2) 1607,1608,1608

IF 1 - ARG ^ 2 < 0 THEN RETURN'gOTO 1607
'  FIPV = arsin(ARG)
'   FIPU = ATN(SIN(ARG) / SQR(1 - SIN(ARG) ^ 2))
  FIPV = ATN(ARG / (SQR(1 - ARG ^ 2)))
  XNDS = NDS%
  'Calculates values used in flow loss calculations and flow integrals
X = 0
DX = 1 / XNDS
NIN% = NDS% + 1

COR = PMAX ^ (XNHT - 2) * DALF ^ (XNHT - 1)
'do 854 I=1,NIN%
FOR I = 1 TO NIN%
'GOSUB PDINT
  CALL pdint(X, XDMW(), XDMC(), RVT, DC(), NDIV%, DMRE, PR(), XINT, DPR(), XI1, XI2, XNHT)
  XINT = XINT / DALF / PMAX
  DMRE = DMRE / PMAX / (2 * PIE)' 6.2832
  XI1 = XI1 * COR / (1.5708 * DMRE) ^ (1 - XNHT)
  XI2 = XI2 * COR / (1.5708 * DMRE) ^ (2 - XNHT)
  XI3 = XI1 / XI2
  GDMS(I) = DMRE
  GINT(I) = XINT
  GI2(I) = XI2
  GI3(I) = XI3
  X = X + DX
NEXT I
'854 continue
'Interpolates flow integrals


FOR I = 1 TO 5
  UIN(I) = plot(GINT(), H(I))
  UDM(I) = plot(GDMS(), H(I))
NEXT I

UI23 = plot(GI2(), H(2))
UI24 = plot(GI2(), H(4))
UI33 = plot(GI3(), H(2))
UI34 = plot(GI3(), H(4))
' **** Calculation of constants *****
'specific for hydrogen gas
HMU = 8.873E-05 + .0000002 * (BTW - 293)
CMU = 8.873E-05 + .0000002 * (BTC - 293)
BTR = (BTW - BTC) / LOG(BTW / BTC)
RMU = 8.873E-05 + .0000002 * (BTR - 293)
CP1 = 14.6
CV1 = 10.46
R2 = 8.23168E+07
R = 4.116
'373 ******** Cold exchanger pressure drop **********
REC = UDM(1) * PMX * SPD * HAC * CTD / (BTC * AFC * CMU * R)

IF (REC <= 2000) THEN
  FC = 16 / REC
ELSE
  FC = EXP(-1.34 - .2 * LOG(REC))
END IF
GLS = CTLL * SPD * SPD * HAC * HAC * FC * UIN(1) / (CTD * AFC * AFC * BTC * R2)
QP = NOC% * SPD * PMX * HAC / (2 * PIE)
QCP = QP * GLS

'********** Hot exchanger pressure drop ***********
REH = UDM(5) * PMX * SPD * HAC * HTD / (BTW * AFH * HMU * R)
'384 if(REH-2000) 1988,1988,1989
IF (REH <= 2000) THEN
  FH = 16 / REH
ELSE
  FH = EXP(-1.34 - .2 * LOG(REH))
END IF
GLH = HTLL * SPD ^ 2 * HAC ^ 2 * BTW * FH * UIN(5) / (HTD * AFH ^ 2 * BTC ^ 2 * R2)
QHP = QP * GLH
'**** Screen---netnet option *********
RER = PMX * HAC * SPD * BWD / (AFR * R)
RE(1) = RER * UDM(2) / (BTC * CMU)
RE(2) = RER * UDM(3) / (BTR * RMU)
RE(3) = RER * UDM(4) / (BTW * HMU)

FOR I = 1 TO 3
 IF (NET%) <= 0 THEN
   IF (RE(I) <= 60) THEN
     FR(I) = EXP(1.73 - .93 * LOG(RE(I)))
   ELSE
     IF (RE(I) <= 1000) THEN
       FR(I) = EXP(.714 - .365 * LOG(RE(I)))
     ELSE
       FR(I) = EXP(.015 - .125 * LOG(RE(I)))
     END IF
   END IF
 ELSE
   FR(I) = 2.73 * (1 + 10.397 / RE(I))
 END IF
NEXT I
'continue
'407 ******* regenerator pressure drop *****
GLR = BRL * SPD * SPD * HAC * HAC / (BWD * AFR * AFR * R2 * BTC)
QR1 = QP * GLR * UIN(2) * FR(1)
QR2 = QP * GLR * UIN(3) * FR(2) * BTR / BTC
QR3 = QP * GLR * UIN(4) * FR(3) * UTR
QRP = (QR1 + QR3 + 4 * QR2) / 6
'Calculates effective hot and cold gas temperatures based upon the
'number of transfer units in the heat exchangers, specific for
'hydrogen
CNTU = .112 * CTLS / (CTD * REC ^ .2)
DTC = WC * (SHR - 1) / (2 * UDM(1) * SHR * (EXP(2 * CNTU) - 1))
BTC = BTC1 * (1 - DTC)
HNTU = .1044 * HTLS / (HTD * REH ^ .2)
DTH = WW * (SHR - 1) / (2 * UDM(5) * SHR * (EXP(2 * HNTU) - 1))
BTW = BTW1 * (1 - DTH)
'Note, temperature ratio is redefined for next iteration
UTR = BTW / BTC
'424 *** Reheat loss ****
RNTU = BRL * 4.37 / (BWD * SQR(PI4 * 2 * RE(1)))
QNTU = BRL * 4.031 / (BWD * SQR(PI4 * 2 * RE(3)))
QNPH = AFR * BRL * .195 / (PI4 * HAC * UDM(2) * (UTR - 1))
QDK = QNPH * (UI33 + UI34 * UDM(2) / UDM(4)) / 2
QLM = (1 + QDK) / (RNTU / UI23 + QNTU * UDM(2) / (UDM(4) * UI24))
QHR = UDM(2) * CP1 * (BTW - BTC) * SPD * PMX * HAC * QLM * NOC% / (R * BTC * 2)
'**** shuttle loss *****
QL1 = 231.2 * SQR(SPD * BRC * BRC)
QB = (2 * QL1 * QL1 - QL1) / (2 * QL1 * QL1 - 1)
QHC = .00146 * BST * (BTW - BTC) * PI4 * BPD * BST * QB * NOC% / (BRC * BPL)
'**** Pumping loss *******
IF (MS% AND 1) = 0 THEN
  QFS = (RP / (BTW / (BTW - 2 * BTC) - BST / BPL)) + (1 / (BTW / ((BTW - 2 * BTC) + BST / BPL)))
ELSE
' Martini writes in the manual the calculation above is wrong !!!!!!!!!!!!!!
' has to be
  X = (BTW + BTC1) / (BTW - BTC1)
  Y = BST / BPL
  QFS = -RP / (X + Y) + 1 / (X - Y)
END IF
QHG = ABS(SPD * PMX * GGV * BST * SHR * QFS * ARG * NOC% / ((SHR - 1) * BPL * RP * 8))
'**** Basic powers *****
'PO = (WW * HAV + WC * HAC) * (4.5) * PMX * SPD * NOC% / PIE
PO = (WW * HAV + WC * HAC) * (.5) * PMX * SPD * NOC% / PIE
'441 ****** Net power ********
POF = PO - QCP - QHP - QRP
'Get ready to report on one iteration and prepare for the next.
'Reset hot end dimensionless heat transfer integral
HTW = 0
'The program tries to keep PMAX=1. This adjustment of the pressure
'and mass does this
'448 do 509 I=1,NDIV%
FOR I = 1 TO NDIV%
  PR(I) = PR(I) / PMAX
  XMCX(I) = XMCX(I) / PMAX
  XMWX(I) = XMWX(I) / PMAX
NEXT I
'452 Dimensionless hot and cold gas temperatures for esch increment.
'if they are less than zero correct to zero.
WI(NDIV1%) = WI(1)
CI(NDIV1%) = CI(1)

FOR I = 1 TO NDIV%
  IF (XMCX(I)) > 0 THEN
    TEC(I) = PR(I) * CI(I + 1) / XMCX(I)
  ELSE
    TEC(I) = 0
  END IF
  IF (XMWX(I)) > 0 THEN
    TEW(I) = PR(I) * WI(I + 1) / XMWX(I)
  ELSE
    TEW(I) = 0
  END IF
NEXT I
'466 Dimensionless average hot and cold gas temperatures for ful cycle.
TEW(NDIV1%) = TEW(1)
TEC(NDIV1%) = TEC(1)
PR(NDIV1%) = PR(1)
XMCX(NDIV1%) = XMCX(1)
XMWX(NDIV1%) = XMWX(1)
TWDM = 0
TCDM = 0
'474 do 573 I=1,NDIV%
FOR I = 1 TO NDIV%
  DMW = XMWX(I + 1) - XMWX(I)
  IF (DMW) < 0 THEN
    TMPW = (TEW(I) + TEW(I + 1)) / 2
    TWDM = TWDM + (TMPW - 1) * DMW
  END IF
  DMC = XMCX(I + 1) - XMCX(I)

  IF (DMC) < 0 THEN
    TMPC = (TEC(I) + TEC(I + 1)) / 2
    TCDM = TCDM + (TMPC - 1) * DMC
  END IF
 
NEXT I
TWDM = TWDM * SHR / (SHR - 1)
TCDM = TCDM * SHR / (SHR - 1)
'Hot end heat transfer integral for full cycle and total gas mass at
'each point in the cycle. Total mass should not change.
'408 do 1021 I=1,NDIV%
FOR I = 1 TO NDIV%
  HTW = HTW + (WI(I + 1) - WI(I)) * (PR(I) + PR(I + 1)) / 2
1021 XMT(I) = XMCX(I) * RVT + XMWX(I) + PR(I) * VD
NEXT I
'Basic heat input, watts
HT = HTW * SPD * PMX * HAV * NOC% / (2 * PIE)
'Specific static condution heat loss for the 4l23 engine
CON = 9680  'put as input
'flow friction credit, watts
FFF = (QHP + .5 * QRP) * (-1)
'Heat to engine, watts
HTE = HT + QHR + QHC + QHG + CON + FFF
'indicated efficiency, %
ZEF = 100 * POF / HTE
'print out results of one iteration
'write(6,12) LUP%
'write(6,3010)PO,HT
'write(6,3020) QHP,QHR
'write(6,1925) QRP,QHC,QCP,QHG,POF,CON,ZEF,FFF,HTE
'write(6,1921) BTW,BTC,RVT,VD
PRINT "Iteration "; lup
PRINT "Basic power    (watts)="; PO,
PRINT "Basic heat     (watts)="; HT
'PRINT
PRINT "Heater windage (watts)="; QHP,
PRINT "Reheat loss    (watts)="; QHR

PRINT "Regen. windage (watts)="; QRP,
PRINT "Shuttle loss (watts)  ="; QHC
PRINT "Cooler windage (watts)="; QCP,
PRINT "Appendix loss  (watts)="; QHG
PRINT "Net power (watts)     ="; POF,
PRINT "Conduction (watts)    ="; CON
PRINT "Indicated efficiency(%)"; ZEF,
PRINT "Flow friction (watts) ="; FFF
PRINT , , , "Heat to engine (watts)="; HTE

PRINT "Effective hot temperature (k)="; BTW,
PRINT "Effective cold temperature(k)="; INT(BTC * 10) / 10
PRINT "Displaced mass ratio  ="; RVT,
PRINT "Reduced dead volume   ="; VD
INPUT I$
GOSUB printresult
'After all losses are taken into account LUP% is indexed. The program
'does 3 iterations with printouts before going into a summary
'LUP% = LUP% + 1
'510 if(LUP%-3) 339,339,1607
NEXT lup

'if input volume NWR% is other than zero the following summary
'Information is printed at the end of the computation
RETURN

printout:
'1607
'IF (NWR%) = 0 THEN gOTO 242
'1613 write(6,51)TWDM,TCDM
PRINT "Dimensionless average bas temp. hot end ="; TWDM
PRINT "Dimensionless average bas temp. cold end="; TCDM
INPUT I$
'Print out each 10 degrees, angle, hot volume, hot gas temperature.
'cold volume, cold gas temperature, total volume, pressure
'1149  write(6,20)

'518 do 3001 I=10,NDIV%,10
xmaximum = 0
vtmaximum = 0
CLS
PRINT
PRINT "   Pressure,MPA-volume,CN3 data for one cylinder"
PRINT " =============================================================================="
PRINT "angle de.    cold vol.cm^3  pressure MPa ";
PRINT "angle de.  cold vol.cm^3   pressure MPa"
PRINT "     hot vol.cm^3    total vol.cm^3";
PRINT "          hot vol.cm^3    total vol.cm^3"

FOR I = 10 TO NDIV% STEP 10
 IF I <= NDIV% / 2 THEN
   hpos = 1
   vpos = 5 + I / 10
 ELSE
   hpos = 41
   vpos = 5 + (I - NDIV% / 2) / 10
 END IF
 X = PR(I) * PMX
 VH = VHD + HAV * WI(I)
 VC = VCD + CI(I) * HAC
 VT = VH + VRD + VC
 IF VT > vtmaximum THEN vtmaximum = VT
 IF X > xmaximum THEN xmaximum = X
 
 LOCATE vpos, hpos
 IF I < 100 THEN PRINT " ";
 PRINT I;
 CALL makestr(VH, 8, 2, printv$): PRINT printv$;
 CALL makestr(VC, 8, 2, printv$): PRINT printv$;
 CALL makestr(VT, 8, 2, printv$): PRINT printv$;
 CALL makestr(X, 8, 2, printv$): PRINT printv$;

' IF I = 140 THEN INPUT I$
NEXT I

PRINT #OUTFILE%,
PRINT #OUTFILE%, "Pressure,MPA-volume,CN3 data for one cylinder"
PRINT #OUTFILE%, "angle de. hot vol.cm^3",
PRINT #OUTFILE%, "cold vol.cm^3",
PRINT #OUTFILE%, "total vol.",
PRINT #OUTFILE%, "pressure MPa"

FOR I = 10 TO NDIV% STEP 10
 X = PR(I) * PMX
 VH = VHD + HAV * WI(I)
 VC = VCD + CI(I) * HAC
 VT = VH + VRD + VC
 IF I < 100 THEN PRINT #OUTFILE%, " ";
 PRINT #OUTFILE%, I;
 CALL makestr(VH, 13, 7, printv$): PRINT #OUTFILE%, printv$,
 CALL makestr(VC, 13, 7, printv$): PRINT #OUTFILE%, printv$,
 CALL makestr(VT, 13, 7, printv$): PRINT #OUTFILE%, printv$,
 CALL makestr(X, 13, 9, printv$): PRINT #OUTFILE%, printv$
NEXT I

PRINT
'524 Starts over with next DATA SET
RETURN


printinput:
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Cold metal temp (k) ="; BTC
PRINT #OUTFILE%, "Hot metal temp  (k) ="; BTW
PRINT #OUTFILE%, "Piston diameter(cm) ="; BPD
PRINT #OUTFILE%, "Piston stroke  (cm) ="; BST
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Hot cap length (cm)  ="; BPL
PRINT #OUTFILE%, "Piston gap     (cm)  ="; BRC
PRINT #OUTFILE%, "Piston end clearence(cm)="; BEC
PRINT #OUTFILE%, "Piston rod diameter (cm)="; PDR
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Regenerator porosity    ="; BRO
PRINT #OUTFILE%, "Regenerator diameter(cm)="; BDR
PRINT #OUTFILE%, "Regenerator wire dia(cm)="; BWD
PRINT #OUTFILE%, "Regenerator length  (cm)="; BRL
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Number of regenerators   ="; MBR%
PRINT #OUTFILE%, "Number of cooler tubes   ="; MCT%
PRINT #OUTFILE%, "Cool tube diameter  (cm) ="; CTD
PRINT #OUTFILE%, "Total cool tube length(cm)="; CTLL
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Cool ct length       (cm)="; CTLS
PRINT #OUTFILE%, "Heater tube diameter (cm)="; HTD
PRINT #OUTFILE%, "Total heater tube length(cm)="; HTLL
PRINT #OUTFILE%, "Heat heater tube length (cm)="; HTLS
PRINT #OUTFILE%, "Number of heater tubes      ="; MHT%
PRINT #OUTFILE%, "Heat transfer exponent      ="; XNHT
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Average pressure MPA ="; PMX1
PRINT #OUTFILE%, "Engine speed (rad/sec)="; SPD
PRINT #OUTFILE%, "NETNET option         ="; NET%
PRINT #OUTFILE%, "Number of cylinders   ="; NOC%
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Cold crank ratio ="; ZZC
PRINT #OUTFILE%, "Hot crank ratio  ="; ZZW
RETURN

printresult:
PRINT #OUTFILE%,
PRINT #OUTFILE%, "Iteration "; lup
PRINT #OUTFILE%, "Basic power    (watts)="; PO,
PRINT #OUTFILE%, "Basic heat     (watts)="; HT
'PRINT
PRINT #OUTFILE%, "Heater windage (watts)="; QHP,
PRINT #OUTFILE%, "Reheat loss    (watts)="; QHR

PRINT #OUTFILE%, "Regen. windage (watts)="; QRP,
PRINT #OUTFILE%, "Shuttle loss (watts)  ="; QHC
PRINT #OUTFILE%, "Cooler windage (watts)="; QCP,
PRINT #OUTFILE%, "Appendix loss  (watts)="; QHG
PRINT #OUTFILE%, "Net power (watts)     ="; POF,
PRINT #OUTFILE%, "Conduction (watts)    ="; CON
PRINT #OUTFILE%, "Indicated efficiency(%)"; ZEF,
PRINT #OUTFILE%, "Flow friction (watts) ="; FFF
PRINT #OUTFILE%, , , , "Heat to engine (watts)="; HTE

PRINT #OUTFILE%, "Effective hot temperature (k)="; BTW,
PRINT #OUTFILE%, "Effective cold temperature(k)="; INT(BTC * 10) / 10
PRINT #OUTFILE%, "Displaced mass ratio  ="; RVT,
PRINT #OUTFILE%, "Reduced dead volume   ="; VD
RETURN

9900 'read variable text strings
RESTORE 9902'variable text
READ ZZC$, ZZW$, XNHT$, PMX1$, SPD$, NET$, NOC$
READ SHR$, NFI$, NDIV$, NWR$, NDS$, DDD$, DLL$
READ BTC$, BTW$, BPD$, BST$, BPL$, BRC$, BEC$
READ PDR$, BRO$, BDR$, BWD$, BRL$, MBR$, MCT$
READ CTD$, CTLL$, CTLS$, HTD$, HTLL$, HTLS$, MHT$
READ CON$
RETURN

9902 'variable text
 DATA "ZZC = (input) Connecting rod length/.5 stroke for cold piston"
 DATA "ZZW = (input) Connecting rod length/.5 stroke for hot piston"
 DATA "XNHT = (input) Value for exponent in heat transfer relation of regenerator matrix"
 DATA "PMX1 = (input) Avg. pressure MPa"
 DATA "SPD = (input) Engine speed, rad/sec"
 DATA "NET% = (input) Regenerator filler option (+5 = metnet, -5 = screen)"
 DATA "NOC% = (input) number of cylinders"

DATA "SHR = (input) Specific heat ratio for working gas"
DATA "NFI% = (input) (phase angle)(NDIV%)/360"
DATA "NDIV% = (input) Number of divisions per crank rotation (must be a multiple of 4)"
DATA "NWR% = (input) Governs printout, zero for overall results only, different from zero added pv data"
DATA "NDS% = (input) Number of divisions in dead space"
DATA "DDD = (input) Cooler duct diameter, cm"
DATA "DLL = (input) Cooler duct length, cm"

DATA "BTC = (input) Effective cold temperature, K"
DATA "BTW = (input) Hot effective gas temperature, K"
DATA "BPD = (input) Piston diameter, cm"
DATA "BST = (input) Piston stroke, cm"
DATA "BPL = (input) Hot cap length, cm"
DATA "BRC = (input) Piston gap, cm"
DATA "BEC = (input) Piston end clearence, cm"

DATA "PDR = (input) Piston rod diameter, cm"
DATA "BRO = (input) Regenerator density factor"
DATA "BDR = (input) Regenerator diameter, cm"
DATA "BWD = (input) Effective regenerator wire diameter, cm"
DATA "BRL = (input) Regenerator length, cm"
DATA "MBR% = (input) Number of regenerators"
DATA "MCT% = (input) Number of cooler tubes per cylinder"

DATA "CTD = (input) Cooler tube inside diameter, cm"
DATA "CTLL = (input) Total cooler tube length, cm"
DATA "CTLS = (input) Cooler cool tube length, cm"
DATA "HTD = (input) Heater tube inside diameter, cm"
DATA "HTLL = (input) Total heater tube length, cm"
DATA "HTLS = (input) Heated heater tube length, cm"
DATA "MHT% = (input) Number of heater tubes for cylinder"
DATA "CON = Conduction loss, watts"


9999 'input
DATA "GM 4L23"
 DATA 5.874, 5.874, .204 , 17.9263, 157.080,     5    , 4
'DATA 5.874, 5.874, .204 ,  9.6516, 209.44 ,     5    , 4
DATA  1.39 , 90   , 360  , 1     , 10    , .76  , 71
 DATA 330  , 1033.15, 10.16, 4.65  , 6.4   , .0406, .0406
'DATA 330  , 1033   , 10.16, 4.65  , 6.4   , .0406, .0406
DATA  4.060, .8   , 3.5  , .0043 , 2.5   , 6    , 312
DATA  .115 , 12.9 , 12.02, .472  , 41.8  , 25.58, 36, 0.0
'CON=specific static conduction heat loss 9680 for the 4L23 engine
DATA 9680
DATA end
'integer data must be right justified. There are 7 data fields per
'line of 7 columns each
'data field layout is as follows

'ZZC, ZZW, XNHT, PMX1, SPD, NET%, NOC%
'SHR, NFI%, NDIV%, NWR%,  NDS%, DDD, DLL
'BTC, BTW, BPD,  BST,  BPL, BRC, BEC
'PDR, BRO, BDR,  BWD,  BRL, MBR%, MCT%
'CTD, CTLL, CTLS, HTD, HTLL,HTLS,MHT%

'Next data set repeats last
'change according to data printout.
'Data must be within control cards(//),

'read(5.5) ZZC,ZZW,XNHT,PMX1,SPD,NET%,NOC%
'read (5+1)  SHR,NFI%,NDIV%,NWR%,NDS%,DDD,DLL
'read (5+1710) BTC,BTW,BPD,BST,BPL,BRC,BEC
'read (5+1720) PDR,BRO,BDR,BWD,BRL,MBR%,MCT%
'read (5+1730) CTD,CTLL,CTLS,HTD,HTLL,HTLS,MHT%

2000 'input variable names
DATA BTC, BTW,BPD,BST,BPL, BRC,BEC,PDR
DATA BRO,BDR,BWD,BRL,MBR%,MCT%,CTD,CTLL
DATA CTLS,HTD,HTLL,HTLS,MHT%,XNHT
DATA PMX1,SPD,NET%, NOC%,ZZC,ZZW
DATA SHR, NDIV%, FI1, NDS%, DDD, DLL
DATA CON

2002 'inpus abreviations
DATA "Cold metal temp (k) =","Hot metal temp  (k) ="
DATA "Piston diameter(cm) =","Piston stroke  (cm) ="
DATA "Hot cap length (cm)  =","Piston gap     (cm)  ="
DATA "Piston end clearence(cm)=","Piston rod diameter (cm)="
DATA "Regenerator porosity    =","Regenerator diameter(cm)="
DATA "Regenerator wire dia(cm)=","Regenerator length  (cm)="
DATA "Number of regenerators   =","Number of cooler tubes   ="
DATA "Cool tube diameter  (cm) =","Total cool tube length(cm)="
DATA "Cool ct length       (cm)=","Heater tube diameter (cm)="
DATA "Total heater tube length(cm)=","Heat heater tube length (cm)="
DATA "Number of heater tubes      =","Heat transfer exponent      ="
DATA "Average pressure MPA =","Engine speed (rad/sec)="
DATA "NETNET option         =","Number of cylinders   ="
DATA "Cold crank ratio =","Hot crank ratio  ="
DATA "Specific heat ratio =","Div. per cycle =","Phase angle (deg.) ="
DATA "Incr. in dead space int.=","Duct diameter (cm) =","Duct length (cm)="
DATA "SPECIFIED CONDUCTION="

2006 'init load inputs abreviationsni$()
RESTORE 2000
FOR var = 1 TO nivar
 READ ni$(var)
NEXT var
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 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, length, decimals, printv$)
'call makestr(printv,length,decimals,printv$)

IF decimals = 0 THEN
  printv$ = RIGHT$("                     " + STR$(INT(printv)), length)
ELSEIF INSTR(STR$(printv), ".") THEN
  IF LEN(STR$(printv)) - INSTR(STR$(printv), ".") > decimals THEN
    printv = INT(printv * 10 ^ decimals) / (10 ^ decimals)
  END IF
  printv$ = STR$(printv) + "               "
  printv$ = LEFT$(printv$, INSTR(printv$, ".") + decimals)
'  IF LEN(printv$) - INSTR(printv$, ".") > decimals THEN
 '   printv$ = printv$ + RIGHT$("           ", 1 + LEN(printv$) - decimals - INSTR(printv$, "."))
  'END IF
  printv$ = RIGHT$("             " + printv$, length)
ELSE
  printv$ = RIGHT$("                  " + STR$(printv), length - decimals - 1)
  printv$ = printv$ + RIGHT$("                            ", decimals + 1)
END IF

END SUB

SUB pdint (X, DMW(), DMC(), RVT, DVC(), NDIV%, DM, PR(), XINT, DPR(), XI1, XI2, XNHT)
'  CALL pdint(X, XDMW(), XDMC(), RVT, DC(), NDIV%, DMRE, PR(), XINT, DPR(), XI1, XI2, XNHT)

'dimension DMW(720), DMC(720), DVC(720), PR(720), DPR(720)
'649 Subroutine for pressure drop integral calculation
DM = 0
XINT = 0
XI1 = 0
EX1 = 1 - XNHT
XI2 = 0
EX2 = 2 - XNHT
FOR I = 1 TO NDIV%
  DMX = DMC(I) - X * (DMW(I) / RVT + DMC(I))
  Y = ABS(DMX)
  DM = DM + Y
  A = DPR(I) * Y ^ EX1
  IF (DMX) < 0 THEN A = -A
  XI1 = XI1 + A
  XI2 = XI2 + Y ^ EX2

101 XINT = XINT + Y * DMX / PR(I) * DVC(I)
NEXT I
'xndiv = NDIV%  'XNDIV KOMT VERDER NIET VOOR

END SUB

FUNCTION plot (X(), H)
'function TO INTERPOLATE
'599 subroutine to .................
'dimension X(20)
'dim X(20)
M% = 10 * H
Z = H - M% / 10
M% = M% + 1
G1 = X(M%)
N% = M% + 1
G2 = X(N%)
plot = Z * G2 + (1 - Z) * G1

END FUNCTION

FUNCTION small (X(), NDIV%)
'subroutine to find smallest of a list
smallnow = X(1)
FOR I = 2 TO NDIV%
 IF smallnow > X(I) THEN
   smallnow = X(I)
 END IF
NEXT I
small = smallnow
END FUNCTION

SUB VOLC (DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(), COFI(), SALF(), CALF(), PDR, ALF)

'subroutine to list cold volumes and derivatives
'CALL VOLC(DALF, NF%, C(), CI(), DC(), DCI(), ZZC, NDIV%, SIFI(),
'COFI(), SALF(), CALF(), PDR, ALF)

NDIV1% = NDIV% + 1

FOR I = 1 TO NDIV1%
  COFI(I) = COS(ALF)
  SALF(I) = SIN(ALF)
  ALF = ALF + DALF
NEXT I

FOR I = 1 TO NDIV%
  CALF(I) = (SALF(I + 1) - SALF(I)) / DALF
NEXT I
CALF(NDIV1%) = CALF(1)'Richard has  NDIVI here
'CALF(NDIVI) = CALF(1)'Richard has NDIVI here
FOR I = 1 TO NDIV%
  SIFI(I) = SALF(I)
  SALF(I) = (SALF(I) + SALF(I + 1)) / 2
NEXT I
COFI(NDIV1%) = COFI(1)
SIFI(NDIV1%) = SIFI(1)
N% = NF% * 4
FOR I = 1 TO N%
  CRC = SQR(ZZC ^ 2 - CALF(I) ^ 2)
'  see note #1
IF (PDR) > 0 THEN
  C(I) = 1 - SALF(I) + CRC - ZZC
  CI(I) = 1 - SIFI(I) + CRC - ZZC
  DC(I) = -CALF(I) * (1 - SALF(I) / CRC)
  DCI(I) = -CALF(I) * (1 - SIFI(I) / CRC)
ELSE
  C(I) = 1 + SALF(I) - CRC + ZZC
  CI(I) = 1 + SIFI(I) - CRC + ZZC
  DC(I) = CALF(I) * (1 - SALF(I) / CRC)
  DCI(I) = CALF(I) * (1 - SIFI(I) / CRC)
END IF
NEXT I

END SUB

SUB VOLW (W(), WI(), DW(), DWI(), CFI, SFI, ZZW, NDIV%, SIFI(), COFI(), SALF(), CALF(), DALF)

'subroutine for pressure drop integral calculation

'subroutine to list hot volumes and derivates
'CALL VOLW(W(), WI(), DW(),DWI(),CFI,SFI,ZZW,NDIV%,SIFI(),COFI(),SALF(),CALF(),DALF)

SIFIP = SIFI(1) * CFI - COFI(1) * SFI

FOR I = 1 TO NDIV%
  SALF1 = SIFI(I + 1) * CFI - COFI(I + 1) * SFI
  SALFP = (SIFIP + SALF1) / 2
  CALFP = (SALF1 - SIFIP) / DALF
  CRW = SQR(ZZW ^ 2 - CALFP ^ 2)
  W(I) = 1 + SALFP - CRW + ZZW
  WI(I) = 1 + SIFIP - CRW + ZZW
  DW(I) = CALFP * (1 - SALFP / CRW)
  DWI(I) = CALFP * (1 - SIFIP / CRW)
  SIFIP = SALF1
NEXT I

END SUB

FUNCTION xlarge (X(), NDIV%)
'subroutine to find largest of a list
xlargenow = X(1)
FOR I = 2 TO NDIV%
 IF xlargenow < X(I) THEN
   xlargenow = X(I)
 END IF
NEXT I
xlarge = xlargenow
END FUNCTION

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