Huang, problem 9.3.  Relativisitic gas.

Problem a.

Clear[i]
$Assumptions = m > 0 && e > 0 ;i[x_, m_, e_] = Integrate[x Sqrt[x^2 - m^2], {x, m, e}]
i2[x_, m_, e_] = Integrate[x Sqrt[x^2], {x, 0, e}]
Integrate[x Sqrt[x^2 - m^2], x]

Problem e.

Clear[ef, d]
ef[c_, h_, rho_, m_] = Sqrt[(3 (c h)^3 rho/(2 Pi)) + (m c^2)^2]
d = ef[3 10^8, 6.62606957 10^(-34), 1/(0.01)^3, 0]
d 6.24150934 10^(18)

Problem b.  volume element change of vars

Clear[p, px, py, pz, v, npSq, j, jm]
$Assumptions = 0 ;
p = {px, py, pz}
npSq = px^2 + py^2 + pz^2 ;
(*Norm[p]^2*)
v = p/Sqrt[ m^2 + npSq/c^2]
jm = D[v, {p}] ;
jm // TraditionalForm
j = Det[ jm] ;
j // TraditionalForm
j// FullSimplify

problem b.  attempt the integrals

i[x_, beta_, mu_, e0_] = (x^2 - 1)^(3/2)/(x^5 (E^(beta e0 x - beta mu) + 1)) ;
i2[x_, beta_, mu_, e0_] = (x^2 - 1)^(1/2)/(x^4 (E^(beta e0 x - beta mu) + 1)) ;
i[x, β, μ, (Subscript[e, 0])] // TraditionalForm
i2[x, β, μ, (Subscript[e, 0])] // TraditionalForm
$Assumptions = β > 0 && μ < 1 &&  (Subscript[e, 0]) > 0 ;
Integrate[ i[x, β, μ, (Subscript[e, 0])], {x, 1, Infinity}]
Integrate[ i2[x, β, μ, (Subscript[e, 0])], {x, 1, Infinity}]

Problem b.  second round of the integrals.

Clear[f1, f2]
f1[a_, b_, y_] = ((y+1)^2 -1)^(3/2)/(a E^(b y) + 1) ;
f2[a_, b_, y_] = (y + 1)((y+1)^2 -1)^(1/2)/(a E^(b y) + 1) ;
f1[E^(β (Subscript[ϵ, 0]))/z, β (Subscript[ϵ, 0]), y] // TraditionalForm
f2[E^(β (Subscript[ϵ, 0]))/z, β (Subscript[ϵ, 0]), y] // TraditionalForm
$Assumptions  = β > 0 && Log[z] < β && (Subscript[ϵ, 0]) > 0 ;
Integrate[f1[E^(β (Subscript[ϵ, 0]))/z, β (Subscript[ϵ, 0]), y], {y, 0, Infinity}] 
Integrate[f2[E^(β (Subscript[ϵ, 0]))/z, β (Subscript[ϵ, 0]), y], {y, 0, Infinity}]

problem d

$Assumptions = efe > 1 ;
i1 = Integrate[ (x^2 - 1)^(3/2), {x, 1, efe}]
i2 = Integrate[ (x-1)x(x^2 - 1)^(1/2), {x, 1, efe}]
i1/i2 /. efe -> Sqrt[ (n/(Subscript[n, 0 ]))^(2/3) + 1] // FullSimplify

problem d.  attempt to sum the series in the expansion of n(zbar)

Clear[k, theta]
$Assumptions = k ∈ Integers && theta > 0 && k ≥ 0 ;

