-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgmosplan.py
211 lines (199 loc) · 10.2 KB
/
gmosplan.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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
from __future__ import division
import numpy
import pyfits
import pylab
import tools
import sys
def create_obstab(prefix, pixscale, objid, ra, dec, x_ccd, y_ccd, mag,
priority, slitsize_x, slitsize_y, slittilt,
slitpos_y='Default'):
'''
'''
# Write the header cards
prihdr = pyfits.Header()
prihdr['COMMENT'] = 'GMOS Object Table for use with Gemini MOS Mask Preparation Software (GMMPS).'
prihdr['PIXSCALE'] = pixscale
prihdr['SLITMASK'] = prefix
prihdu = pyfits.PrimaryHDU(header=prihdr)
# Convert RA in degree to hours
ra_hr = ra/15.0
# populate 1D arrays
n_obj = numpy.size(ra)
one_array = numpy.ones(n_obj)
slitsize_x *= one_array
slittilt *= one_array
if slitpos_y == 'Default':
slitpos_y = numpy.zeros(n_obj)
# For the alignment star objects remove the slit size and alignment info
mask_star = priority == 0
slitsize_x[mask_star] = numpy.nan
slitsize_y[mask_star] = numpy.nan
slittilt[mask_star] = numpy.nan
slitpos_y[mask_star] = numpy.nan
# Define the columns and their formatting and content
col_id = pyfits.Column(name='ID', format='K', array=objid)
col_ra = pyfits.Column(name='RA', format='D', array=ra_hr)
col_dec = pyfits.Column(name='DEC', format='D', array=dec)
col_x_ccd = pyfits.Column(name='x_ccd', format='D', array=x_ccd)
col_y_ccd = pyfits.Column(name='y_ccd', format='D', array=y_ccd)
col_mag = pyfits.Column(name='mag', format='E', array=mag)
col_priority = pyfits.Column(name='priority', format='I', array=priority)
col_slitsize_x = pyfits.Column(name='slitsize_x', format='E', array=slitsize_x)
col_slitsize_y = pyfits.Column(name='slitsize_y', format='E', array=slitsize_y)
col_slittilt = pyfits.Column(name='slittilt', format='E', array=slittilt)
col_slitpos_y = pyfits.Column(name='slitpos_y', format='E', array=slitpos_y)
# define column order
cols = pyfits.ColDefs([col_id, col_ra, col_dec, col_x_ccd, col_y_ccd,
col_mag, col_priority, col_slitsize_x,
col_slitsize_y, col_slittilt, col_slitpos_y])
# populate table hdu
tbhdu = pyfits.BinTableHDU.from_columns(cols)
# create list containing header and data
thdulist = pyfits.HDUList([prihdu, tbhdu])
filename = prefix + '_objtab.fits'
thdulist.writeto(filename, clobber=True)
print 'gmosplan: Finished writing object table fits file.'
def plotcoverage(redshift,lambda_central,filename=None):
'''
Creates a plot of common spectral features in the observed frame. Along with
a spectra coverage region for the DEIMOS 1200 line grating for the given
central wavelength.
INPUT:
redshift = [float] redshift of the cluster
lambda_central = [float] central wavelength for the DEIMOS 1200 line grating
filename = [None or string] if not None then an image will be saved with
filename
OUTPUT:
-- displays a figure
-- if filename given will save figure as a png file
'''
fig = pylab.figure(figsize=(20,4.5))
pylab.xlim((lambda_central-2300, lambda_central+2300))
xl = pylab.xlim()
yl = (0,1)
lcl_low = lambda_central-1300
lcl_up = lambda_central+1300
#Plot the wavelength coverage
pylab.plot((lambda_central,lambda_central),yl,'-b',linewidth=3)
pylab.plot((lcl_low,lcl_low),yl,
'-b',alpha=0.5,linewidth=3)
pylab.plot((lcl_up,lcl_up),yl,
'-b',alpha=0.5,linewidth=3)
pylab.fill_between(pylab.arange(lcl_low,lcl_up+100,100),yl[0],yl[1],
facecolor='blue',alpha=0.25)
#At the ends of the coverage window show the +/- range due to where the
#slit is placed on the mask
pylab.plot((lcl_low+411,lcl_low+411),yl,
'--b',alpha=0.5,linewidth=3)
pylab.plot((lcl_low-411,lcl_low-411),yl,
'-.b',alpha=0.5,linewidth=3)
pylab.plot((lcl_up+411,lcl_up+411),yl,
'--b',alpha=0.5,linewidth=3)
pylab.plot((lcl_up-411,lcl_up-411),yl,
'-.b',alpha=0.5,linewidth=3)
#Plot common spectral lines
x_Lyb = 1025.7*(1+redshift)
x_Lya = 1215.7*(1+redshift)
x_CIV = 1549.1*(1+redshift)
x_AlIII = 1858.7*(1+redshift)
x_FeII = 2600*(1+redshift)
x_MgII = 2799.8*(1+redshift)
x_MgI = 2852*(1+redshift)
x_OII = 3727.61*(1+redshift)
x_CalK = 3933.667*(1+redshift)
x_CalH = 3968.472*(1+redshift)
x_Hd = 4101.74*(1+redshift)
x_Gband = 4305*(1+redshift)
x_Hg = 4340.47*(1+redshift)
x_Hb = 4861.33*(1+redshift)
x_OIII_1 = 4960.3*(1+redshift)
x_OIII_2 = 5008.24*(1+redshift)
x_Mgb = 5176*(1+redshift)
x_FeI = 5269*(1+redshift)
x_NaD = 5893*(1+redshift)
x_NII_1 = 6548.06*(1+redshift)
x_Ha = 6562.799*(1+redshift)
x_NII_2 = 6585.2*(1+redshift)
x_SII = 6725.5*(1+redshift)
# Plot dashed lines at the respective redshifts
pylab.plot((x_Lyb,x_Lyb),yl,'--k')
pylab.plot((x_Lya,x_Lya),yl,'--k')
pylab.plot((x_CIV,x_CIV),yl,'--k')
pylab.plot((x_AlIII,x_AlIII),yl,'--k')
pylab.plot((x_FeII,x_FeII),yl,'--k')
pylab.plot((x_MgII,x_MgII),yl,'--k')
pylab.plot((x_MgI,x_MgI),yl,'--k')
pylab.plot((x_OII,x_OII),yl,'--k')
pylab.plot((x_CalK,x_CalK),yl,'--k')
pylab.plot((x_CalH,x_CalH),yl,'--k')
pylab.plot((x_Hd,x_Hd),yl,'--k')
pylab.plot((x_Gband,x_Gband),yl,'--k')
pylab.plot((x_Hg,x_Hg),yl,'--k')
pylab.plot((x_Hb,x_Hb),yl,'--k')
pylab.plot((x_OIII_1,x_OIII_1),yl,'--k')
pylab.plot((x_OIII_2,x_OIII_2),yl,'--k')
pylab.plot((x_Mgb,x_Mgb),yl,'--k')
pylab.plot((x_FeI,x_FeI),yl,'--k')
pylab.plot((x_NaD,x_NaD),yl,'--k')
pylab.plot((x_NII_1,x_NII_1),yl,'--k')
pylab.plot((x_Ha,x_Ha),yl,'--k')
pylab.plot((x_NII_2,x_NII_2),yl,'--k')
pylab.plot((x_SII,x_SII),yl,'--k')
labeloff = 0.5
pylab.text(lambda_central, labeloff*(yl[0]+yl[1]),
'$\lambda_\mathrm{central}$'+'={0}'.format(lambda_central),
horizontalalignment='right',verticalalignment='center',
rotation='vertical',fontsize=16)
if x_Lyb > xl[0] and x_Lyb < xl[1]:
pylab.text(x_Lyb, labeloff*(yl[0]+yl[1]), 'Ly-beta', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Lya > xl[0] and x_Lya < xl[1]:
pylab.text(x_Lya, labeloff*(yl[0]+yl[1]), 'Ly-alpha', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_CIV > xl[0] and x_CIV < xl[1]:
pylab.text(x_CIV, labeloff*(yl[0]+yl[1]), 'C IV', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_AlIII > xl[0] and x_AlIII < xl[1]:
pylab.text(x_AlIII, labeloff*(yl[0]+yl[1]), 'Al III', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_FeII > xl[0] and x_FeII < xl[1]:
pylab.text(x_FeII, labeloff*(yl[0]+yl[1]), 'Fe II', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_MgII > xl[0] and x_MgII < xl[1]:
pylab.text(x_MgII, labeloff*(yl[0]+yl[1]), 'Mg II', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_MgI > xl[0] and x_MgI < xl[1]:
pylab.text(x_MgI, labeloff*(yl[0]+yl[1]), 'Mg I', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_OII > xl[0] and x_OII < xl[1]:
pylab.text(x_OII, labeloff*(yl[0]+yl[1]), '[O II]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_CalK > xl[0] and x_CalK < xl[1]:
pylab.text(x_CalK, labeloff*(yl[0]+yl[1]), 'Cal K', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_CalH > xl[0] and x_CalH < xl[1]:
pylab.text(x_CalH, labeloff*(yl[0]+yl[1]), 'Cal H', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Hd > xl[0] and x_Hd < xl[1]:
pylab.text(x_Hd, labeloff*(yl[0]+yl[1]), 'Hd', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Gband > xl[0] and x_Gband < xl[1]:
pylab.text(x_Gband, labeloff*(yl[0]+yl[1]), 'G-band', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Hg > xl[0] and x_Hg < xl[1]:
pylab.text(x_Hg, labeloff*(yl[0]+yl[1]), 'Hg', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Hb > xl[0] and x_Hb < xl[1]:
pylab.text(x_Hb, labeloff*(yl[0]+yl[1]), 'Hb', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_OIII_1 > xl[0] and x_OIII_1 < xl[1]:
pylab.text(x_OIII_1, labeloff*(yl[0]+yl[1]), '[OIII]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_OIII_2 > xl[0] and x_OIII_2 < xl[1]:
pylab.text(x_OIII_2, labeloff*(yl[0]+yl[1]), '[OIII]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Mgb > xl[0] and x_Mgb < xl[1]:
pylab.text(x_Mgb, labeloff*(yl[0]+yl[1]), 'Mg I(b)', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_FeI > xl[0] and x_FeI < xl[1]:
pylab.text(x_FeI, labeloff*(yl[0]+yl[1]), 'Fe I', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_NaD > xl[0] and x_NaD < xl[1]:
pylab.text(x_NaD, labeloff*(yl[0]+yl[1]), 'Na I (D)', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_NII_1 > xl[0] and x_NII_1 < xl[1]:
pylab.text(x_NII_1, (labeloff-0.1)*(yl[0]+yl[1]), '[NII]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_Ha > xl[0] and x_Ha < xl[1]:
pylab.text(x_Ha, labeloff*(yl[0]+yl[1]), 'Ha', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_NII_2 > xl[0] and x_NII_2 < xl[1]:
pylab.text(x_NII_2, (labeloff+0.1)*(yl[0]+yl[1]), '[NII]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
if x_SII > xl[0] and x_SII < xl[1]:
pylab.text(x_SII, labeloff*(yl[0]+yl[1]), '[SII]', horizontalalignment='right',verticalalignment='center', rotation='vertical')
pylab.xlim(xl)
frame1 = pylab.gca()
frame1.axes.get_yaxis().set_visible(False)
pylab.xlabel('$\lambda_{observed}$',fontsize=18)
if filename:
pylab.savefig(filename)
pylab.show()