# calm.py = Compressed Air Launcher Modeler by AKBiocca 7/2002
print
print "CALM.py = Compressed Air Launcher Modeler in Python 7/2002 AKB ver 0.4"
print
from math import *

# some important constants
g = 32.2            # gravity ft / sec / sec

def diatoarea(dia):
    "circle diameter to area"
    return pi/4.*dia*dia

def cylvol(dia,len):
    "volume of cylinder"
    return diatoarea(dia)*len

# pvc pipe
pvc4id    = 4.03
pvc4od    = 4.50
pvc3id    = 3.07
pvc3od    = 3.50
pvc2p5id  = 2.45
pvc2p5od  = 2.88
pvc2id    = 2.07
pvc2od    = 2.38
pvc1p5id  = 1.61
pvc1p5od  = 1.90
pvc1p25id = 1.38
pvc1p25od = 1.66
pvc1id    = 1.05
pvc1od    = 1.32
pvcp75id  = 0.82
pvcp75od  = 1.05

## Schedule 40 PVC Data
##
## Nominal Pipe Size (in)
##       O.D.
##              Average I.D.
##                     Min. Wall
##                           Nominal Wt./Ft.
##                                  Max. W.P. PSI  
## 1/8    .405   .261  .068   .045  810  
## 1/4    .540   .354  .088   .081  780  
## 3/8    .675   .483  .091   .109  620  
## 1/2    .840   .608  .109   .161  600  
## 3/4   1.050   .810  .113   .214  480  
## 1     1.315  1.033  .133   .315  450  
## 1-1/4 1.660  1.364  .140   .426  370  
## 1-1/2 1.900  1.592  .145   .509  330  
## 2     2.375  2.049  .154   .682  280  
## 2-1/2 2.875  2.445  .203  1.076  300  
## 3     3.500  3.042  .216  1.409  260  
## 3-1/2 4.000  3.520  .226  1.697  240  
## 4     4.500  3.998  .237  2.006  220  
## 5     5.563  5.017  .258  2.726  190  
## 6     6.625  6.031  .280  3.535  180  
## 8     8.625  7.943  .322  5.305  160  
## 10   10.750  9.976  .365  7.532  140  
## 12   12.750 11.890  .406  9.949  130  


# problem constants
ballweight = 2.0/16.    # ball weight in lb (2.0 oz nom)
boredia  = pvc2p5id     # barrel bore diameter in inches (2.5)

launcherwt = 9.        # lbs

##borelen = 46.0      # barrel bore length in inches (with respect to ball travel)
##borelen  = 12.
boredrag = 10.5     # barrel drag in pounds
##chargepressure = 40 # stored air pressure in pounds per square inch
##
##chargedia = 2.44
##chargelen = 12.
##chargevolume = pi * chargedia * chargedia / 4. * chargelen + 10.

#chargevolume = 530. / 1.   # stored air volume in cubic inches
##deadvolume = 29. / 2.       # misc volume in valve etc in cubic in
#cv = 6.13*6.                 # valve/system conductance 1 psi gpm (.1-.3 typical)
# 5.8801 fits to 2346.1
# 6.13 fits 30-70 to xx
nvalves = 1                 # number of valves in parallel
interval = 0.00001          # modelling time interval in seconds

print "Ball Wt %.1f oz (nom 2)" % (ballweight * 16.)
#print "charge pressure %4.1f vol %4.1f" % (chargepressure,chargevolume)
#print "system Cv",cv

# derived constants
borearea = diatoarea(boredia)   #pi * boredia * boredia / 4. # bore area in sq inches
ballmass = ballweight / g

#print "bore len %4.1f area %3.1f volume %4.1f" % (borelen,borearea,borevolume)

#print "expansion ratio (chamber/rest) %4.2f" % expansionratio

temp = 70.                                  # temp deg f
scfmtocuinps = 12.*12.*12./60.              # scfm to cu in/s
    
