forked from StuartLittlefair/lfit_python
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlimbdark.py
executable file
·83 lines (69 loc) · 3.28 KB
/
limbdark.py
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
'''Uses the output fluxes of a White Dwarf from an MCMC chain, and uses them
to calculate the limb darkening coefficient that should be used when fitting
that body. Data is typically not sensitive enough to this to fit it
independantly, so a few iterations through this script should aid in the
science.
'''
from __future__ import print_function
from builtins import input, zip
import os
import numpy as np
from scipy.interpolate import (RectBivariateSpline, SmoothBivariateSpline,
interp2d)
def ld (band,logg,teff,law='linear'):
assert band in ['u','g','r','i','z']
assert law in ['linear','quad','sqr']
filename = os.path.split(os.path.abspath(__file__))[0]
filename += '/Gianninas13/ld_coeffs_%s.txt' % band
data=np.loadtxt(filename)
x=data[:,0] #logg (grid increments through all teffs at a single logg, then +logg)
y=data[:,1] #teff
t1 = np.unique(y) # unique temps
g1 = np.unique(x) # unique loggs
nt = len(t1)
ng = len(g1)
z0=data[:,2] #linear ld coefficient
z1=data[:,3] # first quad term
z2=data[:,4] # 2nd quad term
z3=data[:,5] # first square-root term
z4=data[:,6] # second square-root term
#func = SmoothBivariateSpline(x,y,z)
#return func(logg,teff)[0]
if law == 'linear':
func = RectBivariateSpline(g1,t1,z0.reshape((ng,nt)),kx=3,ky=3)
return func(logg,teff)[0,0]
elif law == 'quad':
funca = RectBivariateSpline(g1,t1,z1.reshape((ng,nt)),kx=3,ky=3)
funcb = RectBivariateSpline(g1,t1,z2.reshape((ng,nt)),kx=3,ky=3)
return (funca(logg,teff)[0,0],funcb(logg,teff)[0,0])
elif law == 'sqr':
funca = RectBivariateSpline(g1,t1,z3.reshape((ng,nt)),kx=3,ky=3)
funcb = RectBivariateSpline(g1,t1,z4.reshape((ng,nt)),kx=3,ky=3)
return (funca(logg,teff)[0,0],funcb(logg,teff)[0,0])
def main():
logg, gerr = input('> Give log g and error: ').split()
teff, terr = input('> Give eff. temp. and error: ').split()
logg = float(logg); gerr = float(gerr)
teff = float(teff); terr = float(terr)
gvals=np.random.normal(loc=logg,scale=gerr,size=100)
tvals=np.random.normal(loc=teff,scale=terr,size=100)
#ldvals = []
#for g,t in zip(gvals,tvals):
# ldvals.extend( ld('i',g,t) )
print('-------------------')
for band in ['u','g','r','i','z']:
ldvals = [ld(band,g,t) for g,t in zip(gvals,tvals)]
print('%s band LD coeff = %f +/- %f' % (band, np.median(ldvals),np.std(ldvals)))
ldvals = [ld(band,g,t,law='quad') for g,t in zip(gvals,tvals)]
# unpack list of tuples into two lists
a, b = list(zip(*ldvals)) # use splat operator to expand list into positional arguments
print('%s band quad coeff a = %f +/- %f' % (band, np.median(a),np.std(a)))
print('%s band quad coeff b = %f +/- %f' % (band, np.median(b),np.std(b)))
ldvals = [ld(band,g,t,law='sqr') for g,t in zip(gvals,tvals)]
# unpack list of tuples into two lists
a, b = list(zip(*ldvals)) # use splat operator to expand list into positional arguments
print('%s band sqr coeff d = %f +/- %f' % (band, np.median(a),np.std(a)))
print('%s band sqr coeff f = %f +/- %f' % (band, np.median(b),np.std(b)))
print('-------------------')
if __name__ == "__main__":
main()