-
Notifications
You must be signed in to change notification settings - Fork 2
/
compile_transpec.py
66 lines (60 loc) · 2.91 KB
/
compile_transpec.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
import numpy as np
import glob
import argparse
import utils
import pickle
import os
parser = argparse.ArgumentParser()
# This parses in the option file:
parser.add_argument('-ofile',default=None)
parser.add_argument('--CMC', dest='CMC', action='store_true')
parser.set_defaults(CMC=False)
args = parser.parse_args()
CMC = args.CMC
ofile = args.ofile
# Read input file:
datafile, ld_law, idx_time, all_comps, P, Psd, \
a, asd, pmean, psd, b, bsd, t0,\
t0sd, fixed_eccentricity, ecc, eccsd, \
omega, omegasd = utils.read_optfile(ofile)
######################################
target,pfilename = datafile.split('/')
if not CMC:
out_folder = 'outputs/'+datafile.split('.')[0]+'/wavelength'
else:
out_folder = 'outputs/'+datafile.split('.')[0]+'/wavelength-cmc'
out_ofolder = 'outputs/'+datafile.split('.')[0]
if not CMC:
fout = open(out_ofolder+'/transpec.dat','w')
else:
fout = open(out_ofolder+'/transpec_cmc.dat','w')
fout.write('# Wav (Angstroms) \t Rp/Rs \t Rp/RsErrUp \t Rp/RsErrDown \t Depth (ppm) \t Depthup (ppm) \t DepthDown (ppm)\n')
wbinlist = glob.glob(out_folder+'/wbin*')
wis = np.array([])
for wbinfolders in wbinlist:
wis = np.append(wis,int(wbinfolders.split('wbin')[-1]))
wis = wis[np.argsort(wis)]
wis = wis.astype(int)
for wi in wis:
print 'Working on wbin ',wi,'...'
try:
out = pickle.load(open(out_folder+'/wbin'+str(wi)+'/BMA_posteriors.pkl','rb'))
wbins = out['wbin']
p = out['p']
v,vup,vdown = utils.get_quantiles(p)
vd,vupd,vdownd = utils.get_quantiles((p**2)*1e6)
vld,vupld,vdownld = utils.get_quantiles(out['q1'])
if ld_law == 'linear':
try:
fout.write('{0:.2f} \t {1:.2f} \t {2:.10f} \t {3:.10f} \t {4:.10f} \t {5:.10f} \t {6:.10f} \t {7:.10f} \t {8:.10f} \t {9:.10f} \t {10:.10f} \n'.format(wbins[0],wbins[1],v,vup-v,v-vdown,vd,vupd-vd,vd-vdownd,vld,vupld-vld,vld-vdownld))
except:
fout.write('{0:d} \t {1:d} \t {2:.10f} \t {3:.10f} \t {4:.10f} \t {5:.10f} \t {6:.10f} \t {7:.10f} \t {8:.10f} \t {9:.10f} \t {10:.10f} \n'.format(wbins,wbins,v,vup-v,v-vdown,vd,vupd-vd,vd-vdownd,vld,vupld-vld,vld-vdownld))
else:
vld2,vupld2,vdownld2 = utils.get_quantiles(out['q2'])
try:
fout.write('{0:.2f} \t {1:.2f} \t {2:.10f} \t {3:.10f} \t {4:.10f} \t {5:.10f} \t {6:.10f} \t {7:.10f} \t {8:.10f} \t {9:.10f} \t {10:.10f} \t {11:.10f} \t {12:.10f} \t {13:.10f} \n'.format(wbins[0],wbins[1],v,vup-v,v-vdown,vd,vupd-vd,vd-vdownd,vld,vupld-vld,vld-vdownld,vld2,vupld2-vld2,vld2-vdownld2))
except:
fout.write('{0:d} \t {1:d} \t {2:.10f} \t {3:.10f} \t {4:.10f} \t {5:.10f} \t {6:.10f} \t {7:.10f} \t {8:.10f} \t {9:.10f} \t {10:.10f} \t {11:.10f} \t {12:.10f} \t {13:.10f}\n'.format(wbins,wbins,v,vup-v,v-vdown,vd,vupd-vd,vd-vdownd,vld,vupld-vld,vld-vdownld,vld2,vupld2-vld2,vld2-vdownld2))
except:
print 'No data (yet)'
fout.close()