def init():
    global supplypressure, supplypv, initialpv, borepv, boreposition, \
           velocity, time, iteration, borepressure, cvcf, chargevolume, \
           initialvolume, finalvolume
    borevolume = borelen * borearea
    initialvolume = deadvolume + chargevolume
    #finalvolume = initialvolume + borevolume
    #expansionratio = chargevolume / (deadvolume + borevolume+chargevolume)
    supplypressure = chargepressure + 14.7      # supply chamber pressure in psi

    supplypv = supplypressure * chargevolume
    initialpv = 14.7 * initialvolume + supplypv
    borepv = 14.7 * deadvolume                  # initial
    borepressure = 14.7                         # psia
    
    boreposition = boredia/2.                   # bore position inches (ball)
    velocity = 0.                               # velocity in fps
    time = 0.                                   # time in sec
    iteration = 0                               # iteration counter
    #borevolume = 0.                             # vol swept by ball

#print "scfm to cu in/s",scfmtocuinps
#print "iter    x   fps    sec   bp   cp"
#print "chgp iter    x   fps    sec   bp   cp"

def nextposition():
    global boreposition, velocity, time, iteration, borepressure, \
        supplypressure, borepv, supplypv
    
    borevolume = boreposition * borearea - boredia/2.       # cu in
    currentvolume = initialvolume + borevolume              # cu in
    equilibpressure = initialpv / currentvolume             # psi

    p1 = supplypressure
    p2 = borepressure
    if p1 < p2: p1 = p2

    t = temp + 460.                                         # rankin
    
    if p2 < (0.53 * p1):                                    # supersonic
        p2 = 0.53 * p1                                      # limits eff delP
        
    flowrate = 16.05 * cv * sqrt((p1*p1-p2*p2)/t)           # in SCFM

    flowvol = flowrate * scfmtocuinps * interval            # cu in @ stp

    flowpv = flowvol * 14.7 * nvalves                       # convert to pv

    borepv += flowpv                                        # implement flow
    supplypv -= flowpv

    borepressure = borepv / (borevolume + deadvolume)
    supplypressure = supplypv / chargevolume
    
    airforce = (borepressure-14.7) * borearea               # pounds
    
    if velocity > 0.:
        netforce = airforce - boredrag                      # pounds
    else:
        netforce = max(airforce-boredrag,0.)                # bore friction
    
    acceleration = netforce / ballmass                      # ft / sec / sec
    velocity += acceleration * interval                     # ft / sec
    boreposition += velocity*12 * interval                  # inches
    time += interval                                        # sec
    iteration += 1

def prnt():
    print "%4d %4.1f %5.1f %6.4f %4.1f %4.1f"%\
          (iteration,boreposition, velocity, time, \
           borepressure-14.7, supplypressure-14.7)

def prnt2():
    print "%3d %5d %4.1f %5.1f %6.4f %4.1f %4.1f"%\
          (chargepressure, iteration, boreposition, velocity, time, \
           borepressure-14.7, supplypressure-14.7)

def compute():
    "through the bore by time step"
    global velocity,recoile
    init()
    while boreposition < borelen and time < 0.1:
        if (iteration % 250) == 0:                  # print redux
            #prnt2()
            pass
        nextposition()
    #prnt2()
    if boreposition < borelen:                      # no escape
        velocity = 0.
    
    launchermass = launcherwt / g
    recoilv = ballmass * velocity / launchermass
    recoile = 0.5 * launchermass * recoilv * recoilv
    trajectory()

def computecoax():
    "through the bore by time step"
    global velocity,recoile
    modelpiston()

    if boreposition < borelen:                      # no escape
        velocity = 0.
    
    launchermass = launcherwt / g
    recoilv = ballmass * velocity / launchermass
    recoile = 0.5 * launchermass * recoilv * recoilv
    trajectory()

def findcv(p,v):
    #print
    #print "fit cv for given pressure and velocity"
    global cv,chargepressure
    cv = 0.
    chargepressure = p
    delcv = 10.
    stopdel = .01
    #print "chp vel    cv"
    while delcv >= stopdel:
        cv += delcv
        compute()
        if velocity > v:
            cv -= delcv
            delcv *= 0.1
    compute()
    print "%3d %3d %5.2f"%(p,velocity,cv)