(* doesn't converge? *)
Sum[ (1/(k+1))^(k + s + 3/2) Gamma[ k + s + 3/2] (theta/2)^s Gamma[3/2]/(Gamma[s+1] Gamma[3/2-s])(3/2 + s)/(3/2 -s), {s, 0, Infinity}]

Pressure integrals for problem d.

Clear[theta, s, j, k, m, t]
$Assumptions = theta > 0 && s ∈ Integers && m ∈ Integers && s ≥ 0 && m ≥ 0 ;
(*i = Integrate[ x^(m + 1/2) (1 + theta x) Sqrt[1 + theta x/2] Exp[-x(1 + s)], 
{x, 0, Infinity}]*)
j[theta_, s_] = Integrate[  ((theta x + 1)^2-1)^(3/2) Exp[-x(1 + s)], 
{x, 0, Infinity}] ;

j[ θ, s] // TraditionalForm
$Assumptions = θ > 0 && s ∈ Integers && s ≥ 0  ;
Series[j[ θ, s], {θ, 0, 5}] // FullSimplify// TraditionalForm
Clear[t]
t = Plot[{BesselK[2, x], BesselK[2, 1/x]}, {x, 0.3, 3},
AxesLabel->{z} ,
(*, PlotRange -> Full*) 
PlotLegends->Placed[{TraditionalForm[BesselK[2, z]], TraditionalForm[BesselK[2, 1/z]]},{Right, Bottom}]
]

((3 (Sqrt[(π/2)]) (θ^(3/2)))/(((s+1))^(5/2)))+((45 (Sqrt[(π/2)]) (θ^(5/2)))/(8 (((s+1))^(7/2))))+((315 (Sqrt[(π/2)]) (θ^(7/2)))/(128 (((s+1))^(9/2))))-((945 (Sqrt[(π/2)]) (θ^(9/2)))/(1024 (((s+1))^(11/2))))+((31185 (Sqrt[(π/2)]) (θ^(11/2)))/(32768 (((s+1))^(13/2))))+O((θ^6))

problem d.  pressure Bessel exponential plot

Clear[f, theta, s, z, n, θ]
f[theta_, s_, z_] = theta(-1)^s z^(s+1)j[ theta, s] 

f[θ,s, z] // TraditionalForm
fn[theta_, n_, z_] := Sum[ f[theta, s, z], {s, 0, n}]

Manipulate[
Plot[
f[tau, s, z],
{tau, 0.1, 10},
PlotRange -> {-4, 4}
]
,{s,0,10,1}
,{{z, 0.1},0,5}
]

Manipulate[
Plot[
fn[tau, n, z],
{tau, 0.01, 3},
PlotRange -> {0, 4}
]
,{n,0,10,1}
,{{z, 0.02},0,5}
]

Clear[q]
q = Module[{w},
w = {10} ;
Plot[
Evaluate[Map[ fn[theta, #, 1] &, w]],
{theta, 0.1, 100},
AxesLabel -> { θ}
(*PlotRange -> {0, 1.5},*)
]
]
Manipulate[

Plot3D[
fn[theta, 10, z],
{theta, 0, r},
{z,0,r},
AxesLabel->{θ, z}
(*PlotRange -> {-4, 4}*)
],
{r, 5, 100}
]
Clear[r]
r = Plot3D[fn[theta,10,z],{theta,0,5},{z,0,5},AxesLabel->{θ,z}]

Clear[p]
p = Module[{w},
w = {1, 2, 3, 4} ;
Plot[
Evaluate[Map[ f[theta, #, 1] &, w]],
{theta, 0.1, 100},
AxesLabel -> { θ},
(*PlotRange -> {0, 1.5},*)
PlotLegends->Placed[Map[ Text["s ="] Text[#]&, w],{Left, Top}]]
]

Clear[g]
$Assumptions = θ > 0 ;
g[θ_, z_] =Series[ fn[θ, 4, z], {θ, 0, 4}]/Sqrt[ Pi] // FullSimplify
g[θ, z] // TraditionalForm

Clear[h, h2]
(*h[theta_, z_] = Integrate[ ((theta x + 1)^2 - 1)^(3/2)/(z E^(-x) + 1), {x, 0, Infinity}]*)
h2[theta_, z_] = Integrate[ (1 + theta x)((theta x + 1)^2 - 1)^(1/2)/( (1/z)E^(x) + 1), {x, 0, Infinity}]

number density calculations

Clear[nbinomial]
(*$Assumptions = m ∈ Integers  ;*)

nbinomial[b_, k_] = (-1)^k Binomial[-b + k, -b] (-b)/(k-b)
n[theta_, z_, n_] = Sum[ nbinomial[-2, s] (z E^(1/theta))^(s+1)/((s+1)^2) BesselK[2, (s+1)/theta], {s, 0, n}]
$Assumptions = θ > 0 ;
Clear[s]
s = Series[ n[θ, z, 2], {θ, 0, 3}]/ Sqrt[ Pi] // FullSimplify ;
s // TraditionalForm

energy integrals

Clear[theta, s, u]
$Assumptions = theta > 0 && s ∈ Integers && m ∈ Integers && s ≥ 0 && m ≥ 0 ;
(*i = Integrate[ x^(m + 1/2) (1 + theta x) Sqrt[1 + theta x/2] Exp[-x(1 + s)], 
{x, 0, Infinity}]*)
u[theta_, s_] = Integrate[  x ((theta x + 1)^2-1)^(3/2) Exp[-x(1 + s)], 
{x, 0, Infinity}] ;
u[θ, s] // TraditionalForm ;
TraditionalForm[ HoldForm[ Integrate[  x ((θ x + 1)^2-1)^(3/2) Exp[-x(1 + s)] , {x , 0, Infinity}] ]== u[θ, s]]

u[θ, s]

Manipulate[
Plot[
theta^2 nbinomial[-2, s] u[theta, s], 
{theta, 0, 4},
PlotRange -> {-2000, 2000}]
,{s,0,10, 1}]

Manipulate[
Plot[theta^5 HypergeometricU[-(3/2),-4,((2 (1+s))/theta)],
{theta, 0, 4}(*,
PlotRange -> {-2000, 2000}*)
]
,{s,0,10, 1}]

Clear[ph]
ph = Module[{s},
s = 0;
 Plot[
{ HypergeometricU[-(3/2),-4,((2 (1+s))/theta)], theta^5 HypergeometricU[-(3/2),-4,((2 (1+s))/theta)]},
{theta, 0, 2}(*,
PlotLegends->Placed[
{TraditionalForm[ HypergeometricU[-(3/2),-4,((2 (1+s))/θ)]], θ^5 TraditionalForm[HypergeometricU[-(3/2),-4,((2 (1+s))/θ)]]},{Left, Top}]*)
]]

$Assumptions = θ > 0 ;
Series[HypergeometricU[-(3/2),-4,((2 (1+s))/θ)], {θ, 0, 3}] // TraditionalForm

Clear[theta, s, uu, suu]
$Assumptions = theta > 0 && s ∈ Integers && n ∈ Integers && s ≥ 0 && n ≥ 0 ;
(*i = Integrate[ x^(m + 1/2) (1 + theta x) Sqrt[1 + theta x/2] Exp[-x(1 + s)], 
{x, 0, Infinity}]*)
uu[theta_, z_,  n_] = Sum[ nbinomial[-2, s] theta^2 z^(s+1)Integrate[  ((theta x + 1)^2-1)^(3/2) Exp[-x(1 + s)](x-1 - z E^(-x)), 
{x, 0, Infinity}] , {s, 0, n}] ;
uu[θ, z, n] // FullSimplify // TraditionalForm 
$Assumptions = θ > 0 ;
suu= Series[ uu[θ, z, 3], {θ, 0, 3}, {z, 0, 3}] ;
suu// TraditionalForm

(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig3.pdf",p]*)
(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig4.pdf",q]*)
(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig5.pdf",r]*)
(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig2.pdf",t]*)

(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig6sv.eps",ph]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig3sv.eps",p]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig4sv.eps",q]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig5sv.eps",r]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig2sv.eps",t]*)

(*Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig6sv.pdf",ph]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig3sv.pdf",p]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig4sv.pdf",q]

Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig2sv.pdf",t]

Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\figures\\huang93Fig6sv.png",ph]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\figures\\huang93Fig3sv.png",p]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\figures\\huang93Fig4sv.png",q]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\figures\\huang93Fig5sv.png",r]
Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\figures\\huang93Fig2sv.png",t]*)
(*ClearAll[peeters`]*)
<<peeters`

peeters`setGitDir::usage

peeters`exportForLatex[ "huang93Fig6sv",ph ]
peeters`exportForLatex[ "huang93Fig3sv",p ]
peeters`exportForLatex[ "huang93Fig4sv",q ]
peeters`exportForLatex[ "huang93Fig5sv",r ]
peeters`exportForLatex[ "huang93Fig2sv",t ]

(*peeters`setGitDir[\!\(\*"\""\)\!\(\*"blogit"\)\!\(\*"\""\)]*)

SetDirectory["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit"]

Export["C:\\Users\\Peeter\\cygwin_home\\physicsplay\\notes\\blogit\\huang93Fig5pd.pdf",r]

Integrals for high density comparisions

Clear[pHighDensity, uHighDensity, y, θ]
pHighDensity = Integrate[(u^2 - 1)^(3/2), {u, 1, θ y + 1}] // FullSimplify
uHighDensity = Integrate[(u^2 - u)(u^2 - 1)^(1/2), {u, 1, θ y + 1}] // FullSimplify

$Assumptions = y θ > 0 ;
pHighDensity // TraditionalForm

$Assumptions = y θ > 0 ;
uHighDensity // TraditionalForm