-
Notifications
You must be signed in to change notification settings - Fork 0
/
javelin_sdss_quasars_u.py
80 lines (53 loc) · 1.77 KB
/
javelin_sdss_quasars_u.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
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 18 21:16:19 2014
@author: suberlak
Runs javelin on SDSS Quasars, a selected band
(filenames start as u_ , g_, etc.)
from QSO_SDSS_JAV
"""
import numpy as np
from javelin.zylc import get_data
from javelin.lcmodel import Cont_Model
dir_choice=['QSO_try2/', 'QSO_S82/', 'QSO_SDSS_JAV/','QSO_SDSS_chains/']
dir_input=dir_choice[2]
dir_output=dir_choice[3]
band= 'u_band'
names=np.loadtxt(dir_input+band+'.ls',dtype=str)
"""
Restrict to only those lightcurves that have at least 10 observed points
BEGINNING OF CODE
"""
cond_length = np.empty_like(names,dtype=bool)
min_length = 10
for i in range(len(names)): #
test = np.loadtxt(dir_input+names[i])
print names[i], len(test)
if len(test) > min_length:
cond_length[i] = True
else:
cond_length[i] = False
names_upd = names[cond_length]
print 'Test for length 2'
for i in range(len(names_upd)): #
test = np.loadtxt(dir_input+names_upd[i])
print names_upd[i], len(test)
too_short_num = len(np.where(cond_length == False)[0])
print 'We have', too_short_num,' quasars whose lightcurves were shorter than ', min_length
file_short_lc = dir_input+band+'_too_short_lc.ls'
cond_short = -cond_length
files_list = names[cond_short]
np.savetxt(file_short_lc ,files_list, fmt="%s")
print 'Their names were saved to a file', file_short_lc
"""
END OF CODE
"""
for i in range(3118,len(names_upd)): # len(names_upd)
filename=dir_input+names_upd[i]
print '\nWorking file', filename, i+1, ' out of ', len(names_upd)
data = get_data(filename,names=["Continuum"])
cont=Cont_Model(data)
start=len(dir_input)
chain_name = dir_output+'ch_'+filename[start:]+'_chain.dat'
print chain_name
cont.do_mcmc(fchain=chain_name)