(Subscript[H, 2]) heat capacity

B = 60.853; (* rotational constant, in (cm^(-1)) *)
k = 0.695028; (* boltzmann's constant in (cm^(-1))/K *)
(Subscript[T, o]) = 300; (* room temperature (K) *)
(Subscript[σ, o])=(B/(k (Subscript[T, o]))); (* sigma at room temperature *)
(Subscript[dσ, o])=.01; (* small displacement of sigma for approximating derivatives *) 
(Subscript[bf, 0]) = (2j+1)(E^(-σ j(j+1))) ;(* boltzmann factor with degeneracy ignoring proton spin *)
(Subscript[bf, 1]) = (2-(((-1))^j))(Subscript[bf, 0]); (* boltzmann factor w/ deg accounting for proton spin *)
(Subscript[bf, p]) = (1/2)(1+(((-1))^j))(Subscript[bf, 1]); (* boltzmann factor w/ d for para-h2 *)
(Subscript[bf, o]) = (1/2)(1-(((-1))^j))(Subscript[bf, 1]); (* boltzmann factor w/ d for ortho-h2 *)
(Subscript[C, v])[bf_] :=(σ^2)(
Total[
(
Log[
Total[
Table[bf,{j,0,5}]
]
]/.σ->{σ+dσ,σ-dσ,σ}
)*{1,1,-2}
]/((dσ)^2)
)/.dσ->(Subscript[dσ, o]) (* returns expression for heat capacity for given sigma *)
ρ=Divide@@Total[{(Subscript[bf, o]),(Subscript[bf, p])}/.j->Range[0,5],{2}]/.σ-> (Subscript[σ, o]) (* ortho-para ratio at room temp *)
ListPlot[
Append[
Function[
x,Table[{T,x},{T,50,320}]
]/@(
Append[
(Subscript[C, v])/@{(Subscript[bf, 0]),(Subscript[bf, 1])},
((ρ (Subscript[C, v])[(Subscript[bf, o])]+(Subscript[C, v])[(Subscript[bf, p])])/(1+ρ))
]/.σ-> (B/(k T))
),
Import["/home/chris/teaching/thermo/h2heatcap/data.tsv"]
],
Joined->{True,True,True,False},
PlotLegends->{"w/out spin","w/ spin","para + ortho","experiment"},
AxesLabel->{"Temperature (K)","molar heat capacity (Subscript[c, v]) (R)"},
PlotLabel->"(Subscript[H, 2]) heat capacity",
ImageSize->Large
] (* plot the heat capacities *)
Quit[]; (* reset the variables so no funny business on rerun (not sure if this does anything) *)