# Applied Python PRISM
# (PRISM) PaRametrIc System Model
#
# Written by Charlie Taylor <cet@appliedpython.com> 
# Oct,21 2005
from math import *

from prism.MassItem import MassItem
from prism.isp.cea import CEA_Isp
from prism.utils import Constants
from prism.Summary import Summary
from prism.props import Materials
from prism.isp import Nozzle_Eff
from prism.utils import mensuration
from prism.props import Inc_liquid
from prism.engines import ValveSolenoid
from prism.isp import separated_Cf
from prism.pov import POV_Items, POV_Basics

class Engine_FFC( MassItem ):
    
    def __init__(self, name="engine",  mass_lbm=0.0, oxName='N2O4', fuelName='MMH',
        cxw=1.25, Pc=150.0, Fvac=100.0, eps=50.0, mr=1.6, DtInp=None,
        CR=2.5, LoverDt=4.0, LchamMin=1.5, xlnOverLcham=0.5,
        etaERE=0.97, etaNoz=0.99, matlInj="SS", cxwInj=1.0, cxwValves=1.0, isBell=1, pcentBell=80.0,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, etaKinInp=1.0,
        calcEtaNoz=1, thkNozMin=0.02, matlNoz="Cb103", minBipropValveWt=0.2,
        inputIspDel=0, IspDel=300.0, valvesMassInput=None, Pamb=0.0):
        
        MassItem.__init__(self, name, type="inert")
        
        self.Number = Number # number of engines
        self.oxName = oxName
        self.fuelName = fuelName
        self.iprop = oxName + '/' + fuelName
        
        self.Fvac = Fvac
        self.DtInp = DtInp  # if DtInp is input, then Fvac is calculated
        self.Pc = Pc
        self.eps = eps
        self.mr = mr
        self.CR = CR
        self.xlnOverLcham = xlnOverLcham
        self.LoverDt = LoverDt
        self.LchamMin = LchamMin
        self.cxw = cxw
        self.cxwValves = cxwValves
        self.minBipropValveWt = minBipropValveWt
        self.etaERE = etaERE
        self.etaKinInp = etaKinInp
        self.Pamb = Pamb
        
        self.inputIspDel = inputIspDel
        self.IspDel = IspDel
        
        self.isBell = isBell
        self.pcentBell = pcentBell
        self.halfAngDeg = halfAngDeg
        self.etaNoz = etaNoz
        self.calcEtaNoz = calcEtaNoz
        
        self.matlInj = matlInj
        self.thkNozMin = thkNozMin
        self.cxwInj = cxwInj
        self.rhoInj, self.syInj, self.eInj, tmingInj = Materials.getMatlProps(matlInj)
        
        self.matlNoz = matlNoz
        self.rhoNoz, self.syNoz, self.eNoz, self.tmingNoz = Materials.getMatlProps(matlNoz)
        
        # assume storable liquids for the FFC engine
        self.FlObj = Inc_liquid.Inc_liquid(symbol=fuelName,T=None,P=Pc*1.5)
        self.OxObj = Inc_liquid.Inc_liquid(symbol=oxName,T=None,P=Pc*1.5)
        
        # assume identical ox and fuel valves
        self.valvesMassInput = valvesMassInput # if input, then do NOT calculate
        if valvesMassInput:
            self.Valves = None
        else:
            self.Valves = ValveSolenoid.ValveSolenoid( name='biprop valves', Number=2, cxw=cxwValves )
        
        
        self.ispObj = CEA_Isp.CEA_Isp( oxName=oxName, fuelName=fuelName, useFastLookup=useFastCEALookup ) # create isp calculating object
        self.reCalc()
        
    def getPOV_Item(self):
        # be sure to include pov_h, pov_w, and pov_d calcs in reCalc
        if hasattr( self, 'texture'):
            texture = self.texture
        else:
            
            texture = POV_Basics.Texture( colorName="Coral" )
        
        s = POV_Items.TCA_Bell( xlc=self.xlc , xln=self.xln, CR=self.CR, 
            Rt=self.Dt/2., eps=self.eps,  pcentBell=self.pcentBell,  texture=texture,
            isBell=self.isBell, halfAngDeg=self.halfAngDeg)
        
        return s
        
    def reCalc(self, autoCalc=1):
        self.autoCalc = autoCalc
        # set design variables
        
        if self.DtInp:
            self.Dt = self.DtInp
            self.At = pi * self.Dt**2/4.0 # used in calcNozzleEfficiency
        else:
            self.Dt = None
            self.At = None # used in calcNozzleEfficiency
        
        #self.IspODE,self.CstarODE,self.Tc = \
        #    self.ispObj.get_IvacCstrTc(Pc=self.Pc, MR=self.mr, eps=self.eps)

        self.IspODE,self.CstarODE,self.Tc, self.mw, self.gam = \
            self.ispObj.get_IvacCstrTc_ThtMwGam( Pc=self.Pc, MR=self.mr, eps=self.eps)

        self.PcOvPexit = self.ispObj.get_PcOvPe(Pc=self.Pc, MR=self.mr, eps=self.eps)
        self.Pexit = self.Pc / self.PcOvPexit
            
        if self.calcEtaNoz:
            if self.isBell:
                isConical=0
            else:
                isConical=1
            self.etaBL,self.etaDiv,self.etaKin, etaCf = \
                Nozzle_Eff.calcNozzleEfficiency(Pc=self.Pc, Fvac=self.Fvac, eps=self.eps, At=self.At,
                epsAtt=self.eps, isConical=isConical, pcentBell=self.pcentBell, halfAngleDeg=self.halfAngDeg,
                iprop=self.iprop, mr=self.mr, etaKinInp=self.etaKinInp,
                adjBL=1.0, adjKin=1.0, adjDiv=1.0, isRegenCham=0, isRegenNoz=0 )

            if self.etaKinInp != 1.0:
                self.etaKin = self.etaKinInp
                etaCf = self.etaBL * self.etaDiv * self.etaKin
                
            self.etaNoz = etaCf
        
        # if isp delivered is input, correct etaNoz to reflect that input
        if self.inputIspDel:
            self.Isp = self.IspDel
            self.effIsp = self.Isp / self.IspODE
            self.etaNoz = self.effIsp / self.etaERE
        else:
            self.effIsp = self.etaERE * self.etaNoz
            self.Isp = self.IspODE * self.effIsp
            
        # even if Isp is input, use CEA cstar to calculate throat diameter
        self.Cstar = self.CstarODE * self.etaERE
        
        if self.DtInp:
            self.wdotTot = self.Pc * self.At * 32.174 / self.Cstar
            self.Fvac = self.wdotTot * self.Isp
        else:
            self.wdotTot = self.Fvac / self.Isp
            self.At = self.Cstar* self.wdotTot / self.Pc / Constants.gc
            self.Dt = sqrt( self.At / pi ) * 2.0
            
        self.wdotOx  = self.wdotTot * self.mr / (1.0 + self.mr)
        self.wdotFl = self.wdotTot - self.wdotOx 
        
        self.volDotOx = self.wdotOx / self.OxObj.rho
        self.volDotFl = self.wdotFl / self.FlObj.rho
        
        if self.Valves:
            self.Valves.cuInchPerSec = max( self.volDotOx, self.volDotFl )
            self.Valves.reCalc()
            self.WtValves = self.Valves.mass_lbm
        else:
            self.WtValves = self.valvesMassInput
        
        ftPerSec = 30.0 # through valve
        dCharFl = sqrt( 4.0 * self.volDotFl/ pi / ftPerSec) 
        dCharOx = sqrt( 4.0 * self.volDotOx / pi / ftPerSec)
        self.dChar = max( dCharFl, dCharOx)
        
        # weight of two identical valves
        #self.WtValves = 2.0 * 0.268 * (self.dChar/1.279)**2
        #if self.WtValves < self.minBipropValveWt: # minimum of 45 (0.2 lbm) grams per valve (90 total)
        #    self.WtValves = self.minBipropValveWt
        #self.WtValves *= self.cxwValves
            
        self.WtMisc = 0.319 * (self.dChar/1.279)**2  # scale miscellaneous like valve
        #if self.WtMisc < 0.1:
        #    self.WtMisc = 0.1
        
        self.Dcham = self.Dt * sqrt( self.CR )
        self.Lcham = self.Dt * self.LoverDt
        if self.Lcham < self.LchamMin:
            self.Lcham = self.LchamMin
        self.Dexit = self.Dt * sqrt( self.eps )
        
        
        if self.isBell:
            self.Lnoz = (sqrt(self.eps)-1.0)*self.pcentBell*(self.Dt/2.0)/100.0/tan(15.0*pi/180.0)
            # curve fit of ratio to minimum length rao nozzle
            self.ratmlr = (self.pcentBell/100.0) * 1612.1/(self.eps + 1009.0)
            self.SAnoz = self.Dt**2/4.*(3.368*(self.eps+10.875)**1.2606 + \
                self.eps*(self.ratmlr-1.25)*10.75)
        else:
            self.Lnoz = (sqrt(self.eps)-1.0)*(self.Dt/2.0)/tan(self.halfAngDeg*pi/180.0)
            r1 = self.Dt/2.0
            r2 = self.Dexit/2.0
            self.SAnoz = pi * sqrt((r1-r2)**2 + self.Lnoz**2) * (r1+r2)
            
        
        self.thkCham = 0.06 * (self.Pc/1000.0) * (self.Dcham/1.48)
        if self.thkCham < self.thkNozMin: 
            self.thkCham = self.thkNozMin
            
        self.thkNoz = (self.thkCham*0.9 + 3.5*self.tmingNoz) / 4.5
        if self.thkNoz < self.thkNozMin: 
            self.thkNoz = self.thkNozMin
            
        self.WtNoz = self.thkNoz * self.SAnoz * self.rhoNoz
        
        self.xln = self.xlnOverLcham * self.Lcham
        self.xlc = self.Lcham - self.xln
        
        VolChamber = mensuration.cylVol( self.thkCham, self.Dcham, self.xlc ) + \
                     mensuration.coneVol( self.thkCham, self.Dcham, self.Dt, self.xln )
                     
        self.WtChamber = VolChamber * self.rhoNoz # chamber is made of same matl as nozzle
        
        thkMultAcc = 0.05 + 0.25/self.Dcham
        VolAccustic = mensuration.cylVol( thkMultAcc*self.Dcham, self.Dcham, self.Lcham/2.0 )
        self.WtAcoustic = VolAccustic * self.rhoNoz
        
        # the injector ht also gets smaller with increased Dcham
        htMultDCham = min(1.0, 1.5/self.Dcham)
        
        htInj =  self.Dcham*htMultDCham
        VolInj = mensuration.solidCylVol( self.Dcham*1.1 , htInj )
        self.WtInj = self.rhoInj * VolInj * self.cxwInj
           
        self.Lengine = self.Lnoz + self.Lcham + self.Dcham  # use Dcham as an inj face fwd length

        
        # add up parts
        self.mass_lbm = self.WtNoz + self.WtChamber + self.WtInj + self.WtAcoustic  + self.WtMisc
            
        self.mass_lbm = self.Number * ( self.mass_lbm * self.cxw   + self.WtValves )
        
        self.FtoW = self.Fvac * self.Number / self.mass_lbm
        
        self.pov_h = self.Lengine
        self.pov_w = max( self.Dexit, self.Dcham )
        self.pov_d = self.pov_w
        
        # get sea level performance
        #CfSL, CfOverCfvacSL, modeSL = separated_Cf.ambientCf(gam=self.gam, epsTot=self.eps, Pc=self.pc, Pamb=14.7)
        
        self.Cfvac = self.Isp * 32.174 / self.Cstar
        
        if self.Pamb > 0.0:
            CfOvCfvacAtEsep, CfOvCfvac, Cfsep, CfiVac, CfiAmbSimple, CfVac, epsSep, Psep = \
                separated_Cf.sepNozzleCf(self.gam, self.eps, self.Pc, self.Pamb)
    
            if self.Pexit > Psep:
                self.IspAmb = self.Isp - self.Cstar*self.Pamb*self.eps/self.Pc/32.174
                self.CfAmb = self.IspAmb * 32.174 / self.Cstar
            else:
                IspODEepsSep,CstarODE,Tc = \
                    self.ispObj.get_IvacCstrTc(Pc=self.Pc, MR=self.mr, eps=epsSep)
                    
                CfvacAtEsep = self.Cfvac * IspODEepsSep / self.IspODE
                
                self.CfAmb = CfvacAtEsep * CfOvCfvacAtEsep
                self.IspAmb = self.CfAmb * self.Cstar / 32.174
            
            # figure out mode of nozzle operation
            if self.Pexit > Psep:
                if self.Pexit > self.Pamb + 0.05:
                    self.mode = 'UnderExpanded (Pe=%g)'%self.Pexit
                elif self.Pexit < self.Pamb - 0.05:
                    self.mode = 'OverExpanded (Pe=%g)'%self.Pexit
                else:
                    self.mode = 'Pexit = %g'%self.Pexit
            else:
                self.mode = 'Separated (Psep=%g, epsSep=%g)'%(Psep,epsSep)
            
            # calc ambient thrust
            self.Famb = self.wdotTot * self.IspAmb
            
        else:
            self.IspAmb = self.Isp
            self.Famb = self.Fvac

        self.PfaceOvPc =(1. + .25/self.CR**1.8)  # multiply time Pc to get PcFace
        self.PcFace = self.Pc * self.PfaceOvPc


    def buildSummary(self):
        
        summList = []
        summ = Summary(  summName='Bipropellant Engine',
        componentName=self.name, mass_lbm=self.mass_lbm, type=self.type)
        
        summ.addAssumption( 'Propellants : ' + self.oxName + ' / ' + self.fuelName )
        summ.addAssumption( 'NASA CEA Code for ODE performance ')
        summ.addAssumption( 'Physical Weight Model ')
        summ.addAssumption( 'Injector Material is ' + self.matlInj)
        summ.addAssumption( 'Nozzle Material is ' + self.matlNoz)
        if self.isBell:
            summ.addAssumption( 'Bell Nozzle with Percent Bell = %g'%self.pcentBell)
        else:
            summ.addAssumption( 'Conical Nozzle with Half Angle = %g deg'%self.halfAngDeg )
        
        if self.Number>1:
            summ.addAssumption( 'Mass is for %i engines total'%self.Number )


        #summ.addInput(self, label='generic param', value=0.0, units='', format='%g')

        if self.DtInp:
            summ.addInput('Dt', self.Dt, 'in', '%.3f')
        else:
            summ.addInput('Fvac', self.Fvac, 'lbf', '%g')
        summ.addInput('Pc', self.Pc, 'psia', '%.1f')
        summ.addInput('eps', self.eps, '', '%g')
        if self.isBell:
            summ.addInput('%Bell', self.pcentBell, '%', '%.2f')
        else:
            summ.addInput('halfAngDeg', self.halfAngDeg, 'deg', '%.2f')
        summ.addInput('mr', self.mr, '', '%g')
        summ.addInput('CR', self.CR, '', '%g')
        summ.addInput('xlnOverLcham', self.xlnOverLcham, '', '%g')
        
        summ.addInput('LoverDt', self.LoverDt, '', '%g')
        summ.addInput('LchamMin',self.LchamMin, 'in', '%.3f')
        summ.addInput('cxwInj', self.cxwInj, '', '%g')
        summ.addInput('cxwValves', self.cxwValves, '', '%g')
        summ.addInput('cxw', self.cxw, '', '%g')
        summ.addInput('etaERE', self.etaERE, '', '%g')
        if not self.calcEtaNoz:
            summ.addInput('etaNoz', self.etaNoz, '', '%g')

        if self.Pamb > 0.0:
            summ.addInput('Pamb', self.Pamb, 'psia', '%g')
    


        # outputs
        if self.DtInp:
            summ.addOutput('Fvac', self.Fvac, 'lbf', '%g')
        summ.addOutput('Isp', self.Isp, 'sec', '%g')
        summ.addOutput('Cstar', self.Cstar, 'ft/sec', '%.1f')
        if self.inputIspDel:
            summ.addOutput('etaNoz', self.etaNoz, '', '%g')
        
        elif self.calcEtaNoz:
            summ.addOutput('etaBL', self.etaBL, '', '%g')
            summ.addOutput('etaDiv', self.etaDiv, '', '%g')
            summ.addOutput('etaKin', self.etaKin, '', '%g')
            summ.addOutput('etaNoz', self.etaNoz, '', '%g')
        summ.addOutput('effIsp', self.effIsp, '', '%g')
        
        summ.addOutput('IspODE', self.IspODE, 'sec', '%g')
        summ.addOutput('CstarODE', self.CstarODE, 'ft/sec', '%.1f')
        summ.addOutput('Tc', self.Tc, 'degR', '%.1f')
        
        summ.addOutput('PcFace', self.PcFace, 'psia', '%g')
        summ.addOutput('Pexit', self.Pexit, 'psia', '%g')

        summ.addOutput('wdotTot', self.wdotTot, 'lbm/sec', '%g')
        summ.addOutput('wdotOx ', self.wdotOx , 'lbm/sec', '%g')
        summ.addOutput('wdotFl', self.wdotFl, 'lbm/sec', '%g')

        summ.addOutput('rhoFl', self.FlObj.rho, 'lbm/cuin', '%.4f')
        summ.addOutput('rhoOx', self.OxObj.rho, 'lbm/cuin', '%.4f')

        summ.addOutput('volDotOx ', self.volDotOx , 'cuin/sec', '%g')
        summ.addOutput('volDotFl', self.volDotFl, 'cuin/sec', '%g')
        summ.addOutput('DFlow',self.dChar,'in','%.3f')
        
        summ.addOutput('At', self.At, 'sqin', '%g')
        
        if not self.DtInp:
            summ.addOutput('Dt', self.Dt, 'in', '%.3f')
            
        summ.addOutput('Dcham', self.Dcham, 'in', '%.3f')
        summ.addOutput('Dexit', self.Dexit, 'in', '%.3f')
        summ.addOutput('Lcham', self.Lcham, 'in', '%.3f')
        summ.addOutput('xlc', self.xlc, 'in', '%.3f')
        summ.addOutput('xln', self.xln, 'in', '%.3f')

        summ.addOutput('Lnoz', self.Lnoz, 'in', '%.3f')
        summ.addOutput('Lengine', self.Lengine, 'in', '%.3f')

        summ.addOutput('rhoInj', self.rhoInj, 'lbm/cuin', '%.3f')
        summ.addOutput('rhoNoz', self.rhoNoz, 'lbm/cuin', '%.3f')
        
        summ.addOutput('thkCham', self.thkCham, 'in', '%.3f')
        summ.addOutput('thkNoz', self.thkNoz, 'in', '%.3f')

        summ.addOutput('WtNoz', self.WtNoz, 'lbm', '%.3f')
        summ.addOutput('WtChamber', self.WtChamber, 'lbm', '%.3f')
        summ.addOutput('WtInj', self.WtInj, 'lbm', '%.3f')
        summ.addOutput('WtAcoustic',self.WtAcoustic, 'lbm', '%.3f')
        summ.addOutput('WtValves(2)',self.WtValves, 'lbm', '%.3f')
        summ.addOutput('WtMisc',self.WtMisc, 'lbm', '%.3f')
        if self.Number>1:
            summ.addOutput( 'wt/Engine', self.mass_lbm/self.Number, 'lbm', '%.3f' )
            summ.addOutput('F/W', self.Number * self.Fvac/self.mass_lbm, 'lbf/lbm', '%.3f')
        else:
            summ.addOutput('F/W',  self.Fvac/self.mass_lbm, 'lbf/lbm', '%.3f')
            
        if self.Pamb > 0.0:
            summ.addOutput('IspAmb', self.IspAmb, 'sec', '%g')
            summ.addOutput('Famb', self.Famb, 'lbf', '%g')
            summ.addOutput('nozzle condition', self.mode, '', '%s')
        
        summList.append( summ )
        
        # ======== now valves
        if self.Valves:
            summList.append( self.Valves.buildSummary() )
        
        return summList

