Demo entry 6783507

Python 3

   

Submitted by anonymous on Feb 14, 2019 at 19:52
Language: Python 3. Code size: 9.3 kB.

import numpy as np
import pylab as py
from scipy.optimize import fsolve as solve
from IPython.display import display, Math, Latex

###let's define our boatload of constants
c=299792458 
h=6.62607E-34
k=1.3806E-23
sig=5.67E-8
lamb1=np.linspace(50E-9,2000E-9, 100)
#T=5500
x=(h*c)/(lamb1*k*T)

def planck1(lamb , T):
    x=(h*c)/(lamb*k*T)
    u=(2*h*c**2)/(lamb**5)
    return u/(np.exp(x)-1)

yaxis=planck1(lamb1,5500)

def norm(data):
    return data*1E9


#def plot
py.plot(norm(lamb1),yaxis)
py.ion
py.grid('on')
py.xticks(fontsize=8)
py.ylabel('$B_\lambda$',fontsize=15)
py.xlabel('$\lambda$'  '   ' 'in nanometers',fontsize=15)
py.title("Fig 2A Plancks Function vs Wavelength")
py.savefig('a2Q2a.png', format='png')
py.show()



$
    \newcommand{\mqty}[1]{\left[ \begin{aligned}#1\end{aligned} \right]}
$

We want to solve 

$$\frac{dB_\lambda}{d\lambda} = 0 $$


Which after usng the product rule $f'*g+g'*f$ and $x=\frac{hc}{\lambda kT}$ 

we find that 

$$\frac{dB_\lambda}{d\lambda} = \frac{xe^x}{e^x-1}-5 =0
$$

###using the relaxation method to solve the non-linear equation 
###solving for lambdapeak
ylamb=1.0
for i in range(5):
    ylamb=(5*(np.exp(ylamb)-1))/np.exp(ylamb)
    print(ylamb)
    
lpeak=(h*c)/(k*T*4.96)

print('The peak wavelength in nanometers is:', lpeak*1e9,'nm')

###using the relaxation method to solve the non-linear equation 
###solving for lambdapeak
ylamb=1.0
for i in range(5):
    ylamb=(5*(np.exp(ylamb)-1))/np.exp(ylamb)
    print(ylamb)
    
lpeak=(h*c)/(k*T*4.96)

print('The peak wavelength in nanometers is:', lpeak*1e9,'nm')

freq1=np.linspace(0.15E14,30E14,100)

def planck2(freq,T):
    return (2*h*freq**2)/(c**2*np.exp((h*freq)/k*T))

py.plot((freq1),yaxis)
py.ion
py.grid('on')
py.xticks(fontsize=8)
py.ylabel('$B_\mu$',fontsize=15)
py.xlabel('$\mu$'  '   ' 'in Hz',fontsize=15)
py.title("Fig 2B Plancks Function vs Frequency")
py.savefig('a2Q2B.png', format='png')
py.show()


yfreq=1.0

for i in range(5):
    yfreq=(3*(np.exp(yfreq)-1))/np.exp(yfreq)
    print(yfreq)

freqpeak=(yfreq*k*T)/(h)
print('  ')

print('      ')
print("The peak frequency is:", freqpeak," Hz")
lambpeak2=c/freqpeak
print('    ')
print('The peak wavelength correesponding to this frequency is',lambpeak2*1E9,"nm")
Where we have found $\lambda_m= 928.45nm$ using $x=2.819$ as seen on the figure using $\mu=3.23E14Hz$.

Starting with 

$$\lambda = \frac{c}{\mu}$$ 

we have 

$$d\lambda = \frac{c}{\mu^2}d\mu $$

Which shows that 

$$d\lambda \neq d\mu$$

from

$$B\lambda =\frac{dE}{dAdtd\omega d\lambda} $$

and using

$$dln\lambda = \frac{d\lambda}{\lambda} $$

we have 

$$\frac{dE}{dAdtd\omega dln\lambda} = \frac{\lambda dE}{dAdtd\omega d\lambda} = \lambda B_\lambda $$

Similarly,

$$\frac{dE}{dAdtd\omega dln\mu} = \frac{\mu dE}{dAdtd\omega d\mu} = \mu B_\mu $$

Starting with 

$$B\lambda = \frac{2hc^2}{\lambda^5} \frac{1}{e^\frac{hc}{\lambda kt}-1} $$

where

$$\lambda B\lambda = \frac{2hc^2}{\lambda^4} \frac{1}{e^\frac{hc}{\lambda kt}-1} $$

we find 
$$\frac{d}{d\lambda} (\lambda B\lambda) = \frac{-4}{\lambda^5}\frac{1}{e^\frac{hc}{\lambda kt} -1  ^2} + \frac{1}{\lambda^4}\frac{e^\frac{hc}{\lambda kt}}{e^\frac{hc}{\lambda kt}-1} * \frac{hc}{\lambda^2 kt} $$

using 
$$x=\frac{hc}{\lambda kT}$$
we have
$$\frac{d}{d\lambda} (\lambda B\lambda) = \frac{-4}{\lambda^5}\frac{1}{e^x -1 } + \frac{1}{\lambda^4}\frac{e^x}{(e^x -1)^2} * \frac{x}{\lambda}$$

multiplying by
$$\frac{e^x -1}{\lambda ^5}$$

we find 

$$\frac{xe^x}{e^x -1} -4 = 0 $$

which we will solve using the relaxation method
***Same steps were taken in our previous derived non-linear equations

yfreq=1.0
for i in range(6):
    yfreq=(4*(np.exp(yfreq)-1))/np.exp(yfreq)
    print(yfreq)


lambdapeak2 = (h*c)/(k*T*yfreq)
print("       ")
print("For T=5500k and x=3.92 the peak wavelength is ", lambdapeak2*1E9,"nm")

lamb3=np.linspace(20E-9,20000E-9,10000)
lamb4=np.linspace(400E-9,800E-9,1000)

def planck3(lamb , T):
    x=(h*c)/(lamb*k*T)
    u=(2*h*c**2)/(lamb**4)
    return u/(np.exp(x)-1)

yaxis1 = np.log10(planck3(lamb3,3000))
xaxis1 = np.log10(lamb3)
yaxis2 = np.log10(planck3(lamb3,5500))
yaxis3 = np.log10(planck3(lamb3, 30000))
xaxis2= np.log10(lamb4)

py.plot(xaxis1,yaxis1)
py.plot(xaxis1,yaxis2)
py.plot(xaxis1,yaxis3)
py.xlim(-8.0,-4.5)
py.ylim(4,11)
py.xlabel('$log_{10}\lambda$' ,fontsize=13)
py.ylabel('$log_{10} \lambda B_\lambda nm$',fontsize=13)
py.title("Fig 2E.1 $\log$(Plancks Function) vs $\log$Wavelength for T=30000K,5500K,3000K")
py.savefig('a2Q2E1.png', format='png')
py.show()

py.plot(xaxis2,np.log10(planck3(lamb4,3000)))
py.plot(xaxis2,np.log10(planck3(lamb4,5500)))
py.plot(xaxis2,np.log10(planck3(lamb4,30000)))
py.xlim(-6.5,-6)
py.ylim(4,11)
py.xlabel('$log_{10}\lambda$' ,fontsize=13)
py.ylabel('$log_{10} \lambda B_\lambda nm$',fontsize=13)
py.title("Fig 2E.2 $\log$(Plancks Function) vs $\log$Wavelength for T=30000K,5500K,3000K over visible range")
py.savefig('a2Q2E2.png', format='png')
py.show()

$\lambda B_\lambda$ varies considerably between the hottest temperature star (green) (30 000K) and the coolest temperature star (blue) (3000k) as we can see from the graph. The Sun like star has a near constant value of $\lambda B_\lambda$. Since is bolometric flux doesn't vary greatly over wavelength, the sun is white in colour to our eyes (the entire visible spectrum). 

We will start with 
$$\frac{N^{II}}{N_t}=\frac{\frac{N^{II}}{N^I}}{1+\frac{N^{II}}{N^I}} $$

which can be rewritten as:

$$\frac{N^II}{N_t} (1-\frac{N^II}{N_t} = \frac{N^II}{N^I} $$

where

$$\frac{N^II}{N^I}$$

is given by the saha equation:

$$\frac{N^II}{N^I} = \frac{2Q^{II}}{2*n_e Q^I} \frac{2 \pi m_e KT}{h^2}^\frac{3}{2} exp^\frac{-\chi _1}{kT}$$

For hydrogen $Q^I = 2$ and $Q^{II}=1$, as well as $n_e=\frac{N^{II}}{V}$ as each ionized hydrogen contributes an electron. 

We also have

$$N_t=\rho v(m_p+m_e) = \frac{\rho v}{m_p}$$ 

as 

$$m_p >> m_e$$

Combining this we have

$$N_e=\frac{N^{II}\rho}{N_t m_p}$$

and subbing into the saha equation gives 

$$\frac{N^{II}}{N^I} = \frac{N_t m_p}{N^{II}\rho} \frac{2 \pi m_e KT}{h^2}^\frac{3}{2} e^\frac{-\chi _1}{kT}$$

and using our first rewritten equation we get

$$\frac{N^{II}}{N^t} [1 + \frac{N_t m_p}{N^{II}\rho} \frac{2 \pi m_e KT}{h^2}^\frac{3}{2} e^\frac{-\chi _1}{kT}]\frac{N_t m_p}{N^{II}\rho} \frac{2 \pi m_e KT}{h^2}^\frac{3}{2} e^\frac{-\chi _1}{kT} =\frac{N^{II}}{N^I}$$

Which we can rearrange and solve for zero as so:

$$\frac{N^{II}}{N_t}^2 + \frac{N{II}}{N^t}\frac{m_p}{\rho}\frac{2 \pi m_e KT}{h^2}^\frac{3}{2} e^\frac{-\chi _1}{kT} -\frac{m_p}{\rho}\frac{2 \pi m_e KT}{h^2}^\frac{3}{2} e^\frac{-\chi _1}{kT} = 0$$

def tau(dis):
    p = 1E-6
    k = 3.0
    return p*k*dis

s = np.linspace(0,333333,1000)

def T(s):
    Te = 10000
    p = 10**(-6)
    K = 3.0
    return ((3.0/4)*Te**4*(p*K*s+2.0/3))**(1.0/4)




yaxis4=(T(s))

py.plot(s/1000,yaxis4)
py.xlabel('Depth in Km')
py.ylabel('Temperature(K)$')
py.title("Fig 4A Temperature of A0V star as function of depth")
py.savefig('a2Q4A.png', format='png')
py.show()

From $$T=T_E$$ when $$\tau =\frac{2}{3}$$ we have $$s=\frac{2}{3}*\frac{1}{\rho\kappa}$$
$$s=222km$$
$$s=3.2*10^{-4} R_\odot$$

At the surface where $\tau = 0$
$$T_s=2^{\frac{-1}{4}}T_e$$

$$T_s=0.84T_e$$

###plot 4b

def F2(t):
    mp = 1.67e-27
    me = 9.11e-31
    k = 1.38e-23
    hbar = 6.626e-34/(2.*np.pi)
    C = (mp/(1*10**(-6)))*(me * k * t / (2.*np.pi*hbar**2))**(3./2) * np.exp(-13.6*1.602e-19/(k*t))
    NI2_NI1 = 4.0*np.e**(-13.6*(3.0/4)/(8.6173324*10**(-5)*t))
    NI3_NI1 = 9.0*np.e**(-13.6*(8.0/9)/(8.6173324*10**(-5)*t))
    NI4_NI1 = 16.0*np.e**(-13.6*(15.0/16)/(8.6173324*10**(-5)*t))
    NI5_NI1 = 25.0*np.e**(-13.6*(24.0/25)/(8.6173324*10**(-5)*t))
    NI6_NI1 = 36.0*np.e**(-13.6*(35.0/36)/(8.6173324*10**(-5)*t))
    NII_Nt = (-1.0*C/2.0)+(C/2.0)*(np.sqrt(1.0+4.0/C))
    NI_Nt = 1.0-NII_Nt
    N2I_NI = NI2_NI1/(1+NI2_NI1+NI3_NI1++NI4_NI1+NI5_NI1+NI6_NI1)
    return N2I_NI*NI_Nt


def tau_balm(f,s):
    p = 1*10**(-6)
    k = 3.5*10**5
    return f*p*k*s

s = np.linspace(0,3333333,1000)

f1=F2(T(s))


py.plot(s/1000,f1)
py.title("Fig 4B Fraction of hydrogen atoms in first excited state of typical A0V star vs its Depth")
py.xlabel("Depth in km")
py.ylabel("$f_2(s)$")
py.savefig('a2Q4B.png', format='png')
py.show()

### 4c ##
tau1=tau(s)
tau2=tau_balm(f1,s)



py.plot(s/1000,tau1)
py.plot(s/1000,tau2)
py.ylim(0,1)
py.xlim(0,500)
py.title("Fig 4C Optical Depth of Balmer series for A0V Star")
py.ylabel("$\ tau$",fontsize=12)
py.xlabel("Depth in km")
py.savefig('a2Q4C.png', format='png')
py.show()

At $\lambda =656.3nm$ the Balmer photons originate closer to the surface of the star where the temperature is cooler. We would expect the flux from Balmer photons to be less. We can roughly estimate the fraction of photons by taking the ratio of tempertature $\frac{T}{Te}^{1/4}$ where $T\cong 9400K$ and $T_e=10000K$
we then have
$$\frac{9400}{10000}^{1/4} \cong 0.8$$

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).