# Applied Python PRISM
# (PRISM) PaRametrIc System Model
#
# Written by Charlie Taylor <cet@appliedpython.com> 
# Apr 21, 2008

from prism.props.refprop7.n_dll_fluid import n_fluid
from prism.isp.tdk import diskCache

from prism.MassItem import MassItem
from math import *
from prism.Summary import Summary
from prism.press.TankPress import TankPress

#C THIS ROUTINE CALCULATES THE COLD GAS Nitrogen REQUIRED TO PRESSURIZE
#C A STORABLE PROPELLANT TANK
#C

class PressurantInteg( MassItem ):
    
    def __init__(self, name="tank",  mass_lbm=0.0, gas='HE',
        timeProfileL=None, pcentLiqExpelledL=None,
        VpropTnk=1000.0,PGasTnkMEOP=5000.0,PpropNom=350.0,
        PfinGasOvPnom=1.1, 
        tAction=100.0,TminR=500.0,TmaxR=550.0, ullageFrac=0.03,
        PVoW_Bottle=500000., PVoW_Tank=100000.,
        AccGees=1.0,
        Nbottle=1, ellBottle=1.0, LcylOvDBottle=0.0, Cp_effBottle=0.15, # Cp Ti=.125, Al=.2, Monel=.1
        Ntank=2, ellTank=1.414, LcylOvDTank=1.0, Cp_effTank=0.15,
        CdARegMax=None, dPregulator=25.0, NtimeSteps=400, heatExchangerTout=None,
        WtHeLLACS=0.0, QexternalIntoBottle=0.0, velMultTank=1.0, adiabatic=None, TbottleMatlConst=None,
        useDBruns=1, fracExpel=None,  Tullage=None):
        
        MassItem.__init__(self, name, type="pressurant")
        

        self.timeProfileL=timeProfileL
        self.pcentLiqExpelledL=pcentLiqExpelledL

        self.VpropTnk=VpropTnk
        self.ullageFrac = ullageFrac
        self.PGasTnkMEOP=PGasTnkMEOP
        self.PpropNom=PpropNom
        self.PfinGasOvPnom=PfinGasOvPnom
        self.tAction=tAction
        self.TminR=TminR
        self.TmaxR=TmaxR
        self.gas = gas
        self.gasObj = n_fluid(gas,T=530.0,P=1000.0)

        
        self.PVoW_Bottle = PVoW_Bottle
        self.PVoW_Tank = PVoW_Tank
        
        self.AccGees = AccGees
        self.Nbottle = Nbottle
        self.ellBottle = ellBottle
        self.LcylOvDBottle = LcylOvDBottle
        self.Cp_effBottle = Cp_effBottle
        self.Ntank = Ntank
        self.ellTank = ellTank
        self.LcylOvDTank = LcylOvDTank
        self.Cp_effTank = Cp_effTank
        self.CdARegMax = CdARegMax
        self.dPregulator = dPregulator
        self.NtimeSteps = NtimeSteps
        self.heatExchangerTout = heatExchangerTout
        
        self.WtHeLLACS = WtHeLLACS
        
        self.QexternalIntoBottle = QexternalIntoBottle
        self.velMultTank = velMultTank
        self.adiabatic = adiabatic
        
        self.TbottleMatlConst = TbottleMatlConst
        
        self.fracExpel = fracExpel
        self.Tullage = Tullage  # if None, then use TminR as Tullage

        self.useDBruns = useDBruns

        self.reCalc()

    def calcPolytropicValues(self):
        
        Tbi = self.TminR
        Tui = self.TminR
        Tuf = self.TfinalPropGas
        Tbf = self.TfinalGasBot
        
        Pbi = self.PbottleCold
        Pbf = self.PfinalGasBot


        R = log(Tuf/Tui) / log(Pbf/Pbi)
        try:
            gamPolyTank = 1./(1.-R)
        except:
            gamPolyTank = 10.0  # make a silly big number
        #print 'gamPolyTank=',gamPolyTank
        
        R = log(Tbf/Tbi) / log(Pbf/Pbi)
        try:
            gamPolyBottle = 1./(1.-R)
        except:
            gamPolyBottle = 10.0  # make a silly big number
        #print 'gamPolyBottle=',gamPolyBottle
        
        fracIsenBottle = (gamPolyBottle-1.0) / (self.gammaIsen-1.0)
        fracIsenTank = (gamPolyTank-1.0) / (self.gammaIsen-1.0)
        
        return gamPolyBottle, fracIsenBottle, gamPolyTank, fracIsenTank
        
        
    def buildSummary(self):
        
        summ = Summary(  summName='%s Pressurant'%self.gasObj.name,
        componentName=self.name, mass_lbm=self.mass_lbm, type=self.type)
        
        summ.addAssumption( 'Integrated Heat Transfer')
        if self.TbottleMatlConst:
            summ.addAssumption( 'Bottle Material Temperature is Constant = %g degR'%self.TbottleMatlConst)
            
        if self.fracExpel != None:
            summ.addAssumption( 'Blow down after = %g %% Expulsion'%(self.fracExpel*100.,))
			
        if self.heatExchangerTout:
            summ.addAssumption( 'Outlet gas heated to %g degR via Heat Exchanger.'%self.heatExchangerTout )
        
        # add inputs
        summ.addInput( 'VpropTnk', self.VpropTnk, 'cuin', '%g' )
        summ.addInput( 'PGasTnkMEOP', self.PGasTnkMEOP, 'psia', '%g' )
        summ.addInput( 'PpropNom', self.PpropNom, 'psia', '%g' )
        summ.addInput( 'PfinGasOvPnom', self.PfinGasOvPnom, '', '%g' )
        summ.addInput( 'tAction', self.tAction, 'sec', '%g' )
        summ.addInput( 'TminR', self.TminR, 'degR', '%.1f' )
        summ.addInput( 'TmaxR', self.TmaxR, 'degR', '%.1f' )
        if self.fracExpel != None:
            summ.addInput( 'fracExpel', self.fracExpel, '', '%.4f' )
            summ.addInput( 'tRegulated', self.tAction*self.fracExpel, 'sec', '%g' )
        if self.Tullage:
            summ.addInput( 'Tullage', self.Tullage, 'degR', '%.1f' )
        
        # add outputs
        summ.addOutput( 'WGasTotal', self.WtGasTotal, 'lbm', '%.3f' )
        summ.addOutput( 'WtGasResid', self.WtGasResid, 'lbm', '%.3f' )
        if self.WtHeLLACS > 0.0:
            summ.addOutput( 'WtHeLLACS', self.WtHeLLACS, 'lbm', '%.3f' )
        summ.addOutput( 'WGasNon-Resid', self.WGasExpended, 'lbm', '%.3f' )
        summ.addOutput( 'Vbottle', self.Vbottle, 'cuin', '%g' )
        summ.addOutput( 'wdotGas', self.wdotGas, 'lbm/sec', '%g' )
        
        summ.addOutput( 'compressInit', self.compressInit, '', '%g' )
        summ.addOutput( 'PbottleCold', self.PbottleCold, 'psia', '%g' )
        
        summ.addOutput( 'densInitBot', self.densInit, 'lbm/cuft', '%g' )
        summ.addOutput( 'densFinalBot', self.densFinal, 'lbm/cuft', '%g' )
        summ.addOutput( 'densFinalProp', self.densFinalProp, 'lbm/cuft', '%g' )

        summ.addOutput( 'gammaIsen', self.gammaIsen, '', '%g' )
        
        #gamPolyBottle, fracIsenBottle, gamPolyTank, fracIsenTank = self.calcPolytropicValues()
        summ.addOutput( 'gamPolyBottle', self.gamPolyBottle, '', '%g' )
        summ.addOutput( 'gamPolyTank', self.gamPolyTank, '', '%g' )
        summ.addOutput( 'fracIsenBottle', self.fracIsenBottle, '', '%g' )
        summ.addOutput( 'fracIsenTank', self.fracIsenTank, '', '%g' )
        
        summ.addOutput( 'TfinalPropGas', self.TfinalPropGas, 'degR', '%.1f' )
        summ.addOutput( 'TfinalGasBot', self.TfinalGasBot, 'degR', '%.1f' )
        summ.addOutput( 'PfinalGasBot', self.PfinalGasBot, 'psia', '%.1f' )
        
        if self.fracExpel != None:
            summ.addOutput( 'PfinalBlowdown', self.PfinalBlowdown, 'psia', '%.1f' )
            summ.addOutput( 'timeExpel', self.timeExpel, 'sec', '%g' )
            

        return summ


    def reCalc(self):
        if self.useDBruns:
            runSig = diskCache.buildRunSignature( self.gas,
                self.VpropTnk,self.PGasTnkMEOP,self.PpropNom,
                self.PfinGasOvPnom, self.WtHeLLACS,
                self.tAction,self.TminR,self.TmaxR, self.ullageFrac,
                self.PVoW_Bottle, self.PVoW_Tank,
                self.AccGees,
                self.Nbottle, self.ellBottle, self.LcylOvDBottle, self.Cp_effBottle,
                self.Ntank, self.ellTank, self.LcylOvDTank, self.Cp_effTank,
                self.CdARegMax, self.dPregulator, self.NtimeSteps, self.heatExchangerTout,
                self.WtHeLLACS, self.QexternalIntoBottle, self.velMultTank, self.fracExpel)
            if self.adiabatic:
                runSig += 'adiabatic'
            if self.TbottleMatlConst:
                runSig += 'Tbm'+str(self.TbottleMatlConst)
            if self.timeProfileL:
                runSig += 'tPfile'+ str(self.timeProfileL) +'pcpfile'+  str(self.pcentLiqExpelledL)
            if self.Tullage:
                runSig += 'Tullage=%g'%self.Tullage

        
            objDB = diskCache.getSavedRunObj( oxName='CACHE', fuelName='RUNS', runSig=runSig,
                runType="PressInteg")
            if objDB:
                name = self.name + ' '
                diskCache.fillSelfFromDBObj( objDB, selfObj=self )
                self.name = name[:-1]
                return

        
        self.Vullage = self.ullageFrac * self.VpropTnk
        self.Vliq = self.VpropTnk - self.Vullage
        self.vdotLiq = self.Vliq / self.tAction

        #   Start out estimating the WORST CASE Volume with isentropic adiabatic
        #C  CALCULATE Isentropic GAMMA

        #C  FINAL PRESSURE
        PF = self.PfinGasOvPnom*self.PpropNom
        self.PfinalGasBot = PF
        #C  ASSUME MINIMUM TEMPERATURE FOR INITIAL PRESSURE CONDITION
      
        #c... Get starting pressure assuming minimum temperature
        self.gasObj.setTP(T=self.TmaxR, P=self.PGasTnkMEOP)
        WtMolGas = self.gasObj.WtMol

        rhoLoaded = self.gasObj.rho
        self.gasObj.setTD(T=self.TminR,D=rhoLoaded*1728.0)
        self.PbottleCold = self.gasObj.P
        self.densInit = self.gasObj.D
        gammaIsen = self.gasObj.gamma()
        
        self.compressInit = self.gasObj.compressibility()
        
        if self.PpropNom > self.PbottleCold:
        #C  THE FOLLOWING ARE ERROR CONDITION OUTPUTS
            
            print 'ERROR CONDITION IN PressurantInteg.py '
            print 'TANK PRESSURE IS HIGHER THAN Nitrogen TANK PRESSURE'
            print 'self.TmaxR, self.PGasTnkMEOP',self.TmaxR, self.PGasTnkMEOP
            print 'PbottleCold,PropNom,rhoLoaded',self.PbottleCold,self.PpropNom,rhoLoaded

            WGasTotal=(self.PbottleCold-self.PpropNom)*1000.+10000.
            WtGasResid = WGasTotal / 2.0
            WGasExpended = WGasTotal / 2.0
            wdotGas=WGasTotal/self.tAction
            TfinalPropGas = 1.0
            TfinalGasBot = 1.0
            volGasTotal=((self.PbottleCold-self.PpropNom)*100.+100.)*self.VpropTnk
            self.densInit = 0.0
            self.densFinal = 0.0
            self.densFinalProp = 0.0
        else:
            
            VpropTnk = self.Vullage + self.Vliq
            #C  make estimates of min and max bottle volume
            self.gasObj.setTP(T=self.TmaxR, P=PF)
            rhoRT = self.gasObj.rho
            
            VbotMin = max( 0.001, VpropTnk * rhoRT / (rhoLoaded-rhoRT) )
            
            VbotLLACS_min = self.WtHeLLACS / (rhoLoaded-rhoRT)

            #Tmin = self.TminR * (PF/self.PGasTnkMEOP)**( (gammaIsen-1.0)/gammaIsen )
            self.gasObj.setTD(T=self.TminR,D=rhoLoaded*1728.0)
            self.gasObj.constS_newP( PF )
            Tmin = self.gasObj.T
            
            
            self.gasObj.setTP(T=Tmin, P=PF)
            rhoCold = self.gasObj.rho
            VbotMax = VpropTnk * rhoCold / (rhoLoaded-rhoCold)
            
            VbotLLACS_max = self.WtHeLLACS / (rhoLoaded-rhoCold)
            
            print
            print 'Tmin=',Tmin,'  rhoCold=',rhoCold,'  rhoLoaded=',rhoLoaded
            print 'VbotMin=%g'%VbotMin,'   VbotMax=%g'%VbotMax,'  VpropTnk=%g'%VpropTnk

            volGasTotalEst =  (VbotMin + VbotMax + VbotLLACS_min + VbotLLACS_max) / 2.0
            print '   ...Vbottle Est=%g'%volGasTotalEst,

            # Now that 1st estimate has been made, solve for gas volume
            if type(self.Tullage) == type(1.1):
                Tullage = self.Tullage
            else:
                Tullage = self.TminR

        
            self.IntegObj = TankPress(gas=self.gas, Vbottle=volGasTotalEst, Vullage=self.Vullage,
                timeProfileL=self.timeProfileL, pcentLiqExpelledL=self.pcentLiqExpelledL,
                Vliq=self.Vliq, vdotLiq=self.vdotLiq, Pbottle=self.PbottleCold, Ptank=self.PpropNom, PVoW_Bottle=self.PVoW_Bottle, 
                PVoW_Tank=self.PVoW_Tank, Tbottle=self.TminR, Tullage=Tullage, initTullage=1,
                AccGees=self.AccGees, Nbottle=self.Nbottle, ellBottle=self.ellBottle, LcylOvDBottle=self.LcylOvDBottle, 
                Cp_effBottle=self.Cp_effBottle, Ntank=self.Ntank, ellTank=self.ellTank, LcylOvDTank=self.LcylOvDTank, 
                Cp_effTank=self.Cp_effTank, CdARegMax=self.CdARegMax, dPregulator=self.dPregulator, NtimeSteps=self.NtimeSteps, 
                heatExchangerTout=self.heatExchangerTout, WtHeLLACS=self.WtHeLLACS, QexternalIntoBottle=self.QexternalIntoBottle, 
                velMultTank=self.velMultTank, adiabatic=self.adiabatic, TbottleMatlConst=self.TbottleMatlConst)
            
            
            self.IntegObj.solveForVbottle( finalPbotOvPtank=self.PfinGasOvPnom, 
                vMax=VbotMax*2.0+VbotLLACS_max, vMin=VbotMin/2.0+VbotLLACS_min, fracExpel=self.fracExpel)
                
            self.TfinalGasBot = self.IntegObj.TbottleFinal
            self.TfinalPropGas = self.IntegObj.TpropFinal
            
            self.PfinalBlowdown = self.IntegObj.PpropFinal
            self.timeExpel = self.IntegObj.timeExpel
            
            TfinalPropGas = self.TfinalPropGas
            TfinalGasBot = self.TfinalGasBot

            volGasTotal = self.IntegObj.Vbottle
            print ' Vbottle Soln=%g'%volGasTotal,' for Pbottle =%g'%self.PGasTnkMEOP

            self.gasObj.setTP(T=self.TfinalGasBot, P=self.PfinalGasBot)
            self.densFinal = 1728.0 * self.gasObj.rho
            WtGasResid = self.gasObj.rho * self.IntegObj.Vbottle
            
            self.gasObj.setTP(T=self.TmaxR, P=self.PGasTnkMEOP)            
            WGasTotal = self.IntegObj.Vbottle * self.gasObj.rho
            WGasExpended = WGasTotal - WtGasResid

            self.gasObj.setTP(T=self.TfinalPropGas, P=self.PpropNom)
            self.densFinalProp = self.gasObj.D
        
        


        #C  AVERAGE GAS FLOW RATE
        wdotGas = WGasExpended/self.tAction
        if wdotGas < 1.0E-10:wdotGas=1.0E-10
                
        self.WtGasTotal = WGasTotal
        self.WtGasResid = WtGasResid
        self.WGasExpended = WGasExpended
        self.Vbottle = volGasTotal
        #print '   ...Vbottle =%g'%self.Vbottle,' for Pbottle =%g'%self.PGasTnkMEOP
        print '   ...TfinalPropGas=%g'%TfinalPropGas,'   TfinalGasBot=%g'%TfinalGasBot
        self.wdotGas = wdotGas
        self.TfinalPropGas=    TfinalPropGas
        self.TfinalGasBot=    TfinalGasBot
        self.gammaIsen = gammaIsen
        
        self.mass_lbm = WGasTotal
        
        self.gamPolyBottle, self.fracIsenBottle, self.gamPolyTank, self.fracIsenTank = self.calcPolytropicValues()

        if self.useDBruns:
            #  will save whether or not run was good (i.e. self.ranOK=0 or 1)
            objDB = diskCache.buildResultObjCopy( self )
            objDB = diskCache.saveResultObj( oxName='CACHE', fuelName='RUNS', runSig=runSig,
                runType="PressInteg",
                resultObj=objDB)
                


