<jupyter_start><jupyter_text>03_StellarDynamics_HW_Chap06   Question (1) Low Surface Density Part<jupyter_code># initialize the symbols
kappa, kappa_0, P, T, F_rad, F_conv, a, b, mu = sp.symbols(
    'kappa kappa_0 P T F_rad F_conv a b mu')
# kappa for opacity
exprK = sp.Eq(kappa, kappa_0 * P**a * T**b)
exprK
# eq (6.121), pressure for photosphere with optical depth = 3/2
eqPh = sp.Eq(sp.simplify('P_ph'), 
             sp.simplify('2/3 * (1+a) * g / (kappa)'))
eqPh
# approx expression for F_rad in (6.123)
exprFr = sp.Eq(F_rad, sp.simplify('a * c / 4 * T^4'))
exprFr
# (6.124)
exprFc = sp.Eq(F_conv, sp.simplify('1/2 * rho * c_P * v * T'))
exprFc
# approx the radial-velocity fluctuation 
v = sp.Symbol('v')
gamma = sp.Symbol('gamma')
exprV = sp.simplify('1 / 2 * ((r * k * T)/(mu * m_H))^(1/2)').subs(
    sp.sympify('r'), gamma)
sp.Eq(v, exprV)
# (6.126) approx ρcT to internal energy
exprInternal = sp.simplify('3 / 2 * rho * (k * T / mu/ m_H)') * gamma 
eqInternal = sp.Eq(sp.simplify('rho * c_P * T'), exprInternal)
eqInternal
# ideal gas law
eqIdeal = sp.Eq(P, sp.simplify('rho * k * T'))
eqIdeal = sp.solve(eqIdeal, sp.simplify('rho'))[0]
sp.Eq(sp.simplify('rho'), eqIdeal)
# substitute (6.124)
exprFc = exprFc.subs(
    eqInternal.lhs, eqInternal.rhs).subs(
    v, exprV).subs(sp.simplify('rho'), eqIdeal)
exprFc
# This is eq (6.127)
# for energy conservation, the radiative flux which enters the
# sphere shold equals to the convective flux leaving the 
# stellar interior
sp.Eq(exprFc.lhs, exprFr.lhs)
# So we have
eqFcFr = sp.Eq(exprFc.rhs, exprFr.rhs)
eqP = sp.Eq(sp.simplify('P_ph'), sp.solve(eqFcFr, P)[0])
eqP
# So, the effective temperature for the fully convective star
# in low surface density is
## here we substitute kappa and P_ph
eqP = eqP.subs(eqPh.lhs, eqPh.rhs).subs(exprK.lhs, exprK.rhs).subs(P, eqP.rhs)

## here we symbolically solve for T
sp.solve(eqP, T)<jupyter_output><empty_output><jupyter_text>The above solution is awful.  But we can roughly see that the exponents of, e.g., $\kappa_0$ is the same as `eq (6.129)`.  There is come inconsistent between above solution and `eq (6.129)`. However, since `eq. (6.129)` is just a approximation, let's just take the approximation made by the textbook.<jupyter_code># eq (6.129), effective temperature T for low surface density
eqTeff = sp.Eq(sp.simplify('T_eff'), 
    sp.simplify('(M / mu^((a+1)/2) / kappa_0 / R^2)^(1/(b+3.5*(1+a)))'))
eqTeff<jupyter_output><empty_output><jupyter_text>- Population 1<jupyter_code># Population I effective temperature 
eqTeff.subs(a, 0.7).subs(b, 5.3).subs(kappa_0, 6.9*10**-26)<jupyter_output><empty_output><jupyter_text>- Population 2<jupyter_code># Population II effective temperature 
eqTeff.subs(a, 0.6).subs(b, 9.4).subs(kappa_0, 6.1*10**-40)<jupyter_output><empty_output><jupyter_text>Question (2) Low Surface Density Part<jupyter_code>eqL = sp.Eq(sp.simplify('L'), 
            sp.simplify('4 * pi * sigma * R^2 * T_eff^4'))
# equaiton of luminosity 
eqL = sp.Eq(sp.simplify('L'), 
            sp.simplify('4 * pi * sigma * R^2 * T_eff^4'))