def modelpiston():
    "model coaxial piston system"
    global velocity, borepressure, boreposition, supplypressure, time, recoile

    pcv = 2.8       # piston volume exhaust valve cv
    pistonwt = 8    # ozs

    #boredia = pvc2p5id

    borearea = diatoarea(boredia)
    ballvolume = 4./3.*pi * boredia*boredia*boredia / 8.
    
    pistonvolume = cylvol(chargedia,2)+cylvol(pvc1p5id,8)
    supplypressure = chargepressure + 14.7      # supply chamber pressure in psi
    pistonpressure = supplypressure
    pistonpv = pistonpressure*pistonvolume
    borearea = diatoarea(boredia)
    pistonarea = diatoarea(pvc4id) - borearea   # ring outside bore
    borecircum = pi * boredia
    #print "pistonvolume",pistonvolume,"pistonarea",pistonarea,"borearea",borearea

    pistontravel = 0.75             # in
    pistonmass = pistonwt/16./g     # lb/g
    pistonpos = 0.                  # in
    pistonvel = 0.                  # fps
    pistonaccel = 0.                # fpsps
    t = temp + 460.                 # rankin

    chargearea = diatoarea(chargedia) - diatoarea(innerdia)
    chargevolume = chargelen * chargearea
    #borevolume = borelen * borearea
   # initialvolume = deadvolume + chargevolume
    #finalvolume = initialvolume + borevolume
    #expansionratio = chargevolume / (deadvolume + borevolume+chargevolume)
    
    supplypv = supplypressure * chargevolume
    #initialpv = 14.7 * initialvolume + supplypv

    boreposition = boredia/2. + initialpos      # bore position inches (ball)
    deadvolume = cylvol(boredia,boreposition)-ballvolume/2.
    borepressure = 14.7                         # psia
    borepv = borepressure * deadvolume          # initial

    velocity = 0.                               # velocity in fps
    ecv = 0.
    time = 0.
    iteration = 0

    if details:
        print
        print "  time ppsi     pf  ppos   pvl vel  ecv  pos   bp   cp"

    while boreposition < borelen and time < 0.1:
        p1 = pistonpressure                                     # vent piston
        p2 = 14.7
        if p1 < p2: p1 = p2
        if p2 < (0.53 * p1):                                    # supersonic
           p2 = 0.53 * p1                                       # limits eff delP 
        flowrate = 16.05 * pcv * sqrt((p1*p1-p2*p2)/t)          # in SCFM
        flowvol = flowrate * scfmtocuinps * interval            # cu in @ stp
        flowpv = flowvol * 14.7                                 # to pv
        pistonpv -= flowpv                                      # implement flow
        pistonswept = cylvol(chargedia,pistonpos)
        pistonpressure = pistonpv / (pistonvolume - pistonswept)
        pistonforce = (supplypressure-pistonpressure)*pistonarea - \
                      (pistonpressure-borepressure) * borearea
        if pistonpos <= 0:
            pistonvel = 0.
            pistonpos = 0.
        if pistonpos > pistontravel:
            pistonvel = 0.
            pistonpos = pistontravel

        pistonaccel = pistonforce/pistonmass    # a=f/m
        pistonvel += pistonaccel * interval     # fps
        pistonpos += pistonvel * interval * 12  # in

        effopen = pistonpos*borecircum          # sq in open cylinder
        ecv = cvpi * effopen                    # effective cv wild guess

        p1 = supplypressure                             # calc flow into bore
        p2 = borepressure
        if p1 < p2: p1 = p2
        if p2 < (0.53 * p1):                                    # supersonic
           p2 = 0.53 * p1                                       # limits eff delP
        flowrate = 16.05 * ecv * sqrt((p1*p1-p2*p2)/t)          # in SCFM
        flowvol = flowrate * scfmtocuinps * interval            # cu in @ stp
        flowpv = flowvol * 14.7                                 # convert to pv
        borepv += flowpv                                        # implement flow
        supplypv -= flowpv

        borevolume = boreposition * borearea - ballvolume/2.
        borepressure = borepv / borevolume                      # psi
        supplypressure = supplypv / chargevolume                # psi
        
        airforce = (borepressure-14.7) * borearea               # pounds
    
        if velocity > 0.:
            netforce = airforce - boredrag                      # pounds
        else:
            netforce = max(airforce-boredrag,0.)                # bore friction
        
        acceleration = netforce / ballmass                      # ft / sec / sec
        velocity += acceleration * interval                     # ft / sec
        boreposition += velocity*12 * interval                  # inches

        time += interval                                        # sec
        iteration += 1

        if (iteration % 100) == 0 and (pistonpos > 0) and details:
            pass            
            print "%6.4f %4.1f %6.1f %4.3f %5.1f %3d %4.1f %4.1f %4.1f %4.1f"%\
                  (time,pistonpressure-14.7,pistonforce,\
                                      pistonpos,pistonvel,velocity,ecv,\
                   boreposition,borepressure-14.7,supplypressure-14.7)

    if details:
        print "%6.4f %4.1f %6.1f %4.3f %5.1f %3d %4.1f %4.1f %4.1f %4.1f"%\
                  (time,pistonpressure-14.7,pistonforce,\
                                      pistonpos,pistonvel,velocity,ecv,\
                   boreposition,borepressure-14.7,supplypressure-14.7)