if __name__ == "__main__":  #self test

    print "Actual IR&D Style Engine =", 2.071,"lbm (1.815 w/o valves)"
    
    h = Engine_FFC(name="Ref Axial Engine",  mass_lbm=0.0, oxName='N2O4', fuelName='MMH',
        cxw=1.27, Pc=500.0, Fvac=714.0, eps=17.0, mr=1.08, CR=2.06, LoverDt=1.2,
        etaERE=0.97, etaNoz=0.99, matlInj="Ti", cxwInj=1.0, isBell=1, pcentBell=84.58,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, calcEtaNoz=1,  matlNoz="Cb103")
    print "Calculated =   ",h.getMassStr()
    print "w/o valves= %.3f lbm"%(h.mass_lbm - h.WtValves,)
    print
    print h.getSummary()

    print

    print "Actual R-1E 25 lbf Engine =", 4.4,"lbm (2.4 w/o valves)"
    
    h = Engine_FFC(name="R-1E 25 lbf Engine",  mass_lbm=0.0, oxName='N2O4', fuelName='MMH',
        cxw=2.7, Pc=106.0, Fvac=25.0, eps=100.0, mr=1.65, CR=2.06, LoverDt=4.0,
        etaERE=0.90,  matlInj="SS", cxwInj=1.0, isBell=1, pcentBell=80.0,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, calcEtaNoz=1,  
        matlNoz="Cb103", thkNozMin=0.04, minBipropValveWt=2.0)
    print "Calculated =   ",h.getMassStr()
    print "w/o valves= %.3f lbm"%(h.mass_lbm - h.WtValves,)
    print
    print h.getSummary()
    
    print
    

    print "Actual R4-D Engine =", 7.5,"lbm (5.5 w/o valves)"
    
    h = Engine_FFC(name="R4-D Engine",  mass_lbm=0.0, oxName='N2O4', fuelName='MMH',
        cxw=2.63, Pc=108.0, Fvac=110.0, eps=44.0, mr=1.65, CR=2.06, LoverDt=2.5,
        etaERE=0.97, etaNoz=0.99, matlInj="SS", cxwInj=1.0, isBell=1, pcentBell=80.0,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, calcEtaNoz=1,  
        matlNoz="Cb103", thkNozMin=0.04, minBipropValveWt=2.0, inputIspDel=1, IspDel=300.0)
    print "Calculated =   ",h.getMassStr()
    print "w/o valves= %.3f lbm"%(h.mass_lbm - h.WtValves,)
    print
    print h.getSummary()

    

    print "MKV-R ACS Engine"
    
    h = Engine_FFC(name="ACS Engine",  mass_lbm=0.0, oxName='N2O4', fuelName='N2H4',
        cxw=2.63, Pc=800.0, Fvac=10.0, eps=15.0, mr=0.9, CR=2.5, LoverDt=2.5, LchamMin=1.2,
        etaERE=0.7, etaNoz=0.99, matlInj="SS", cxwInj=1.0, isBell=1, pcentBell=80.0,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, calcEtaNoz=1,  
        matlNoz="Cb103", thkNozMin=0.04, minBipropValveWt=2.0, inputIspDel=0, IspDel=300.0, valvesMassInput=0.8)
    print "Calculated =   ",h.getMassStr()
    print "w/o valves= %.3f lbm"%(h.mass_lbm - h.WtValves,)
    print
    print h.getSummary()

    
    h = Engine_FFC(name="Phoenix AKV Engine",  mass_lbm=0.0, oxName='N2O4', fuelName='M20_NH3',
        cxw=1.27, Pc=500.0, Fvac=714.0, eps=17.0, mr=1.5, CR=2.06, LoverDt=1.2,
        etaERE=0.97, etaNoz=0.99, matlInj="Ti", cxwInj=1.0, isBell=1, pcentBell=84.58,
        halfAngDeg=15.0, useFastCEALookup=0, Number=1, calcEtaNoz=1,  matlNoz="Cb103", Pamb=14.7)
    print "Calculated =   ",h.getMassStr()
    print "w/o valves= %.3f lbm"%(h.mass_lbm - h.WtValves,)
    print
    print h.getSummary()

    print