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

#Let's plot TTau in microns and
#First column is wavelength, second is flux density, third is something we ignore.
TTmic = loadtxt('/home/kmiskove/Documents/TTau_merged_microns.txt', usecols=(0,1))
TTw, TTf = (TTmic[:,0],TTmic[:,1])
M2V = loadtxt('/home/kmiskove/Documents/M2V_Gl806.txt', skiprows=45, usecols=(0,1))
M2Vw, M2Vf = (M2V[:,0],M2V[:,1])

#Let's normalize the fluxes:
TTfn = TTf/median(TTf)
M2fn = M2Vf/median(M2Vf)

#This is where we specify our Av, T, and r values.
Av = 0.3 #veiling
T = 1500 #temp
r = 0.7 #scaling
#And let's define some constants:
h = 6.26e-34
c = 3.0e8
k = 1.38e-23
omega = 5.67e-8

#The following defines the planck function
def planck(wave, T):
    """
    Input: Wavelength values
    Returns: B_nu(wave,temp), the Planck distribution function in SI units
    """
    numer = 2.*h*(c)**2/(wave*1.e-6)**5
    denom = exp(h*c/(k*T*(wave*1.e-6)))-1.
    return numer/denom

pl = planck(M2Vw,T)
pln = r*planck(M2Vw,T)/median(planck(M2Vw,T))
ext = 0.302*Av/M2Vw**2

#our star with the added planck function, normalized
star = M2fn + pln
redd = star/10**(ext/2.5)
reddn = redd/median(redd)

#Now we plot:import numpy as np
from pylab import *
from pylab import loadtxt
from scipy import *
import matplotlib.pyplot as plt
import scipy.constants as const
import scipy.stats as stats
from scipy import optimize

#Let's plot TTau in microns and
#First column is wavelength, second is flux density, third is something we ignore.
TTmic = loadtxt('/home/kmiskove/Documents/TTau_merged_microns.txt', usecols=(0,1))
TTw, TTf = (TTmic[:,0],TTmic[:,1])
M2V = loadtxt('/home/kmiskove/Documents/M2V_Gl806.txt', skiprows=45, usecols=(0,1))
M2Vw, M2Vf = (M2V[:,0],M2V[:,1])

#Let's normalize the fluxes:
TTfn = TTf/median(TTf)
M2fn = M2Vf/median(M2Vf)

#This is where we specify our Av, T, and r values.
Av = 0.3 #veiling
T = 1500 #temp
r = 0.7 #scaling
#And let's define some constants:
h = 6.26e-34
c = 3.0e8
k = 1.38e-23
omega = 5.67e-8

#The following defines the planck function
def planck(wave, T):
    """
    Input: Wavelength values
    Returns: B_nu(wave,temp), the Planck distribution function in SI units
    """
    numer = 2.*h*(c)**2/(wave*1.e-6)**5
    denom = exp(h*c/(k*T*(wave*1.e-6)))-1.
    return numer/denom

pl = planck(M2Vw,T)
pln = r*planck(M2Vw,T)/median(planck(M2Vw,T))
ext = 0.302*Av/M2Vw**2

#our star with the added planck function, normalized
star = M2fn + pln
redd = star/10**(ext/2.5)
reddn = redd/median(redd)

#Now we plot:
plt.plot(TTw,TTfn)
plt.plot(M2Vw, reddn)
plt.plot(M2Vw, pln)
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux density')
plt.show(_)

#Calculating velocity or smth...c = 3.0e8
yeet = 0.302*.3/2.1660**2
print(yeet)
#^ This is for lambda sub b

#Now we do our extinction of 0.0193 = -2.5log(F) and solve for F. I'll just do that in my calculator.
flux = 0.982381
#So we're losing 0.0176 of our light. So we take our flux from IRAF, which was 3.821E-13, and multiply it by 1.0176
#Our new, corrected flux is 3.888E-13! erg/s/cm^2
#When we use figure 3, our  mass accretion rate ends up being E-7, so we'd need E7 years to form the sun. That's 10,000,000 years. Ten million years. That sounds fine.

#Now for part 3, we measure our lamba sub b to be 10820 and our lambda not to be 10834.
lb = 10820
ln = 10834
c = 3.0e8

v = c*((ln - lb)/ln)
print(v)

#So that's 387km/s. The escape velocity of the sun from Earth's distance is about 42km/s, so we can assume this material would escape and be ejected.
plt.plot(TTw,TTfn)
plt.plot(M2Vw, reddn)
plt.plot(M2Vw, pln)
plt.xlabel('Wavelength (microns)')
plt.ylabel('Flux density')
plt.show(_)

#Calculating velocity or smth...c = 3.0e8
yeet = 0.302*.3/2.1660**2
print(yeet)
#^ This is for lambda sub b

#Now we do our extinction of 0.0193 = -2.5log(F) and solve for F. I'll just do that in my calculator.
flux = 0.982381
#So we're losing 0.0176 of our light. So we take our flux from IRAF, which was 3.821E-13, and multiply it by 1.0176
#Our new, corrected flux is 3.888E-13! erg/s/cm^2
#When we use figure 3, our  mass accretion rate ends up being E-7, so we'd need E7 years to form the sun. That's 10,000,000 years. Ten million years. That sounds fine.

#Now for part 3, we measure our lamba sub b to be 10820 and our lambda not to be 10834.
lb = 10820
ln = 10834
c = 3.0e8

v = c*((ln - lb)/ln)
print(v)

#So that's 387km/s. The escape velocity of the sun from Earth's distance is about 42km/s, so we can assume this material would escape and be ejected.