# test
##borelen=20
##deadvolume = cylvol(pvc2p5id,1.5)
##chargelen = borelen-2
##chargevolume = cylvol(pvc4id,chargelen)-cylvol(pvc2p5od,chargelen)
##chargepressure = 40
##modelpiston()

def piston():
    "simple force piston travel model"
    travel = 0.75                   # in
    pistonmass = pistonwt/16./g     # lb/g
    pistonpos = 0.                  # in
    pistonvel = 0.                  # fps
    pistonaccel = 0.                # fpsps
    iter = 0
    time = 0.
    
    while pistonpos < travel:
        time += interval
        iter += 1
        pistonforce = 40. * 6.                  # f=pressure*area
        pistonaccel = pistonforce/pistonmass    # a=f/m
        pistonvel += pistonaccel * interval     # fps
        pistonpos += pistonvel * interval * 12  # in
        
        if (iter % 100) == 0:
            pass
            #print "%3.1f %3.1f %5.3f %4.1f"%(pistonwt,time*1000.,pistonpos,pistonvel)
    
    print "%3d %3.1f %5.3f %4.1f"%(pistonwt,time*1000.,pistonpos,pistonvel)

def pistons():
    global pistonwt
    print
    print "Piston Motion"
    print
    print "pwt  ms   pos  fps"
    for pistonwt in (3,4,5,6,12,16,20):    # oz
        piston()

#pistons()

def trajectory():
    "simple vertical trajectory based on D*v^2 drag and terminal v"
    global peak,ttime
    iter = 0
    tinterval = interval * 100.
    vposition = 0.                  # vertical position
    vvelocity = velocity            # vertical velocity
    terminalv = 70./60./60.*5280.   # 70 mph terminal from www.phys.washington.edu/~young208A/ball-air.html
##    print "terminal velocity %d fps"%terminalv
##    criticalspeed = 53.             # m/s if smooth (relovent?)
    ttime = 0.                      # trajectory time
    ballarea = diatoarea(2.56)
    peak = 0.
##    airdensity = 1.
##    #Cd = 0.5
##    Cd = 0.45                       # [Daish]

    # terminal velocity: ballwt = drag: 2./16. = D * tv * tv
    D = 2./16./(terminalv*terminalv)
##    print "D",D
##    print
##    print " iv time  pos  vel"
    while vposition >= 0. :
        iter += 1
        ttime += tinterval
        vposition += vvelocity * tinterval
        if peak < vposition:
            peak = vposition
        vacceleration = 0.
##        airviscosity = 1.8e-5 # kg/(s*m)
##        reynoldsnum = airdensity * velocity * balldiameter / airviscosity
##        if reynoldsnum < 1.4e5:     # [de Mestre]
##            Cd = 0.45
##        else:
##            Cd = 0.2
        #drag = 0.5 * airdensity * vvelocity * vvelocity * Cd * ballarea
            
        # perhaps use the fact that dragforce = D * v * v, and
        # the terminal velocity to make a simple calculation of D
        drag = D * vvelocity * vvelocity
        if vvelocity > 0. :
            drag = -drag
        drag *= 1.0  # adjust for altitude 0.7 for 8k?
        vacceleration = drag/ballmass
        vvelocity += (vacceleration - g) * tinterval
        
#        if (iter % 50000) == 0 :
#            print "%5.3f %4d %4d"%(ttime,vposition,vvelocity)
            
    #print "%3d %6.3f %4d %4d %4d"%(velocity,ttime,vposition,vvelocity,peak)

def trajectories():
    "trajectory table vs velocity"
    global velocity
    print
    print " iv   time  pos  vel peak"
    for velocity in range(100,451,50):
        trajectory()

#trajectories()
        