if __name__ == "__main__":  #self test
    
    
    from prism.press import PressurantHe
    # validate against Transtage
    h = PressurantInteg(name="General Test", gas='HE',
        VpropTnk=114000.0,PGasTnkMEOP=10000.0,PpropNom=80.0,
        Nbottle=4,
        PfinGasOvPnom=2.0, heatExchangerTout=None,
        tAction=15.0,TminR=510.0,TmaxR=550.0, useDBruns=0, TbottleMatlConst=None)
    print h.getMassStr()
    print
    print h.getSummary()

    hSimple = PressurantHe.PressurantHe(name="tank",  mass_lbm=0.0,
        VpropTnk=h.VpropTnk,PHeTnk=h.PGasTnkMEOP,PpropNom=h.PpropNom,
        PfinHeOvPnom=h.PfinGasOvPnom, wtHeACS=0.0,
        tAction=h.tAction,TminR=h.TminR,TmaxR=h.TmaxR,
        tPolyCorr=240.0, gamPolyCorr=1.66, gamLimPolyCorr=1.0,
        gammaPolyInp=None,
        THeTnkHX=None)
    
    hSimple.calibrate( TfBottle=h.TfinalGasBot, TfProp=h.TfinalPropGas )
    print
    print hSimple.getSummary()
    
    if 1:#h.IntegObj:
        h.IntegObj.makeExcelPlots(title='Testing 1,2,3')