o6stack/Code for continum fitting at main · RomeoPallikkara/o6stack · GitHub

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
# mount drive
from google.colab import drive
drive.mount('/content/gdrive')

#importing all important libraries for the code

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 *
from scipy.stats import norm
import seaborn as sns
#setting the path to the folder containg the fits files of the systems
path ="/content/gdrive/MyDrive/OVI PROJECT DATA/FITS DATA/Pure OVI systems"

#using the os.walk  command in python to get only the fits files inside the directory and subdirectory

fitsfiles = []
for root, dirs, files in os.walk(path):
	for file in files:
		fitsfiles.append(os.path.join(root,file))
#now the list contains the path to each fits files

#opening the fits file by using hdul command

hdul = fits.open (fitsfiles[298])
cont= hdul[1].data['CONT']
wave =hdul[1].data['WAVE']
flux= hdul[1].data['FLUX']
err=hdul[1].data['ERR']
errc=hdul[1].data['ERRC']
econt=hdul[1].data['ECONT']
expt=hdul[1].data['EXPT']
#The value of fwhm(b=50km/s) upto which in the region of OVI chunk we need to remove the pixels
fwhm_value=3
o1_wave = 1031.912
o2_wave = 1037.613
dw1 = (1.665*50*fwhm_value*o1_wave)/(3*(10**5))
dw2 = (1.665*50*fwhm_value*o2_wave)/(3*(10**5))

a=(o1_wave-(dw1/2))
b=(o1_wave+(dw1/2))
c=(o2_wave-(dw2/2))
d=(o2_wave+(dw2/2))

pair = list(zip(wave, flux))
pair_1 = list(zip(wave, flux))


#defined a function to remove the OVI chunk regions
def filter_ovi(x):
  y = []
  for e in x:
    if not (a < e[0] < b or c < e[0] < d):
      y.append(e)

  return y


#created the new pair of wavelength and flux with the pixels in region of 3fwhm of ovi removed.
pair_new = filter_ovi(pair_1)
flux_new = [pair[pair.index(i)][1] for i in pair]

#defined a function to remove the flux below a certain level of sigma and create new pair
def filter_sigma(x, slm):
  x2 = []
  for e in x:
    i = x.index(e)
    if x[i][1] >= slm:
      x2.append(e)
  return x2

number_of_ovi_deletions=len(pair)-len(pair_new)
print("number_of_ovi_deletions = ",number_of_ovi_deletions)
#creating a while loop in which we are looping untill aperticular number of pixels are removed
scatter_pair=[]
sigma_value=4

N=number_of_ovi_deletions +10 # setting up a value of sigma value deletions (set it up as 40 pixels to be removed)
while len(scatter_pair) < N and sigma_value > 0:

    mean_flux,std_flux=norm.fit(flux_new)

    sigma_flux =sigma_value*std_flux
    sigma_length_min =mean_flux-sigma_flux

    pair_new = filter_sigma(pair_new, sigma_length_min)

    for e in pair:
      if e not in pair_new and e not in scatter_pair:
        scatter_pair.append(e)

    sigma_value-=0.1

    flux_new = [pair_new[pair_new.index(i)][1] for i in pair_new]

    print(sigma_value, len(scatter_pair))
#defined a function to fit the curve

def fit_polynomial(x, y, n):
  # Create a matrix of powers of x, from 0 to n
  X = np.vander(x, n+1)

  # Use the least squares method to find the coefficients of the polynomial
  coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)

  # Return the coefficients of the polynomial
  return coeffs

#------------------------------------------------PLOTTING CODE


#seperating out the wavelength and flux from the zip we obtained
scatter_wave, scatter_flux = list(zip(*scatter_pair))
new_wave, new_flux = list(zip(*pair_new))

#calling the polynomial fitting function
coeffs = fit_polynomial(new_wave, new_flux, 5)
poly = np.poly1d(coeffs)
new_x = new_wave
new_y = poly(new_x)


#plotting the results
fig = plt.figure()
fig.set_size_inches(20, 5)
plt.plot(wave,flux,color='black')
plt.axvline(x=(1031.93) ,color='green', linestyle='--' , label = 'OVI line(1031.93)' , ymin=0, ymax=0.99)
plt.axvline(x=(1037.62) , color='red', linestyle='--' , label ='OVI line(1037.62' , ymin=0, ymax=0.99 )
plt.plot(wave,cont,color='red',label='Danfoth continum ')
plt.plot(new_x,new_y,color='blue',label='Obtained continum')
plt.scatter(scatter_wave, scatter_flux,color='blue',label='pixels removed')
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.title("Continum fitted plot")
plt.legend(bbox_to_anchor=(1.04,1), borderaxespad=0)