def fitericscvs():
    "calc cvs fm erics data"
    global borelen,deadvolume,chargevolume,nvalves
    borelen = 47.
    deadvolume = 29
    chargevolume = 530./2.
    nvalves = 2
    print "chp vel    cv"
    for chargepressure,measuredv in ((20,78),(30,173),(40,237),(50,284),\
                                     (60,323),(70,362)):
        findcv(chargepressure,measuredv)
        
#fitericscvs()
        # 20-5.63, 30-5.98, 40-6.41, 50-6.67, 60-6.89, 70-7.30

def computerange():
    "check against erics data"
    global chargepressure,errsq
    errsq = 0.
    for chargepressure,measuredv in ((20,78),(30,173),(40,237),(50,284),\
                                     (60,323),(70,362)):
        if (measuredv > 10):
            compute()
            err1 = (measuredv - velocity)
            errsq += err1 * err1

#computerange()

def showrange():
    global cvcf,cv
    for i in range(-12,12):
        cv = (6.13 + (i) * 0.001)*2
        computerange()
        print "err %8.3f, cv %8.4f"%(errsq,cv/2.)
        #print "err %8.3f, cvcf %8.4f"%(errsq,cvcf)

#showrange()

def explore1():
    ""
    global chargepressure,chargelen,borelen
    print "chgp chgv bbl fps sec xbp xcp"
    for chargepressure in range(10,101,5):
        for chargelen in range(16,17,2):
            for borelen in range(chargelen-2,chargelen-1,2):
                compute()
                if velocity > 10.:
                    print "%3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)

#explore1()

def erics2():
    print
    print "erics mark 2 dual straight valve u chamber long barrel"
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume,nvalves
    # eric's data
    # psi fps
    # 20  78
    # 30 173
    # 40 237  
    # 50 284
    # 60 323
    # 70 362  
    #
    # bore len 47" - 1"
    # bore drag 10.5 lb
    # chamber vol 529.5 cu in
    # misc vol 29 cu in
    
    # 20-5.63, 30-5.98, 40-6.41, 50-6.67, 60-6.89, 70-7.30
    cv = "variable"
    nvalves = 2
    deadvolume = 29.
    print "system Cv",cv,"for each of",nvalves,"valves"
    print
    print "chp  cv bbl  fps   sec  xbp  xcp  ht tofl"
    print
    for chargepressure,cv in ((20,5.63),(30,5.98),(40,6.41),(50,6.67),(60,6.89),(70,7.30)):
        for chargelenx in (47,):
            for borelen in range(47,47+1,2):
                chargedia = pvc3id
                chargearea = diatoarea(chargedia)
                chargelen = 529.5 / chargearea / 2.
                chargevolume = chargearea * chargelen
                #print chargelen,chargevolume
                compute()
                if velocity > 10.:
                    print "%3d %3.1f %3d %4.0f %5.3f %4.1f %4.1f %3d %4.1f"%\
                          (chargepressure, cv, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7, peak, ttime)

#erics2()


# TAC Tennis ball Air Cannon

# TBL tennis ball launcher

# AC boredia-chambercuin

def sltest1():
    print
    print "straight line system"
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume
    cv = 24.
    print "system Cv",cv
    valvelen = 15      # valve and misc
    deadvolume = cylvol(pvc1p5id,4) #more helps?
    chargedia = pvc4id
    print "chp oal chl bbl  fps   sec  xbp  xcp"
    for chargepressure in (20,40,60):
        for oal in (48,):
            for borelen in range(12,24,2):
                chargelen = oal - borelen - valvelen
                chargevolume = cylvol(chargedia,chargelen)
                compute()
                if velocity > 30.:
                    print "%3d %3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, oal, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)
sltest1()

def ous():
    print
    print "over and under straight valve short bore"
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume
    cv = 6.
    print "system Cv",cv
    deadvolume = 3. #more helps?
    chargedia = pvc3id
    print "chp chl bbl  fps   sec  xbp  xcp"
    for chargepressure in (20,30,40,50,60):
        for oal in (36,):
            chargelen = oal - .5 - 2. - 3.
            for borelen in range(18,oal - 2. - 6. - 3.,2):
                chargearea = pi/4. * (chargedia*chargedia)
                chargevolume = chargearea * chargelen
                compute()
                if velocity > 10.:
                    print "%3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)
#ous()

