DA_Pollution_Tables/InterpDiffusion.py at master · evbauer/DA_Pollution_Tables · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

import numpy as np
import cgs_const as cgs
from scipy import interpolate

# Build interpolators for cvz mass and diffusion timescales
# for Ca based on the master tables.


# Note: column 0 for Teff in these files just represents the label of the run,
# column 1 is the actual model Teff
# at the end of the run that started at the temperature of column 0.

master = np.genfromtxt('diffusion_timescales/WD_M_0.9.data')
fm85 = interpolate.interp1d(master[:,1],master[:,3],fill_value='extrapolate') # log(M
def Mcvz85(Teff):
    return (10**fm85(Teff))*0.9*cgs.msun # in grams
ft85 = interpolate.interp1d(master[:,1],master[:,10],fill_value='extrapolate') # log(taudiff/yr) as a function of Teff
def tau85(Teff):
    return (10**ft85(Teff))*cgs.year # in years
logg85 = interpolate.interp1d(master[:,1],master[:,2],fill_value='extrapolate')
# needed because loggs vary slightly from 7.5,8.0,8.5,
# but this can be easily corrected by just interpolating in the actual logg of the model at the Teff being used.

master = np.genfromtxt('diffusion_timescales/WD_M_0.6.data')
fm80 = interpolate.interp1d(master[:,1],master[:,3],fill_value='extrapolate') # log(M
def Mcvz80(Teff):
    return (10**fm80(Teff))*0.6*cgs.msun # in grams
ft80 = interpolate.interp1d(master[:,1],master[:,10],fill_value='extrapolate') # log(taudiff/yr) as a function of Teff
def tau80(Teff):
    return (10**ft80(Teff))*cgs.year # in years
logg80 = interpolate.interp1d(master[:,1],master[:,2],fill_value='extrapolate')
master = np.genfromtxt('diffusion_timescales/WD_M_0.38.data')

fm75 = interpolate.interp1d(master[:,1],master[:,3],fill_value='extrapolate') # log(M
def Mcvz75(Teff):
    return (10**fm75(Teff))*0.38*cgs.msun # in grams
ft75 = interpolate.interp1d(master[:,1],master[:,10],fill_value='extrapolate') # log(taudiff/yr) as a function of Teff
def tau75(Teff):
    return (10**ft75(Teff))*cgs.year # in years
logg75 = interpolate.interp1d(master[:,1],master[:,2],fill_value='extrapolate')

# Rescale Ca mdot to a total based on Bulk Earth
from EarthComposition import *

def mdot038(XCa,Teff):
    return (XCa*Mcvz75(Teff)/tau75(Teff))/Earth_MassFrac_McD['Ca']
def mdot06(XCa,Teff):
    return (XCa*Mcvz80(Teff)/tau80(Teff))/Earth_MassFrac_McD['Ca']
def mdot09(XCa,Teff):
    return (XCa*Mcvz85(Teff)/tau85(Teff))/Earth_MassFrac_McD['Ca']

def log_mdot(XCa,Teff,logg):
    lgmdot = interpolate.interp1d([logg75(Teff), logg80(Teff), logg85(Teff)],
                                  [np.log10(mdot038(XCa,Teff)), np.log10(mdot06(XCa,Teff)), np.log10(mdot09(XCa,Teff))],
                                  fill_value='extrapolate')
    return lgmdot(logg)

def log_mdots(XCa,Teff,logg):
    logmdot = []
    for i in range(len(Teff)):
        logmdot.append(log_mdot(XCa[i],Teff[i],logg[i]))
    return np.array(logmdot)