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

from math import *
import sys

from prism.tanks import tank_eles
from prism.tanks import EllDome
from prism.utils.mensuration import ellipseVol
from prism.utils.Goal import Goal
from prism.MassItem import MassItem
from prism.fortran import tank_eles
from prism.props import Materials
from prism.Summary import Summary
from math import *
from prism.pov import POV_Items, POV_Basics



# comments from the ELES tank routine
#C INPUT
#C   makeCompositeTank - Make Composite Tank (0=No, 1=Yes)
#C   KALMOD-CALCULATION MODE (0=INPUT RCYLTD, 1=INPUT RINSID)
#C   SY - design strength (psi)
#C   E - material elastic modulus (psi)
#C   RHO - tank material density (lbm/in^3)
#C   TMING - minimum gauge (in)
#C   VFREE-total volume Usabel Volume (ullage plus propellant)
#          (does NOT include Residual Propellant)
#C   ELL-ellipse ratio
#C   ptank-pressure in tank
#C   SF-safety factor
#C   CXW-tank weight multiplier (welds, bosses, etc.)
#C   ITHCYL-tank wall thickness flag(0=variable thickness,1=constant)
#C   KACQUI-kind of acquisition device flag(0=none,1=transverse collapsin
#C          Alum. bladder,2=full bonded rolling diaphram(BRD) Al,3=half
#C          BRD Al,4=full BRD stainless steel,5=half BRD ss,6=capillary
#C          device, 7=elastomeric bladder, 8=Titanium Bladder)
#C   INPEX-input expulsion efficiency flag(0=calculate eff,1=input)
#C   EXPEFI-input expulsion efficiency
#C   TBLAD-bladder thickness(in)(0<KACQUI<6, or =7,8)
#C   TBOND-bond thickness(in)(1<KACQUI<6 or =7,8)
#C   TTRSPC-transverse collapsing bladder initial space distance(in)
#C          (KACQUI=1 or 7)
#C   RHOBND-density of bonding material (lbm/in**3)(1<KACQUI<6 or =7,8)
#C   RHOACQ-density of acquisition device(lbm/in**3)(KACQUI=6 or 7)

acq_dict = {0:'no acq. device',1:'transverse collapsing Alum. bladder',
    2:'full bonded rolling diaphram(BRD) Al',3:'half BRD Al',
    4:'full BRD stainless steel',5:'half BRD ss',6:'capillary device',
    7:'elastomeric transverse bladder (dP=0)',
    8:'Titanium transverse collapsing'}

ith_dict = {0:'variable thickness',1:'constant thickness'}
    
exp_dict = {0:'calculate expulsion eff',1:'input expulsion eff'}