def ouu():
    print
    print "over and under u valve"
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume,nvalves
    cv = 3.
    nvalves = 1
    print "system Cv",cv
    deadvolume = 40.
    print "chp chl bbl  fps   sec  xbp  xcp"
    for chargepressure in (20,30,40,50,60,80,100):
        for oal in (36,):#48,47+3+1+6):
            chargelen = oal - 6. - 1. - 2.
            for borelen in range(18,chargelen,4):
                chargedia = pvc3id
                chargearea = diatoarea(chargedia)
                chargevolume = chargearea * chargelen
                compute()
                if velocity > 50.:
                    print "%3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)
#ouu()

uvalvelen = 6.5 # u valve plus pipe

def alansouudb():
    print
    print "alans over and under u valve dual dia bore"
    print
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume,nvalves
    cv = 3.+2
    nvalves = 1
    chargedia = pvc3id
    chargearea = diatoarea(chargedia)
    deaddia = pvc2id
    print "system Cv",cv,"chargearea %.1f"%chargearea,"borearea %.1f"%borearea
    print
    print "chp oal chl bbl  dl  fps   sec  xbp  xcp  ht tofl"
    print
    for chargepressure in (30,40,50,60,70,80):
        #print
        for oal in (36,):#48,47+3+1+6):
            chargelen = oal - uvalvelen - 1. - 2. - 2.
            for borelen in range(20,chargelen+3-4,20):
                chargevolume = chargearea * chargelen
                deadlen = oal - uvalvelen - borelen -1 -1 -1
                deadvolume = cylvol(deaddia, deadlen) + 6.
                compute()
                if velocity > 50.:
                    print "%3d %3d %3d %3d %3d %4.0f %5.3f %4.1f %4.1f %3d %4.1f"%\
                          (chargepressure, oal, chargelen, borelen, deadlen,\
                           velocity, time, \
                           borepressure-14.7, supplypressure-14.7, peak, ttime)
#alansouudb()

def dualu():
    print
    print "dual u valve"
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume,nvalves
    
    cv = 3.
    nvalves = 2
    deadvolume = 30.
    oal = 36.
    chargelen = oal - 6. -1. - 2.
    chargedia = 2.98
    chargearea = pi/4. * (chargedia*chargedia)
    chargevolume = chargearea * chargelen * 2.
    print "system Cv",cv,"Chargevolume %d"%chargevolume
    #print chargelen,chargevolume
    print "chp chl bbl  fps   sec  xbp  xcp"
    for chargepressure in (20,30,40,50):
        for chargelenx in (47,): 
            for borelen in (12,18,20,22):
                compute()
                if velocity > 10.:
                    print "%3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)
#dualu()

def ouds():
    print
    print "over and under dual straight valve short bore"
    print
    global chargepressure,chargelen,borelen,cv,deadvolume,chargevolume,nvalves
    cv = 3.
    nvalves = 2
    print "system Cv",cv,"for each of",nvalves,"valves"
    deadvolume = cylvol(pvc1id,14)
    chargedia = pvc3id
    print
    print "chp chl bbl  fps   sec  xbp  xcp"
    for chargepressure in (20,30,40,50,60,70,80,100):
        for oal in (36,):
            chargelen = 24. #oal - .5 - 2. - 3.
            for borelen in (20,):  #range(14,oal - 2. - 6. - 3. - 5.,6):
                chargevolume = cylvol(pvc3id,chargelen) + cylvol(pvc1id,(4*5+4*2))
                #print chargevolume
                compute()
                if velocity > 50.:
                    print "%3d %3d %3d %4.0f %5.3f %4.1f %4.1f"%\
                          (chargepressure, chargelen, borelen,\
                           velocity, time, borepressure-14.7, supplypressure-14.7)
#ouds()

