-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgacrux_extraction.py
executable file
·142 lines (114 loc) · 5.65 KB
/
gacrux_extraction.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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
"""A script testing the extraction pipeline of RHEA
Steps
1) Initialise Format, Extractor and RadialVelocity
2) Define file paths for science, flat and dark frames
3) Extract/import spectra
4) Create/import reference spectra
5) Calculate radial velocities
6) Plot radial velocities
"""
import numpy as np
try:
import pyfits
except:
import astropy.io.fits as pyfits
import pymfe
import glob
from astropy.coordinates import SkyCoord
from astropy import units as u
#===============================================================================
# Parameters/Constants/Variables/Initialisation
#===============================================================================
# Constants/Variables
do_bcor = True
med_cut = 0.6
# Specified header parameters
xbin = 2
ybin = 1
exptime = 1200
badpixel_mask= pyfits.getdata('/priv/mulga1/jbento/rhea2_data/badpix.fits')
badpix=np.where(badpixel_mask==1)
# Initialise objects
rhea2_format = pymfe.rhea.Format()
rhea2_extract = pymfe.Extractor(rhea2_format, transpose_data=False,
badpixmask=badpix)
xx, wave, blaze = rhea2_format.spectral_format()
rv = pymfe.rv.RadialVelocity()
#===============================================================================
# File paths (Observations, Flats and Darks, save/load directories)
#===============================================================================
# Science Frames
star = "gammaCrucis"
base_path = "/priv/mulga1/jbento/rhea2_data/gammaCrucis/"
# Find all Gacrux files and sort by observation date in MJD
all_files = np.array(glob.glob(base_path + "2015*/*" + star + "*.fit*"))
sorted = np.argsort([pyfits.getheader(e)['MJD'] for e in all_files])
all_files = all_files[sorted]
files = []
# Only consider files that have the same exposure time and correct binning
for f in all_files:
fits = pyfits.open(f)
header = fits[0].header
x_head = header["XBINNING"]
y_head = header["YBINNING"]
exp_head = header["EXPTIME"]
if x_head == xbin and y_head == ybin and exp_head == exptime:
files.append(f)
fits.close()
# Flats and Darks
dark_path = base_path + "Dark frames/Masterdark_target_" + str(exptime) +".fits"
star_dark = pyfits.getdata(dark_path)
flat_dark = pyfits.getdata(base_path + "Dark frames/Masterdark_flat.fit")
# Note: this particular flat was chosen as it has an exposure time of 3 seconds,
# the same length as the flat dark that will be used to correct it
flat_path = base_path + "20150527/20150527_Masterflat.fit"
flat_files = [flat_path]*len(files)
# Extracted spectra output
#out_path = "/priv/mulga1/arains/Gacrux_Extracted_" + str(exptime) + "/"
out_path = "/priv/mulga1/arains/Gacrux_Extracted_Both/"
extracted_files = np.array(glob.glob(out_path + "*" + star + "*extracted.fits"))
# Sort to account for files not being labelled with MJD
sorted = np.argsort([pyfits.getheader(e)['JD'] for e in extracted_files])
extracted_files = extracted_files[sorted]
# RV csv output
base_rv_path = out_path + star
ref_path = "/priv/mulga1/arains/Gacrux_Extracted_Both/reference_spectrum_10_gammaCrucis.fits"
#===============================================================================
# Extract and save spectra/load previously extracted spectra
#===============================================================================
# Extract spectra ("wave" removed)
# OPTION 1: Extract and save spectra
#fluxes, vars, bcors, mjds = rv.extract_spectra(files, rhea2_extract,
# star_dark=star_dark,
# flat_files=flat_files,
# flat_dark=flat_dark,
# do_bcor=do_bcor,
# ra_dec_hr=True)
# Save spectra (Make sure to save "wave" generated from rhea2_format)
#rv.save_fluxes(files, fluxes, vars, bcors, wave, mjds, out_path)
# OPTION 2: Load previously extracted spectra
#extracted_files = extracted_files[7:44]
fluxes, vars, wave, bcors, mjds = rv.load_fluxes(extracted_files)
#===============================================================================
# Create and save/import reference spectrum
#===============================================================================
# Number of frames to use for reference spectrum
n = 10
# OPTION 1: Create and save a new reference spectrum
#wave_ref, ref_spect = rv.create_ref_spect(wave, fluxes[:n,:,:], vars[:n,:,:],
# bcors[:n], med_cut=med_cut)
#rv.save_ref_spect(files[:n], ref_spect, vars[:n,:,:], wave_ref, bcors[:n],
# mjds[:n], out_path, star)
# OPTION 2: Import a pre-existing reference spectrum
ref_spect, vars_ref, wave_ref, bcors_ref, mjds_ref = rv.load_ref_spect(ref_path)
#===============================================================================
# Calculate, save and plot radial velocities
#===============================================================================
# Calculate RVs
rvs, rv_sigs = rv.calculate_rv_shift(wave_ref, ref_spect, fluxes, vars, bcors,
wave)
nf = fluxes.shape[0]
nm = fluxes.shape[1]
bcor_rvs = rvs + bcors.repeat(nm).reshape( (nf, nm) )
# Save RVs
rv.save_rvs(rvs, rv_sigs, bcors, mjds, bcor_rvs, base_rv_path)