class Tank_Conical( MassItem ):
    
    def __init__(self, name="tank",  mass_lbm=0.0, coneAngleDeg=10.0, RbigOvrRsml=1.2,
        makeCompositeTank=0,kalmod=0, matlName="Ti",  Cp_eff=0.15, # Cp Ti=.125, Al=.2, Monel=.1
        tMinGaugeUser=0.0,
        vfree=1000.0,ell=1.414,ptank=350.0,
        sf=1.5,cxw=1.25,  NumExtraBaffles=0,
        ithcyl=1,kacqui=0,inpex=0,expefi=0.99, Number=1,
        inpTblad=1, tblad=0.030,tbond=0.030,ttrspc=0.010,
        rhobnd=0.04,rhoacq=0.28,tliner=0.0,rholiner=0.1):
            
        
        MassItem.__init__(self, name, type="inert", mass_lbm=mass_lbm)
        
        self.makeCompositeTank = makeCompositeTank
        self.kalmod = kalmod
        self.matlName = matlName
        
        self.coneAngleDeg = coneAngleDeg
        self.RbigOvrRsml = RbigOvrRsml
        
        self.rho, self.sy, self.e, self.tming = Materials.getMatlProps(matlName)
        self.sy *= cos( self.coneAngleDeg * pi / 180.0 )  # account for conical section
        
        self.Cp_eff = Cp_eff  # effective Cp of all tank materials combined

        if tMinGaugeUser> self.tming:
            self.tming = tMinGaugeUser
        
        self.Number = Number  # each tank has vfree
        self.vfree = vfree
        self.ell = ell
        self.ptank = ptank
        self.sf = sf
        self.cxw = cxw
        self.ithcyl = ithcyl
        self.kacqui = kacqui
        self.NumExtraBaffles = NumExtraBaffles
        self.inpex = inpex
        self.expefi = expefi
        self.inpTblad = inpTblad
        if self.kacqui in [7,8] and not inpTblad:
            print 'WARNING... inpTblad should be equal to 1 for kacqu=7 or 8'
            self.inpTblad = 1
        
        self.tblad = tblad
        self.tbond = tbond
        self.ttrspc = ttrspc
        self.rhobnd = rhobnd
        self.rhoacq = rhoacq
        self.tliner = tliner
        self.rholiner = rholiner

        self.reCalc()
        
    def getPOV_Item(self):
        if hasattr( self, 'texture'):
            texture = self.texture
        else:
            
            texture = POV_Basics.Texture( colorName="Gray50" )
        s = POV_Items.Conical_Tank( Rbig=self.OR, Rsmall=self.Rsmall + (self.OR-self.Rbig), 
            ellRatio=self.ell, coneAngleDeg=self.coneAngleDeg,
            coneHt=self.Hcone, texture=texture)
        
        return s
        
    def setDsmallID(self, Dsmall=2.0):
        
        self.Rsmall = Dsmall / 2.0
                
        # solve for conical tank geometry so that CSE tank model can approximate it.
        def conicalTankGeom_v2( Rb ):
            RbORs = Rb / self.Rsmall
            
            # 0 deg at center, 90 deg at Rcyl  (i.e. zero degrees is zero radius)
            Eb = EllDome.EllDome( Rcyl=Rb, ellRatio=self.ell )
            Es = EllDome.EllDome( Rcyl=self.Rsmall, ellRatio=self.ell )
            
            # ellipse sections match cone at coneAngleDeg
            Vb = Eb.volume() + Eb.volumeAtPhiDeg( 90.0 - self.coneAngleDeg )[0]
            Vs = Es.volumeAtPhiDeg( 90.0 - self.coneAngleDeg )[1]
            
            # radii of ellipse sections at coneAngleDeg
            rb = Eb.getRatPhiDeg( 90.0 - self.coneAngleDeg  )
            rs = Es.getRatPhiDeg( 90.0 - self.coneAngleDeg  )
            
            # heights of ellipse sections at coneAngleDeg
            hb = rb * sin( self.coneAngleDeg * pi / 180.0 )
            hs = rs * sin( self.coneAngleDeg * pi / 180.0 )
            
            # get conical properties
            Hcone = (rb-rs)/sin( self.coneAngleDeg * pi / 180.0 )
            Vcone = pi * Hcone * (rb**2 + rb*rs + rs**2) / 3.0
            
            IHtotal = Eb.hDome + Es.hDome + hb - hs + Hcone
            Vtotal = Vcone + Vb + Vs
            #print 'for Vtotal=',Vtotal
            #print ' ==> Vcone , Vb , Vs',Vcone , Vb , Vs
            
            return RbORs,Hcone,IHtotal,Vcone,Vtotal
            
        def solveTankGeom_v2( V ):
            def getVol( Rb ):
                RbORs,Hcone,IHtotal,Vcone,Vtotal = conicalTankGeom_v2( Rb )
                #print 'At Rbig=',Rb,'Volume =',Vtotal
                return Vtotal
        
            # calc Rmax for Goal seek
            Rmax = ( V * self.ell * 3. / 4. / pi )**(1./3.)
            #print 'For Rmax=',Rmax
            
            G = Goal(goalVal=V, minX=Rmax/20., maxX=Rmax, 
                    funcOfX=getVol, tolerance=1.0E-5, maxLoops=40, failValue=Rmax)
            self.Rbig, ierror = G()
            
            self.RbigOvrRsml,self.Hcone,self.IHtotal,self.Vcone,Vtotal = conicalTankGeom_v2( self.Rbig )
            
            #print 'Rsmall,Hcone,IHtotal,Vtotal=',self.Rsmall,self.Hcone,self.IHtotal,Vtotal
        
        volMin = ellipseVol(Dsmall, self.ell) * 2.0 # 2 domes = full tank volume
        if self.vfree > volMin:
            solveTankGeom_v2( self.vfree )
        else:
            self.RbigOvrRsml = 1.000001 # let tank be oversized
        self.reCalc()
        
    def setToOH(self, Hgoal=25.0):
        
        self.RbigOvrRsml = 1.01
        self.reCalc()
        
        print 'self.OH=%s, Hgoal=%s'%(self.OH, Hgoal)
        if self.OH < Hgoal:
            rorMin = 1.01
            rorMax = 3.
        
            def testHt( ror ):
                self.RbigOvrRsml = ror
                self.reCalc()
                print 'ror=%s, OH=%s'%(ror, self.OH)
                return self.OH
            
            G = Goal(goalVal=Hgoal, minX=rorMin, maxX=rorMax, 
                    funcOfX=testHt, tolerance=1.0E-5, maxLoops=40, failValue=self.vfree)
            ror, ierror = G()
            self.RbigOvrRsml = ror
            self.reCalc()
        
    def reCalc(self):
        
        default_cxw = 1.0
        kacqui = self.kacqui
        if self.kacqui in [7]:
            kacqui = 1
        if self.kacqui in [8]:
            kacqui = 2 # like Al BRD
        
        tliner_eff = self.tliner 
        rholiner_eff = self.rholiner
        
        # solve for conical tank geometry so that CSE tank model can approximate it.
        def conicalTankGeom_v1( Rb ):
            Rs = Rb / self.RbigOvrRsml
            
            # 0 deg at center, 90 deg at Rcyl  (i.e. zero degrees is zero radius)
            Eb = EllDome.EllDome( Rcyl=Rb, ellRatio=self.ell )
            Es = EllDome.EllDome( Rcyl=Rs, ellRatio=self.ell )
            
            # ellipse sections match cone at coneAngleDeg
            Vb = Eb.volume() + Eb.volumeAtPhiDeg( 90.0 - self.coneAngleDeg )[0]
            Vs = Es.volumeAtPhiDeg( 90.0 - self.coneAngleDeg )[1]
            
            # radii of ellipse sections at coneAngleDeg
            rb = Eb.getRatPhiDeg( 90.0 - self.coneAngleDeg  )
            rs = Es.getRatPhiDeg( 90.0 - self.coneAngleDeg  )
            
            # heights of ellipse sections at coneAngleDeg
            hb = rb * sin( self.coneAngleDeg * pi / 180.0 )
            hs = rs * sin( self.coneAngleDeg * pi / 180.0 )
            
            # get conical properties
            Hcone = (rb-rs)/sin( self.coneAngleDeg * pi / 180.0 )
            Vcone = pi * Hcone * (rb**2 + rb*rs + rs**2) / 3.0
            
            IHtotal = Eb.hDome + Es.hDome + hb - hs + Hcone
            Vtotal = Vcone + Vb + Vs
            #print 'for Vtotal=',Vtotal
            #print ' ==> Vcone , Vb , Vs',Vcone , Vb , Vs
            
            return Rs,Hcone,IHtotal,Vcone,Vtotal
            
        def solveTankGeom_v1( V ):
            def getVol( Rb ):
                Rs,Hcone,IHtotal,Vcone,Vtotal = conicalTankGeom_v1( Rb )
                #print 'At Rbig=',Rb,'Volume =',Vtotal
                return Vtotal
        
            # calc Rmax for Goal seek
            Rmax = ( V * self.ell * 3. / 4. / pi )**(1./3.)
            #print 'For Rmax=',Rmax
            
            G = Goal(goalVal=V, minX=Rmax/20., maxX=Rmax, 
                    funcOfX=getVol, tolerance=1.0E-5, maxLoops=40, failValue=Rmax)
            self.Rbig, ierror = G()
            
            self.Rsmall,self.Hcone,self.IHtotal,self.Vcone,Vtotal = conicalTankGeom_v1( self.Rbig )
            
            #print 'Rsmall,Hcone,IHtotal,Vtotal=',self.Rsmall,self.Hcone,self.IHtotal,Vtotal
        
        solveTankGeom_v1( self.vfree )
        
        LcylOvrDiamEst = self.Hcone / (self.Rsmall + self.Rbig)

        self.rinsid,hinsid,self.cyl,self.wacqui,\
            self.vacqui,self.dpacq,self.pullag,self.vresid,\
            self.vtank,self.thkcyl,self.thkend,self.thkBladOut,self.wliner,self.wtank= \
            tank_eles.csetnk(self.makeCompositeTank,self.kalmod,
            self.sy,self.e,self.rho,self.tming,self.vfree,
            self.ell,LcylOvrDiamEst,self.ptank,self.sf,default_cxw,
            self.ithcyl,kacqui,self.inpex,self.expefi,
            self.inpTblad,self.tblad,self.tbond,self.ttrspc,
            self.rhobnd,self.rhoacq,tliner_eff,rholiner_eff)

        if self.kacqui==7: # user defined bladder (elastomeric)
            self.dpacq = 0.0
            self.wacqui = self.wacqui * self.rhoacq * self.tblad / self.thkBladOut / 0.1
            self.vacqui = self.vacqui * self.tblad / self.thkBladOut
            self.thkBladOut = self.tblad

        if self.kacqui==8: # Titanium Bladder
            self.dpacq = 2.0 * self.dpacq # estimate Ti is 4x Al
            self.wacqui = self.wacqui * 0.16 * self.tblad / self.thkBladOut / 0.1
            self.vacqui = self.vacqui * self.tblad / self.thkBladOut
            #self.thkBladOut = self.tblad

        self.dinsid = 2.0*self.rinsid
        self.pov_h = self.cyl + self.dinsid/self.ell
        self.pov_w = self.dinsid
        self.pov_d = self.dinsid

        self.Across = pi * self.rinsid**2

        ECC=sqrt(1.-(1./self.ell)**2)
        HEADCO=pi/2.
        if self.ell>1.0:
            HEADCO=(pi/4.)*(1.+log((1.+ECC)/(1.-ECC))/(2.*ECC*self.ell**2))
        HSA = HEADCO*self.dinsid**2
        self.SAinsid = pi*self.dinsid*self.cyl + 2.*HSA
        
        if self.NumExtraBaffles > 0:
            self.WtExtraBaffles = pi * self.rinsid**2 * self.tming * self.rho * self.NumExtraBaffles
        else:
            self.WtExtraBaffles = 0.0
        
        self.mass_lbm = (self.wtank + self.wacqui ) * self.cxw   + self.WtExtraBaffles 
        
        if self.thkcyl <= self.tming:
            self.Pburst = self.tming * self.sy / self.rinsid 
        else:
            self.Pburst = self.sf * self.ptank
            
        self.OD = 2.0* (self.Rbig + self.thkcyl + self.thkBladOut + self.tliner)
        self.OR = self.OD / 2.0
        
        self.OH = self.IHtotal +  2.0* ( self.thkend + self.thkBladOut + self.tliner)

        # if more than one tank, increase mass to reflect that
        self.mass_lbm *= self.Number

        
        try:
            self.PVoverW = self.ptank*self.vtank*self.Number/self.mass_lbm
            self.PburstVoverW = self.Pburst*self.vtank*self.Number/self.mass_lbm
        except:
            self.PVoverW = 0.0
            self.PburstVoverW = 0.0
    
    def minGaugeStr(self, t):
        if t<= self.tming:
            return 'Min Gauge'
        else:
            return ''
        
    def buildSummary(self):
        
        summ = Summary(  summName='Conical Tank with Elliptical Domes',
        componentName=self.name, mass_lbm=self.mass_lbm, type=self.type)
        
        
        if self.makeCompositeTank:
            summ.addAssumption('Composite Tank Algorithm')
        else:
            summ.addAssumption('Metal Tank Algorithm')
            summ.addAssumption('Structural Material : ' + self.matlName )



        if self.NumExtraBaffles > 0:
            summ.addAssumption("%i Extra Baffles are added to tank"%self.NumExtraBaffles)
            
        if self.kacqui in [1,2,3,4,5,7,8]:
            if self.inpTblad:
                summ.addAssumption('Bladder Thickness is input at %.3f in'%self.tblad)
            else:
                summ.addAssumption('Bladder Thickness is calculated at %.3f in'%self.tblad)
            
        summ.addAssumption('kalmod = ' + str( self.kalmod ) )
            
        if self.Number>1:
            summ.addAssumption( 'Mass is for %i Tanks total'%self.Number )
        
        
        # add inputs
        summ.addInput('vfree', self.vfree, 'cuin', '%g')
        summ.addInput('vfreeTotal', self.vfree*self.Number, 'cuin', '%g')
        
        
        summ.addInput('ell', self.ell, '', '%g')
        summ.addInput('RbigOvrRsml', self.RbigOvrRsml, '', '%g')
        summ.addInput('coneAngleDeg', self.coneAngleDeg, 'deg', '%.2f')
        summ.addInput('ptank', self.ptank, 'psia', '%g')
        summ.addInput('sf', self.sf, '', '%g')
        summ.addInput('cxw', self.cxw, '', '%g')
        summ.addInput('ithcyl', self.ithcyl, '', '%g')
        summ.addInput('kacqui', self.kacqui, acq_dict[self.kacqui], '%s')
        summ.addInput('inpex', self.inpex, '', '%g')
        summ.addInput('expefi', self.expefi, '', '%g')
        summ.addInput('tblad', self.tblad, 'in', '%.3f')
        summ.addInput('tbond', self.tbond, 'in', '%.3f')
        summ.addInput('ttrspc', self.ttrspc, 'in', '%.3f')
        summ.addInput('rhobnd', self.rhobnd, 'lbm/cuin', '%g')
        summ.addInput('rhoacq', self.rhoacq, 'lbm/cuin', '%g')
        summ.addInput('tliner', self.tliner, 'in', '%.3f')
        summ.addInput('rholiner', self.rholiner, 'lbm/cuin', '%g')
        
        # add outputs

        summ.addOutput( 'IRbig', self.Rbig, 'in', '%.3f' )
        summ.addOutput( 'IDbig', self.Rbig*2., 'in', '%.3f' )
        summ.addOutput( 'IRsmall', self.Rsmall, 'in', '%.3f' )
        summ.addOutput( 'IDsmall', self.Rsmall*2., 'in', '%.3f' )
        summ.addOutput( 'OR', self.OR, 'in', '%.3f' )
        summ.addOutput( 'OD', self.OD, 'in', '%.3f' )
        summ.addOutput( 'OH', self.OH, 'in', '%.3f' )
        summ.addOutput( 'IHtotal', self.IHtotal, 'in', '%.3f' )
        summ.addOutput( 'SAinsid', self.SAinsid, 'sqin', '%.3f' )
        summ.addOutput( 'Hcone', self.Hcone, 'in', '%.3f' )
        summ.addOutput( 'wacqui', self.wacqui, 'lbm', '%.3f' )
        summ.addOutput( 'vacqui', self.vacqui, 'cuin', '%g' )
        summ.addOutput( 'dpacq', self.dpacq, 'psig', '%g' )
        summ.addOutput( 'pullag', self.pullag, 'psia', '%g' )
        summ.addOutput( 'vresid', self.vresid, 'cuin', '%g' )
        summ.addOutput( 'Vcone', self.Vcone, 'cuin', '%g' )
        summ.addOutput( 'vtank', self.vtank, 'cuin', '%g' )

        summ.addOutput( 'tming', self.tming, 'in', '%.3f' )
        if self.thkcyl <= self.tming:
            summ.addOutput( 'thkcyl', self.thkcyl, 'in (use tming)', '%.3f' )
        else:
            summ.addOutput( 'thkcyl', self.thkcyl, 'in', '%.3f' )
            
        if self.thkend <= self.tming:
            summ.addOutput( 'thkend', self.thkend, 'in (use tming)', '%.3f' )
        else:
            summ.addOutput( 'thkend', self.thkend, 'in', '%.3f' )
        
        summ.addOutput( 'thkBladOut', self.thkBladOut, 'in', '%.3f' )
        if self.NumExtraBaffles > 0:
            summ.addOutput( 'WtExtraBaffles', self.WtExtraBaffles, 'lbm', '%.3f' )


        if self.wliner > 0.0:
            summ.addOutput( 'wliner', self.wliner, 'lbm', '%.3f' )
            summ.addOutput( 'wtank(+liner)', self.wtank, 'lbm', '%.3f' )
        else:
            summ.addOutput( 'wtank', self.wtank, 'lbm', '%.3f' )
            
        summ.addOutput( 'rho', self.rho, 'lbm/cuin', '%g' )
        if not self.makeCompositeTank:
            summ.addOutput( 'sy', self.sy, 'psi', '%g' )
            summ.addOutput( 'e', self.e, 'psi', '%g' )
        summ.addOutput( 'PmeopVoverW', self.PVoverW, 'lbf-in/lbm', '%g' )
        summ.addOutput( 'Pburst(est.)', self.Pburst, 'psia', '%.1f' )
        summ.addOutput( 'PburstVoverW', self.PburstVoverW, 'lbf-in/lbm', '%g' )

        return summ