eqL
# substitute for effective temperature
eqL = eqL.subs(eqTeff.lhs, eqTeff.rhs)
eqL<jupyter_output><empty_output><jupyter_text>- Population 1<jupyter_code># Population I luminosity
eqL.subs(a, 0.7).subs(b, 5.3).subs(kappa_0, 6.9*10**-26)<jupyter_output><empty_output><jupyter_text>- Population 2<jupyter_code># Population II luminosity
eqL.subs(a, 0.6).subs(b, 9.4).subs(kappa_0, 6.1*10**-40)<jupyter_output><empty_output><jupyter_text>Question (1) High Surface Density PartFor high surface density, the convective transport is efficient. The temperature gradient is adiabatic. So we may have...$$P \propto \rho^{\gamma} \propto \rho T$$$$P \propto T^{\gamma / (\gamma - 1)}$$<jupyter_code>K = sp.Symbol('K')
eqAdiabatic = sp.Eq(P, K * T**(gamma/(gamma-1)))
eqAdiabatic<jupyter_output><empty_output><jupyter_text>So, we would have the relation between $P_{ph}$ and $P_{c}$<jupyter_code>Pc, Tc = sp.symbols('P_c T_c')
P_ph, T_eff = sp.symbols('P_ph T_eff')
eqAdiabaticHigh = sp.Eq(
    eqAdiabatic.subs(P, Pc).lhs / eqAdiabatic.subs(T, Tc).rhs * K, 
    eqAdiabatic.subs(P, P_ph).lhs / eqAdiabatic.subs(T, T_eff).rhs * K
     ).subs(gamma, 5/3)
eqAdiabaticHigh
# P_ph Substitute with (6.121) 
eqHighSurface = eqAdiabaticHigh.subs(P_ph, eqPh.rhs)
eqHighSurface
# substitute kappa 
eqHighSurface = eqHighSurface.subs(kappa, exprK.rhs).subs(T, T_eff).subs(P, P_ph)
eqHighSurface
# it is awuful that I need to replace P_ph with Pc/Tc again 
eqHighSurface = eqHighSurface.subs(P_ph, 
    sp.solve(eqAdiabaticHigh, P_ph)[0])
eqHighSurface<jupyter_output><empty_output><jupyter_text>不太確定為什麼課本不給一下 core 的 profile ，但是老闕上課好像有給。$$\frac{P_c}{T_c^{2.5}} \propto M^{-0.5} R^{1.5}$$<jupyter_code># substitute core profile with above relation provided by 老闕
M, R = sp.symbols('M R')
eqHighSurface = eqHighSurface.subs(Pc, M**-0.5 * R**1.5 * Tc**2.5)
eqHighSurface
# it seems Tc should have zero exponent, 
eqHighSurface = eqHighSurface.subs(Tc, 1)
eqHighSurface
## subsitute for g = M / R^2
g = sp.Symbol('g')
eqHighSurface = eqHighSurface.subs(g, sp.simplify('M / R^2'))
eqHighSurface
## solve for T_eff
eqHighSurface = sp.Eq(T_eff, sp.solve(eqHighSurface, T_eff)[0])
eqHighSurface<jupyter_output><empty_output><jupyter_text>- Population 1<jupyter_code># Population 1 effective temperature 
eqHighSurface.subs(
    a, 0.7).subs(b, 5.3).subs(kappa_0, 6.9*10**-26)<jupyter_output><empty_output><jupyter_text>這個結果也太好笑了XD  不知道為什麼程式不把指數項爆開，但是我們可以試著幫他算一下：- exponent of M : 3.7 * 0.0523 = 0.1935, the same as (6.135)- exponent of R : 9.1 * 0.0523 = 0.4759, not the same as (6.135)R 的指數是錯的，有可能老闕給的那條式子是錯的或者我用錯地方，但是 M 的指數是對的 - Population 2<jupyter_code># Populatoin II effective temperature 
eqHighSurface.subs(
    a, 0.6).subs(b, 9.4).subs(kappa_0, 6.1*10**-40)<jupyter_output><empty_output><jupyter_text>這次是兩個都錯了。好像不是很確定哪裡算錯了。  Question (2) High Surface Density Part<jupyter_code># eq of luminosity
eqLHighSurface = sp.Eq(
    sp.simplify('L'), sp.simplify('4 * pi * sigma * R^2 * T_eff^4'))
eqLHighSurface
eqLHighSurface = eqLHighSurface.subs(T_eff, eqHighSurface.rhs)
eqLHighSurface<jupyter_output><empty_output><jupyter_text>- Population 1<jupyter_code># population I luminosity
eqLHighSurface.subs(
    a, 0.7).subs(b, 5.3).subs(kappa_0, 6.9*10**-26)
