import numpy as np
from pylab import *
from pylab import loadtxt
from scipy import *
from astropy.io import fits
import matplotlib.pyplot as plt
import scipy.constants as const
import scipy.stats as stats
from scipy import optimize

#Open our fits and turn them into 2D numpy arrays
FPtau = fits.open("/home/kmiskove/Documents/lab-5/FP_Tau.cleaned.fits")
FPt_w = FPtau[0].data

#Make arrays for wavelengths and flux density and plot them
FPt_wl = FPt_w[0]
FPt_fd = FPt_w[1]
plt.plot(FPt_wl, FPt_fd)
plt.ylabel('Flux density')
plt.xlabel('Wavelength (microns)')

#But our plot isn't that coherent, so let's try plotting less of it:
plt.axis([2.3, 2.31, 9.4743346e-15, 1.3214704e-14])
#plt.show()

#We pick a little dip and try to fit a gaussian to it now
plt.axis([2.3025, 2.3035, 9.4743346e-15, 1.2e-14,])

#Define our gaussian function
def gaussian(x,a,b,c):
    e = np.e
    t = (-(x-b)**2)
    s = (2*c**2)
    m = e**(t/s)
    g = (1-(a*m))
    return g

#By manually fiddling with the following and looking for when it outputs values closest to 2.3025 and 2.3035, we can manually set our limits:
#print(FPt_wl[24759:24872])
#we found our indices should be from 24759 to 24872.

#Now let's define our a, b, and c values and plot our gaussian
a = 0.1
b = 2.30297
c = 0.00015

guess = gaussian(FPt_wl[24759:24872],a,b,c)
d = 1.12e-14

plt.plot(FPt_wl[24759:24872], FPt_fd[24759:24872])
plt.plot(FPt_wl[24759:24872], guess*d)
plt.ylabel('Flux density')
plt.xlabel('Wavelength (microns)')
plt.show()

#Now let's do the standard deviation velocity thing...

c = 3e8
velocity = c*(2.3035-2.3025/2.3035)
print(velocity)
print("This is the velocity of HK Tau in m/s")

#Sun's diameter: 10.3e9
d = 10.3e9
#convert to circumference:
cir = 32358404332
rot = cir/velocity
print(rot)
print("This is the star's rotational period in seconds.")
print("The sun's rotational period is 2114208 seconds, so this is absolutely insane.")