def minicoaxial():
    
    global chargepressure,borelen,deadvolume,chargevolume,cv,launcherwt,nvalves
    cv = 6.*4
    nvalves = 1
    deadvolume = 4.
    handlelen = 4.
    launcherwt = 4.
    pistonlen = 2.  # incl travel es stopper
    chargedia = pvc4id
    innerdia = pvc2p5od
    chargearea = diatoarea(chargedia) - diatoarea(innerdia)
    print
    print "mini coaxial design, %4.2f in coaxial chamber, %4.2f in bore"%\
          (chargedia,boredia)
    print
    print "system Cv",cv,"chargearea %.1f"%chargearea,"borearea %.1f"%borearea
    print
    print "chp oal chl bbl  fps   sec  rec  xbp  xcp  ht tofl"
    for chargepressure in (15,20,30,40,60):
        for borelen in range(16,17,2):
            oal = borelen + handlelen + pistonlen + 2.
            for chargelen in range(borelen-2,borelen-1,1):
                chargevolume = chargearea * chargelen
                compute()
                if velocity > 50.:
                    print "%3d %3d %3d %3d %4.0f %5.3f %4.1f %4.1f %4.1f %3d %4.1f"%\
                          (chargepressure, oal, chargelen, borelen,\
                           velocity, time, recoile,\
                           borepressure-14.7, supplypressure-14.7,peak,ttime)

#minicoaxial()

def simplecvcoaxial():
    
    global chargepressure,borelen,deadvolume,chargevolume,cv,nvalves
    cv = 6./(pi*1*1/4)*4.7
    nvalves = 1
    deadvolume = cylvol(pvc2p5id,1)
    handlelen = .25+4.5+2.75+2+3.25
    pistonlen = .125+.125+.75+1.25  # incl travel es stopper
    chargedia = pvc4id
    innerdia = pvc2p5od
    chargearea = diatoarea(chargedia) - diatoarea(innerdia)
    launcherwt = 11.
    print
    print "simple cv coaxial design"
    print "%4.2f in coaxial chamber, %4.2f bore %.1f handle"%\
          (chargedia,boredia,handlelen)
    print
    print "system Cv",cv,"chargearea %.1f"%chargearea,"borearea %.1f"%borearea
    print
    print "chp  oal cv chl bbl  fps   sec  rec  xbp  xcp  ht tofl"
    print
    for chargepressure in (15,20,30,40,60,80):
        #print
        for xoal in (36,):
            #borelen = oal - handlelen - 2. - pistonlen
            for chargelen in range(17,19,2):
              #for cv in (2,4,8,16,32,64):
                chargevolume = chargearea * chargelen
                borelen = chargelen + 2 # effective due to start pos
                oal = 0.125+borelen+pistonlen+1+.25+handlelen
                compute()
                if velocity > 50.:
                    print "%3d %4.1f %2d %3d %3d %4.0f %5.3f %4.1f %4.1f %4.1f %3d %4.1f"%\
                          (chargepressure, oal, cv, chargelen, borelen,\
                           velocity, time, recoile,\
                           borepressure-14.7, supplypressure-14.7,peak,ttime)

#simplecvcoaxial()
                    
def coaxial():
    "eval coaxial design w piston model"
    global chargepressure,boredia,borelen,chargelen,chargedia,innerdia,cvpi,\
           initialpos,details

    details = 0
    cvpi = 6    # cv per square inch opening
    boredia = pvc2p5id
    chargedia = pvc4id
    innerdia = pvc2p5od

    initialpos = 1  # inches bore start to rear ball

    handlelen = .25+4.5+2.75+2+3.25
    pistonlen = .125+.125+.75+1.25  # incl travel es stopper
    
    #chargearea = diatoarea(chargedia) - diatoarea(innerdia)
    launcherwt = 11.                # lbs
    print
    print "coaxial design, piston modeled"
    print "%4.2f in coaxial chamber, %4.2f bore %.1f handle"%\
          (chargedia,boredia,handlelen)
    print
    print "chp  oal chl bbl ip  fps   sec  rec  xbp  xcp  ht tofl"
    #print
    for chargepressure in (15,20,30,40,60,80):
        print
        for initialpos in range(1,12,8):
            #borelen = oal - handlelen - 2. - pistonlen
            for chargelen in range(17,21,12):
              #for cv in (2,4,8,16,32,64):
                #chargevolume = chargearea * chargelen
                borelen = chargelen + 2
                oal = 0.125+borelen+pistonlen+1+.25+handlelen
                computecoax()
                if velocity > -50.:
                    print "%3d %4.1f %3d %3d %2d %4.0f %5.3f %4.1f %4.1f %4.1f %3d %4.1f"%\
                          (chargepressure, oal, chargelen, borelen,\
                           initialpos, velocity, time, recoile,\
                           borepressure-14.7, supplypressure-14.7,peak,ttime)

#coaxial()

#eof