# population II luminosity
eqLHighSurface.subs(
    a, 0.6).subs(b, 9.4).subs(kappa_0, 6.1*10**-40)<jupyter_output><empty_output><jupyter_text>Question (1) Ionization Zones課本這邊寫得看不太懂，不過大概是叫我們硬解吧。Use the same adiabatic condition.<jupyter_code>eqAdiabatic<jupyter_output><empty_output><jupyter_text>The textbook says that the ionization zone is convective and the entropy across it is constant. So that $S(x=1) - S(x=0) = 0$, where x denotes the fraction of hydrogen ionized.Also, argued in the textbook, - x = 1 : the center of the star- x = 0 : at the base of the photosphereYeah ... so, basically we would have $$S(x=1, P=P_c) = S(x=0, P=P_{ph})$$<jupyter_code># entropy per unit mass
exprEntropy = sp.simplify(
    '''
    X * k / m_H * 
    (5 / 2 * (1 + x + delta) + 
    chi / k / T + 
    ln(2 * pi * m_H / h^2)^(3/2) + 
    delta * ln(2 * pi * m_He / h^2)^(3/2) + 
    x * ln(2 * pi * m_e / h^2)^(3/2) + 
    (1 + x + delta) * ln(k * T)^(5/2) * (1 + x + delta) / P)
    ''')
S = sp.Symbol('S')
epEntropy = sp.Eq(S, exprEntropy)
epEntropy
# lol this is ugly<jupyter_output><empty_output><jupyter_text>Get the relation here,  $$S(x=1, P=P_c) = S(x=0, P=P_{ph})$$<jupyter_code># S(x=1, P=P_c) = S(x=0, P=P_{ph})
x, delta = sp.symbols('x delta')
eqBoundary = sp.Eq(epEntropy.subs(x, 1).subs(P, Pc).subs(T, Tc).rhs, 
    epEntropy.subs(x, 0).subs(P, P_ph).subs(T, T_eff).rhs)

eqBoundary
# lol 超級複雜，用手算會瘋掉
# substitute delta = Y / (4 X)
eqBoundary = eqBoundary.subs(delta, sp.simplify('Y / (4 * X)'))
eqBoundary<jupyter_output><empty_output><jupyter_text>And then, substitute with $$\frac{P_c}{T_c^{2.5}} \propto M^{-0.5} R^{1.5}$$<jupyter_code>eqBoundary = eqBoundary.subs(Pc, M**-0.5 * R**1.5 * Tc**2.5)
eqBoundary<jupyter_output><empty_output><jupyter_text>天啊，感覺 $T_c$ 是消不掉的。用 `simplify` 來試試看好了。註： `simplify` 是 `sympy` 用來簡化方程式的懶人指令，他自動會幫你跑他自己覺得比較簡化的結果。<jupyter_code>eqBoundary = sp.simplify(eqBoundary)
eqBoundary<jupyter_output><empty_output><jupyter_text>都做到這裡了也只能硬著頭皮做下去了：- 用 (6.121) 替代掉 $P_{ph}$- 替代掉 $\kappa$- 再用 adiabatic equation 替代掉 多餘的 $P_{ph}$，留下$T_{eff}$ 可解<jupyter_code># 用 (6.121)
eqBoundary = eqBoundary.subs(P_ph, eqPh.rhs)
eqBoundary
# 換掉 kappa
eqBoundary = eqBoundary.subs(kappa, exprK.rhs.subs(P, P_ph).subs(T, T_eff))
eqBoundary
# adiabatic condition Pc / Tc^2.5 = P_ph / T_eff^2.5
# 但是這邊就直接把 Pc / Tc^2.5 換成  M^{-0.5} R^{1.5} 好了
eqAdiabaticHigh.subs(eqAdiabaticHigh.lhs, M**-0.5 * R**1.5)
eqBoundary = eqBoundary.subs(
    P_ph,
    sp.solve(
        eqAdiabaticHigh.subs(eqAdiabaticHigh.lhs, M**-0.5 * R**1.5),
        P_ph
        )[0]
    )
eqBoundary<jupyter_output><empty_output><jupyter_text>感覺是很嚴重的大失敗。  但是還是讓他把答案寫完好了。<jupyter_code>sp.solve(eqBoundary, T_eff)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1626205520452278, 'fraction_numerical': 0.03469681853453054, 'mean_word_length': 4.850843060959792, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '38777', 'license': 'MIT', 'n_tokens_mistral': 4286, 'n_tokens_neox': 4098, 'n_words': 814, 'path': '03_StellarDynamics_HW_Chap06.ipynb', 'repo_name': 'jibanCat/StellarDynamicsHW', 'size': '9021'}