-
Notifications
You must be signed in to change notification settings - Fork 0
/
code for mean and median stacked spectra
87 lines (67 loc) · 2.93 KB
/
code for mean and median stacked spectra
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
#The code to find the mean and median stacked spectra for the PURE OVI systems
#importing all the important libraries
from astropy.io import fits
from astropy.table import Table
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import os
from statistics import *
# mount drive
from google.colab import drive
drive.mount('/content/gdrive')
#getting all the fits files of the pure systems from the path by using os.path
path="/content/gdrive/MyDrive/OVI PROJECT DATA/FITS DATA/Pure OVI systems"
fitsfiles = []
for root, dirs, files in os.walk(path):
for file in files:
fitsfiles.append(os.path.join(root,file))
#setting a value for the wavelength range
lower_bound = 1028.986328125
upper_bound = 1040.813232421875
w_int = np.arange(lower_bound,upper_bound,.03)
#code for rebinning
rebinned_flux=[]
for i in fitsfiles:
new_flux=[]
hdul = fits.open (i)
flux= hdul[1].data['FLUX']
waves =hdul[1].data['WAVE']
cont=hdul[1].data['CONT']
err= hdul[1].data['ERR']
f = interp1d(waves,((flux/cont)),fill_value ='extrapolate')
f_int = f(w_int)
new_flux = list(f_int)
rebinned_flux.append(new_flux)
#we need the first elements of all members of rebinned_flux to be stacked.The following code does this
arrays = [np.array(x) for x in rebinned_flux]
mean_flux = [np.mean(k) for k in zip(*arrays)]
median_flux = [np.median(k) for k in zip(*arrays)]
#----------------------------------------PLOTTING CODES
o1_wave = 1031.93
o2_wave = 1037.62
dw1 = (1.665*50*3*o1_wave)/(3*(10**5))
dw2 = (1.665*50*3*o2_wave)/(3*(10**5))
#plotting the figures
fig = plt.figure()
fig.set_size_inches(20, 5)
plt.title('Medain stacked spectra')
plt.axvline(x=(1031.93) ,color='blue', linestyle='--' , label = 'OVI line(1031.93)' , ymin=0, ymax=0.99)
plt.axvline(x=(1037.62) , color='blue', linestyle='--' , label ='OVI line(1037.62' , ymin=0, ymax=0.99 )
plt.axhline(y=1, color='green', linestyle='-')
plt.axvspan(1031.93-(dw1/2),1031.93+(dw1/2), facecolor='green', alpha=0.5,label='3fwhm range of OVI (b=50km/s)')
plt.axvspan(1037.62-(dw2/2),1037.62+(dw2/2), facecolor='red', alpha=0.5,label='3fwhm range of OVI (b=50km/s)')
plt.text(1033,1.01,"median stacked spectra ",size='x-large',color ='black')
plt.legend()
plt.plot(w_int,(median_flux),color='black')
fig = plt.figure()
fig.set_size_inches(20, 5)
plt.title('mean stacked spectra')
plt.axvline(x=(1031.93) ,color='blue', linestyle='--' , label = 'OVI line(1031.93)' , ymin=0, ymax=0.99)
plt.axvline(x=(1037.62) , color='blue', linestyle='--' , label ='OVI line(1037.62' , ymin=0, ymax=0.99 )
plt.axhline(y=1, color='green', linestyle='-')
plt.axvspan(1031.93-(dw1/2),1031.93+(dw1/2), facecolor='green', alpha=0.5,label='3fwhm range of OVI (b=50km/s)')
plt.axvspan(1037.62-(dw2/2),1037.62+(dw2/2), facecolor='red', alpha=0.5,label='3fwhm range of OVI (b=50km/s)')
plt.text(1033,1.01,"mean stacked spectra ",size='x-large',color ='black')
plt.legend()
plt.plot(w_int,(mean_flux),color='black')