def load_elements(self):

        #Set atomic data 
        #atomicData.setDataFile('he_i_rec_Pal12-Pal13.fits')
        atomicData.setDataFile('s_iii_coll_HRS12.dat')
 
        #Default: 's_iii_atom_PKW09.dat' 
        'S3: All energy and A values: Podobedova, Kelleher, and Wiese 2009, J. Phys. Chem. Ref. Data, Vol.'
        'S3: collision strengths: Tayal & Gupta 1999, ApJ, 526, 544'

        #New Atomic data s_iii_coll_HRS12.dat
        'S3: All energy and A values: Podobedova, Kelleher, and Wiese 2009, J. Phys. Chem. Ref. Data, Vol.'
        'S3: collision strengths: Hudson, Ramsbottom & Scott 2012, ApJ, 750, 65'
        
        #Declare ions
        self.S2_atom            = Atom('S', 2)
        self.S3_atom            = Atom('S', 3)
        self.Ar3_atom           = Atom('Ar', 3)
        self.Ar4_atom           = Atom('Ar', 4)
        self.N2_atom            = Atom('N', 2)
        self.O2_atom            = Atom('O', 2)
        self.O3_atom            = Atom('O', 3)
        self.H1_atom            = RecAtom('H', 1)
        self.He1_atom           = RecAtom('He', 1)
        self.He2_atom           = RecAtom('He', 2)
        
        #Pyneb objects
        self.diags              = Diagnostics()
        
        #Ohrs 2016 relation for the OI_SI gradient
        self.logSI_OI_Gradient  = random.normal(-1.53,  0.05, size = self.MC_array_len) # random.normal(-1.78,  0.03, size = self.MC_array_len)
        self.OI_SI              = power(10, -self.logSI_OI_Gradient)
        
        #Theoretical ratios
        self.S3_ratio           = self.S3_atom.getEmissivity(10000, 100, wave = 9531) / self.S3_atom.getEmissivity(10000, 100, wave = 9069)
        self.S3_9000_ratio      = random.normal(self.S3_atom.getEmissivity(10000, 100, wave = 9531) / self.S3_atom.getEmissivity(10000, 100, wave = 9069),  0.01, size = self.MC_array_len)
        self.N2_6000_ratio      = self.N2_atom.getEmissivity(10000, 100, wave = 6584) / self.N2_atom.getEmissivity(10000, 100, wave = 6548)
        self.O3_5000_ratio      = self.O3_atom.getEmissivity(10000, 100, wave = 5007) / self.O3_atom.getEmissivity(10000, 100, wave = 4959)
        
        #Factors to speed calculations
        self.lines_factors = {}
        self.lines_factors['S3_9069A'] = 1 + self.S3_ratio
        self.lines_factors['S3_9531A'] = 1 + 1/self.S3_ratio
        
        #Cloudy models for the SIV contribution
        
        self.m_SIV_correction   = random.normal(1.1628,  0.00559, size = self.MC_array_len)
        self.n_SIV_correction   = random.normal(0.0470,  0.0097, size = self.MC_array_len)
        #self.m_SIV_correction   = random.normal(1.109,  0.01, size = self.MC_array_len)
        #self.n_SIV_correction   = random.normal(0.135,  0.0173, size = self.MC_array_len)
 
        #CHAOS relation TNII-TSIII
        #T[SIII]  = 1.312(+-0.075)T[NII]-0.313(+-0.058)
                #TNII     = (0.762+-0.044)*TSIII  + 0.239+-0.046
        self.m_TNII_correction   = random.normal(0.762,  0.044, size = self.MC_array_len)
        self.n_TNII_correction   = random.normal(0.239,  0.046, size = self.MC_array_len)       
        
        
        #Truncated gaussian for the density
        lower_trunc, upper_trunc = (1.0 - 50.0) / 25.0, (100 - 50) / 25.0
        self.Truncated_gaussian  = truncnorm(lower_trunc,  upper_trunc, loc = 50, scale = 25)
       
        print '-Elements loaded\n'
       
        return
S3: All energy and A values: Podobedova, Kelleher, and Wiese 2009, J. Phys. Chem. Ref. Data, Vol.