if __name__ == "__main__":  #self test

    print "Actual FFC Prop Tank =", 3.703,"lbm"
    print "Calculated =   ",
    oxFFC = Tank_Conical(name="FFC Propellant Tank", mass_lbm=3.703,
        makeCompositeTank=1, kalmod=0,  
        matlName="grEpox",   vfree=486.0,ell=1.767,
        ptank=1400.0,sf=1.5,cxw=1.5,
        ithcyl=1,kacqui=1,inpex=1,expefi=0.98,
        tblad=0.030,tbond=0.030,ttrspc=0.010,
        rhobnd=0.04,rhoacq=0.098,tliner=0.03,rholiner=0.098)
    print oxFFC.getMassStr()
    print
    print oxFFC.getSummary()

    print '===================================================================='
    print "Actual Mars PMD Prop Tank =", 8.5,"lbm (PSI SN 80353-1)"
    print "Calculated Mars =   ",
    pmd = Tank_Conical(name="Mars Observer Propellant Tank", mass_lbm=8.5,
        makeCompositeTank=0, kalmod=0,  
        matlName="Ti",   vfree=3490.0,ell=1.0,
        ptank=400.0,sf=1.5,cxw=1.07, NumExtraBaffles=3,
        ithcyl=1,kacqui=6,inpex=1,expefi=0.98,
        tblad=0.0,tbond=0.0,ttrspc=0.0,rhoacq=0.16)
    print pmd.getMassStr()
    print
    print pmd.getSummary()
    print '===================================================================='
    print "Spherical  Prop Tank =", 12.04,"lbm w/o diaphram or PK baffles (PSI SN 80274-1)"
    print "Calculated Mass =   ",
    pmd = Tank_Conical(name="PSI 80274-1 Tank", mass_lbm=12.04, tMinGaugeUser=0.035,
        makeCompositeTank=0, kalmod=0,  
        matlName="Ti",   vfree=3660.0,ell=1.0,
        ptank=377.0,sf=1.5,cxw=1.32, NumExtraBaffles=3,
        ithcyl=1,kacqui=6,inpex=1,expefi=0.98,
        tblad=0.0,tbond=0.0,ttrspc=0.0,rhoacq=0.16)
    print pmd.getMassStr()
    print
    print pmd.getSummary()
    
    print '===================================================================='
    print "Actual FFC He Tank=",3.157,"lbm (single tank)"
    print "Calculated =",
    HeFFC = Tank_Conical(name="FFC Helium Tank", mass_lbm=3.157,
        makeCompositeTank=1, kalmod=0,  Number=1,
        matlName="grEpox",   vfree=160.0, ell=1.764,
        ptank=10000.0,sf=2.0,cxw=1.25,
        ithcyl=1,kacqui=0,inpex=0,
        inpTblad=0, tblad=0.030,tbond=0.030,ttrspc=0.010,
        rhobnd=0.04,rhoacq=0.1,tliner=0.04,rholiner=0.1)
    print HeFFC.getMassStr()
    print
    print HeFFC.getSummary()
    #sys.exit()

    #HeFFC.setToMaxID( IDmax=4.0)
    #print
    #print HeFFC.getSummary()
    HeFFC.setToOH( Hgoal=6.0 )
    print
    print HeFFC.getSummary()
    
    print '+'*66
    print 'Testing setDsmallID'
    HeFFC.setDsmallID( Dsmall=5.5)
    print HeFFC.getSummary()
    
    if 0:
        HeFFC.setToLcyl( 7.0 )
        #HeFFC.setToLength( L=10.0)
        print
        print HeFFC.getSummary()