Python metallicity.printsafemulti函数代码示例

 def calcNIISII(self):
     if self.hasS2 and self.hasN2:
         self.N2S2=self.N26584/self.S26717+self.dustcorrect(k_N2,k_S2,flux=True)#0.4*self.mds['E(B-V)']*(k_N2-k_S2) 
         self.logN2S2=np.log10(self.N26584/self.S26717)+self.dustcorrect(k_N2,k_S2)#0.4*self.mds['E(B-V)']*(k_N2-k_S2) 
         self.hasN2S2=True
     else: 
         printsafemulti( "WARNING: needs SII6717 and NII6584 to calculate NIISII: did you run setN2() and setS?",self.logf,self.nps)

    def setSII(self,S26717,S26731,S39069,S39532):
        if S26717 is not None and sum(S26717>0)>0:
            self.S26717=S26717
            self.hasS2=True

            if self.hasHa:
                self.logS2Ha=np.log10(self.S26717/self.Ha)+self.dustcorrect(k_S2,k_Ha)               
            else: 
                printsafemulti( "WARNING: needs SII6717 and Ha to calculate SIIHa: did you run setHab() and setS()?",self.logf,self.nps)
        if S26731 is not None and sum(S26731>1e-9)>0:
            self.S26731=S26731
            self.hasS26731=True
        if S39069 is not None and sum(S39069>1e-9)>0:
            self.S39069=S39069
            self.hasS39069=True
        if S39532 is not None and sum(S39532>1e-9)>0:
            self.S39532=S39532
            self.hasS39532=True
        if self.hasS2 :
            
            if self.hasN2 and self.NII_SII is None and self.hasS26731:
                self.NII_SII=np.log10(self.N26584/(self.S26717+self.S26731))#+self.dustcorrect(k_N2,k_O2,flux=True) 
                #lines are very close: no dust correction            
            if self.hasO3  and self.OIII_SII is None and self.hasS26731:
                self.OIII_SII=np.log10(self.O35007/(self.S26717+self.S26731)+self.dustcorrect(k_O3,k_S2,flux=True) )

 def calclogq(self,Z):
     if not self.hasO3O2:
         printsafemulti(  "WARNING: needs O3,O2,Hb to calculate logq properly.",self.logf,self.nps)
         return -1
     if self.logO3O2sq is None:
         self.logO3O2sq=self.logO3O2**2
     return (32.81 -1.153*self.logO3O2sq + Z*(-3.396 -0.025*self.logO3O2 + 0.1444*self.logO3O2sq))/(4.603-0.3119*self.logO3O2-\
             0.163*self.logO3O2sq+ Z*(-0.48 + 0.0271*self.logO3O2+ 0.02037*self.logO3O2sq)) 

 def calcS23(self):
     printsafemulti(  "calculating S23",self.logf,self.nps)
     #the original code here uses S267176731, 
     #which is however set to 6717 as default
     #Vilchez & Esteban (1996)
     if  self.hasS2 :
         if self.hasS39069 and self.hasHb:
             self.logS23=np.log10((self.S26717/self.Hb)*self.dustcorrect(k_S2,k_Hb,flux=True) + (self.S39069/self.Hb)*self.dustcorrect(k_S3,k_Hb,flux=True))                                 

 def calcD02(self):
     # [NII]/Ha Denicolo, Terlevich & Terlevich (2002), MNRAS, 330, 69
     #FED:added uncertainties
     printsafemulti(  "calculating D02",self.logf,self.nps)
     
     e1=np.random.normal(0,0.05,self.nm)
     e2=np.random.normal(0,0.1,self.nm)
     if self.hasN2 and self.hasHa:
         self.mds['D02'] = 9.12+e1+(0.73+e2)*self.logN2Ha
     else:
         printsafemulti(  "WARNING: need N2Ha to do this. did you run setHab and setNII",self.logf,self.nps)

 def calcP(self):
     if self.P is None:
         if self.logR23 is None:
             printsafemulti(  "WARNING: Must first calculate R23",self.logf,self.nps)
             self.calcR23()
             if self.logR23 is None:
                 printsafemulti(  "WARNING: Cannot compute this without R23",self.logf,self.nps)
                 return -1
         #R3=10**self.logO349595007Hb
         #R2=10**self.logO2Hb
         #P = R3/(R2+R3)
         self.P=self.R3/self.R23

    def calcC01_ZR23(self):
        # C01 = Charlot, S., & Longhetti, M., 2001, MNRAS, 323, 887
        # Charlot 01 R23 calibration: (case F) ##        
        # available but deprecated
        printsafemulti(  "calculating C01",self.logf,self.nps)

        if self.hasO3 and self.hasO2 and self.hasO3Hb:
            x2=self.O2O35007/1.5
            x3=(10**self.logO3Hb)*0.5
            self.mds['C01_R23']=np.zeros(self.nm)+float('NaN')        
            self.mds['C01_R23'][self.O2O35007<0.8]=np.log10(3.78e-4 * (x2[self.O2O35007<0.8])**0.17 * x3[self.O2O35007<0.8]**(-0.44))+12.0    
         
            self.mds['C01_R23'][ self.O2O35007 >= 0.8]=np.log10(3.96e-4 * x3[self.O2O35007 >= 0.8]**(-0.46))+12.0   
        else:
            printsafemulti('''WARNING: need [OIII]5700, [OII]3727, and Ha to calculate calcC01_ZR23, 
did you set them up with  setOlines()?'''        ,self.logf,self.nps)

        # Charlot 01 calibration: (case A) based on [N2]/[SII]##
        # available but deprecated
        if not self.hasN2S2:
            printsafemulti(  "WARNING: trying to calculate logNIISII",self.logf,self.nps)
            self.calcNIISII()
        if self.hasN2S2 and self.hasO3 and self.hasO2 and self.hasO3Hb:
            self.mds['C01_N2S2']=np.log10(5.09e-4*(x2**0.17)*((self.N2S2/0.85)**1.17))+12
        else:
            printsafemulti('''WARNING: needs [NII]6584, [SII]6717, [OIII]5700, [OII]3727, and Ha to calculate calcC01_ZR23, 
did you set them up with  setOlines() and ?''',self.logf,self.nps)        

    def calcKD02_N2O2(self):
        ##  Kewley & Dopita (2002) estimates of abundance 
        ##  KD02
        # KD02 [N2]/[O2] estimate (can be used for whole log(O/H)+12 range, 
        # but rms scatter increases to 0.11 rms for log(O/H)+12 < 8.6
        # rms = 0.04 for
        # log(O/H)+12 > 8.6
        # uses equation (4) from KD02 paper
        # FED: i vectorized the hell out of this function!!! 
        # from a 7 dimensional if/for loop to 1 if and 1 for :-D
        #vectorizing makes fed happy ...

        printsafemulti(  "calculating KD02_N2O2",self.logf,self.nps)

        if self.hasN2 and self.hasO2 and self.hasHa and self.hasHb:         
            self.mds['KD02_N2O2']=np.zeros(self.nm)+float('NaN')
            if not self.hasN2O2:
                printsafemulti(  "WARNING: must calculate logN2O2 first",self.logf,self.nps)
                self.calcNIIOII()
            if  not self.hasN2O2 or self.N2O2_roots  is None or sum(np.isnan(self.N2O2_roots.flatten())) == len(self.N2O2_roots.flatten()):
                printsafemulti(  "WARNING:  cannot calculate N2O2",self.logf,self.nps)
                return -1
            roots=self.N2O2_roots.T
            for k in range(4):
                indx=(abs(roots[k]) >= 7.5) * (abs(roots[k]) <= 9.4) * (roots[k][:].imag ==  0.0 )
                self.mds['KD02_N2O2'][indx]=abs(roots[k][indx]) 
        else:
            printsafemulti(  "WARNING: need NII6584 and OII3727 and Ha and Hb to calculate this. did you run setO() setHab() and setNII()?",self.logf,self.nps)
        return 1

 def calcDP00(self):
     # Diaz, A. I., & Perez-Montero, E. 2000, MNRAS, 312, 130 
     # As per KD02: DP00 diagnostic systematically underestimates the
     # abundance relative to the comparison abundance. 
     # A term is added to improve the fit according to KD02 Eq. 6
     # AVAILABLE BUT DEPRECATED
     printsafemulti(  "calculating DP00",self.logf,self.nps)
     
     if self.logS23  is None:
         self.calcS23()
         if self.logS23 is None:
             printsafemulti(  "WARNING: Cannot compute this without S23",self.logf,self.nps)
             return -1
     self.mds['DP00'] = 1.53*self.logS23+8.27+1.0/(2.0-9.0*self.logS23**3)

    def calcR23(self):
        printsafemulti( "calculating R23",self.logf,self.nps)

        #R23 NEW Comb, [NII]/Ha: KK04 = Kobulnicky & Kewley, 2004, submitted'
        if  self.hasO3   and self.hasO2 and self.hasHb:
            self.R2=(self.O23727/self.Hb)*self.dustcorrect(k_O2,k_Hb, flux=True) 
            self.R3=(self.O34959p5007/self.Hb)*self.dustcorrect(k_O3,k_Hb, flux=True) 
            self.R23=self.R2+self.R3
            self.logR23=np.log10(self.R23)
            self.mds['logR23']=self.logR23
            #note that values of logR23 > 0.95 are unphysical. 
            #you may choose to uncomment the line below
            #self.logR23[self.logR23>0.95]=0.95
        else:
            printsafemulti( "WARNING: need O3, O2, Hb",self.logf,self.nps)

 def setNII(self,N26584):
     if N26584 is not None and sum(N26584>0):
         self.N26584=N26584
         self.hasN2=True
         if self.hasHa :
             self.logN2Ha=np.log10(self.N26584/self.Ha)#+self.dustcorrect(k_N2,k_Ha,flux=True) 
             #lines are very close: no dust correction
             #Note: no dust correction cause the lies are really close!
         else: 
             printsafemulti( "WARNING: needs NII6584 and Ha to calculate NIIHa: did you run setHab()?",self.logf,self.nps)
         if self.hasS2 and self.hasS26731 and self.hasN2:
             self.NII_SII=np.log10(self.N26584/(self.S26717+self.S26731))#+self.dustcorrect(k_N2,k_S2,flux=True) 
             #lines are very close: no dust correction
         if self.hasO2 and self.hasN2:
             self.NII_OII=np.log10(self.N26584/self.O23727+self.dustcorrect(k_N2,k_O2,flux=True) )

    def fz_roots(self,coef): 
        if len(coef.shape)==1:
            coef[~(np.isfinite(coef))]=0.0
            rts= np.roots(coef[::-1])
            if rts.size==0:
                printsafemulti( 'WARNING: fz_roots failed',self.logf,self.nps)
                rts=np.zeros(coef.size-1)
            return rts

        else:
            rts=np.zeros((coef.shape[0],coef.shape[1]-1),dtype=complex)
            coef[~(np.isfinite(coef))]=0.0

            for i in range(coef.shape[0]):
                rts[i]= np.roots(coef[i][::-1])#::-1][0])
            return rts

    def calcM13(self):
        #Marino+ 2013
        printsafemulti(  "calculating M13",self.logf,self.nps)

        if not self.hasHa  or not self.hasN2:
            printsafemulti(  "WARNING: need O3, N2, Ha and Hb, or at least N2 and Ha",self.logf,self.nps)
            return -1
        else:
            e1=np.random.normal(0,0.027,self.nm)
            e2=np.random.normal(0,0.024,self.nm)
            self.mds["M13_N2"] = 8.743+e1 - (0.462+e2)*self.logN2Ha
            if   self.hasHb and self.hasO3:
                e1=np.random.normal(0,0.012,self.nm)
                e2=np.random.normal(0,0.012,self.nm)
                O3N2=self.logO3Hb-self.logN2Ha
                self.mds["M13_O3N2"] = 8.533+e1 - (0.214+e1)*O3N2

    def calcZ94(self):
        ### calculating z from Kobulnicky,Kennicutt,Pizagno (1998)
        ### parameterization of Zaritzky et al. (1994)
        ###Z94 = Zaritsky, D., Kennicutt, R. C., & Huchra, J. P. 1994,
        ###ApJ, 420, 87
        ### only valid on the upper branch of R23 (KE08 A2.4)

        #printsafemulti(  "calculating Z94",self.logf,self.nps)
        if self.logR23 is None:
            printsafemulti(  "WARNING: Must first calculate R23",self.logf,self.nps)
            self.calcR23()
            if self.logR23 is None:
                printsafemulti(  "WARNING: Cannot compute this without R23",self.logf,self.nps)
                return -1
        self.mds['Z94']=nppoly.polyval(self.logR23, [9.265,-0.33,-0.202,-0.207,-0.333])
        self.mds['Z94'][(self.logR23 > 0.9)]=None

    def calcKK04_R23(self):
        # Kobulnicky & Kewley 2004
        # calculating upper and lower metallicities for objects without
        # Hb  and for objects without O3 and/or O2

        printsafemulti(  "calculating KK04_R23",self.logf,self.nps)

        #this is in the original code but not used :( 
        #if self.hasN2 and self.hasHa:
        #logq_lims=[6.9,8.38]
        #logN2Ha=np.log10(self.N26584/self.Ha) CHECK!! why remove dust correction??
        #Z_new_N2Ha_lims= np.atleast_2d([1.0,1.0]).T*nppoly.polyval(self.logN2Ha,[7.04, 5.28,6.28,2.37])-
        #np.atleast_2d( logq_lims).T*nppoly.polyval(self.logN2Ha,[-2.44,-2.01,-0.325,0.128])+
        #np.atleast_2d(logq_lims).T*(10**(self.logN2Ha-0.2)*(-3.16+4.65*self.logN2Ha)) 
        # R23 diagnostics from Kobulnicky & Kewley 2004

        Zmax=np.zeros(self.nm)
        # ionization parameter form logR23
        if not self.hasO3O2:
            logq=np.zeros(self.nm)
        else:
            if self.Z_init_guess is None:
                self.initialguess()
            Z_new=self.Z_init_guess.copy()
            if self.logR23 is None:
                printsafemulti(  "WARNING: Must first calculate R23",self.logf,self.nps)
                self.calcR23()
            if self.logR23 is None:
                printsafemulti(  "WARNING: Cannot compute this without R23" ,self.logf,self.nps)
            else:
                logqold,convergence,ii=np.zeros(self.nm)+100,100,0
                tol=1e-4
                #3 iterations are typically enought to achieve convergence KE08 A2.3
                while convergence>tol and ii<100:
                    Zmax=Zmax*0.0
                    ii+=1
                    logq=self.calclogq(Z_new)
                    Zmax[(logq >= 6.7) * (logq < 8.3)]=8.4
                    # maximum of R23 curve:               
                    Z_new=nppoly.polyval(self.logR23,[9.72, -0.777,-0.951,-0.072,-0.811])-\
                           logq*nppoly.polyval(self.logR23,[0.0737,  -0.0713, -0.141, 0.0373, -0.058])
                    indx=self.Z_init_guess<=Zmax
                    Z_new[indx]=nppoly.polyval(self.logR23[indx], [9.40 ,4.65,-3.17])-\
                                 logq[indx]*nppoly.polyval(self.logR23[indx],[0.272,0.547,-0.513])
                    convergence=np.abs((logqold-logq).mean())
                    logqold=logq.copy()
                if ii>=100:
                    printsafemulti(  "WARNING: loop did not converge" ,self.logf,self.nps)
                    Z_new=np.zeros(self.nm)+float('NaN')
                Z_new_lims=[nppoly.polyval(self.logR23,[9.40, 4.65,-3.17])-\
                            logq*nppoly.polyval(self.logR23,[0.272,0.547,-0.513]),
                            nppoly.polyval(self.logR23,[9.72, -0.777,-0.951,-0.072,-0.811])-\
                            logq*nppoly.polyval(self.logR23,[0.0737,  -0.0713, -0.141, 0.0373, -0.058])]
                Z_new[(Z_new_lims[0]>Z_new_lims[1])]=None
                self.mds['KK04_R23']=Z_new