Reduced_Chi_Square

# 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(596 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))

#opening the fits file by using hdul command

hdul = fits.open (fitsfiles[155])
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']


fwhm_value=3
sigma_value=1.75


#removing the pixels which are within 3fwhm of both ovi region
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 ,err)) 
pair_1 = list(zip(wave, flux ,err)) 


def filter_o6(p):
  p2 = []
  for e in p:
    if not (a < e[0] < b or c < e[0] < d):
      p2.append(e)

  return p2

pair_new = filter_o6(pair)


#filtering out the pixels that are below a certain level of sigma
#setting up the sigma_value(can find from the KDE plot , upto which value we need to remove the pixels) 

mean_flux,std_flux=norm.fit(flux)
sigma_flux =sigma_value*std_flux
sigma_length_min =mean_flux-sigma_flux

def filter_sigma(p):
  p2 = []
  for e in p:
    i = p.index(e)
    if p[i][1] >= sigma_length_min:
      p2.append(e)
  return p2
pair_new2 = filter_sigma(pair_new)


#the pair of removed pixels 
scatter_pair=[]
for e in pair:
  if e not in pair_new2:
    scatter_pair.append(e)

sw, sf ,se = list(zip(*scatter_pair)) # unpaired the removed pixel to corresponding wavelength flux and error
w, f ,e = list(zip(*pair_new2))   # unpaired the new set of values  to corresponding wavelength flux and error


#fitting the newly obtained pixels 
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


#Defining a function to find the reduced chi square values
def lsq(x,y,n):
  lsq=[]
  coeffs = fit_polynomial(x,y,n)
  poly = np.poly1d(coeffs)
  new_x = x
  new_y = poly(new_x)
  predicted = new_y
  observed = y
  residuals = observed - predicted
  lsq.append(sum(residuals**2)) 
  return lsq 

def rcsq(x,y,z,n):
  reduced_chi_sq = []
  coeff = fit_polynomial(x, y, n)
  poly = np.poly1d(coeff)
  new_x = x
  predicted = poly(new_x)
  observed = y
  uncertainties = z
  residuals = observed - predicted
  chi_square = sum((residuals/uncertainties)**2)
  #print("chi square" ,chi_square)
  degrees_of_freedom= len(observed)-n
  reduced_chi_square = chi_square / degrees_of_freedom
  reduced_chi_sq.append(reduced_chi_square)
  return reduced_chi_sq

#calculating the reduced chi square and least square values for n=1 to n=10
least_square_value=[]
reduced_chi_square=[]
for i in range(1,11):
  least_square_value.append(lsq(w,f,i))
  reduced_chi_square.append(rcsq(w,f,e,i))



#-----------------------------------------------


#plotting the value of least_square_value with different n
N=[1,2,3,4,5,6,7,8,9,10]
fig = plt.figure()
fig.set_size_inches(10, 5)
plt.scatter(N,least_square_value,color="black")
plt.title("Least square value vs n  ")
plt.xlabel("n")
plt.ylabel("Least square value ")
plt.xticks(N)
m= min(least_square_value)
index = (least_square_value.index(m))
for i in least_square_value:
  if i == m:
    plt.scatter(N[index],least_square_value[index],color="red",label="Minimum value of least square value")
plt.legend(bbox_to_anchor=(1.04,1), borderaxespad=0)
plt.savefig('/content/leastsquare.png',bbox_inches='tight')


# Plotting the reduced chi square value vs n 
fig = plt.figure()
fig.set_size_inches(10, 5)
plt.scatter(N,reduced_chi_square,color="black")
plt.title("Reduced chi  square value vs n  ")
plt.xlabel("n")
plt.ylabel("Reduced chi  square value ")
plt.axhline(y=1)
plt.xticks(N)
l= min(reduced_chi_square)
index = (reduced_chi_square.index(l))
for k in reduced_chi_square:
  if k == l:
    plt.scatter(N[index],reduced_chi_square[index],color="red",label="Minimum value of reduced chi square value")
plt.legend(bbox_to_anchor=(1.04,1), borderaxespad=0)
plt.savefig('/content/Reduced chi square.png',bbox_inches='tight')

'''
#just to chcek the plot of the continum and spectra
coeffs1 = fit_polynomial(w, f,8)
poly1 = np.poly1d(coeffs1)
new_x1 = w
new_y1 = poly1(new_x1)
fig = plt.figure()
fig.set_size_inches(20, 5)
plt.plot(wave,flux,color='black')
plt.plot(new_x1,new_y1,"blue",label='Newly fitted continum ')
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.scatter(sw,sf,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 , n=8 ")
plt.legend(bbox_to_anchor=(1.04,1), borderaxespad=0)
plt.savefig('/content/continum fitted plot.png',bbox_inches='tight')
'''