stot[m_,n_]:=Log[Binomial[m,n]]

s12[m_,n_,n1_]:=Log[Binomial[m/2,n1]]+Log[Binomial[m/2,n-n1]]

al=0
facstir[n_]:=n Log[n]-n+al Log[Sqrt[2 Pi n]]

sstir[m_,n_]:=facstir[m]-facstir[n]-facstir[m-n]
s12stir[m_,n_,n1_]:=sstir[m/2,n1]+sstir[m/2,n-n1]

Plot[s12stir[10000.,1000.,n1],{n1,0,1000},PlotRange->All]

Needs["PlotLegends`"]

Manipulate[
Show[
al=If[fstir,1,0];
Graphics[{Red,Dashed,Thick,Line[{{0,stot[m,n]},{n,stot[m,n]}}],
Dashing[{}],Blue,Line[{{0,sstir[m,n]},{n,sstir[m,n]}}],
}],
ListPlot[Table[{n1,s12[m,n,n1]},{n1,0,n}],PlotStyle->Red],
Plot[s12stir[m,n,n1],{n1,0,n},PlotStyle->Blue],
Graphics[{RGBColor[0,0,0.8],
Text["S(N,M) [Stirl.]",Offset[{5,0},{0,sstir[m,n]}],{-1,0},Background->White],
RGBColor[0.8,0,0],
Text["S(N,M) [Ex.]",Offset[{-5,0},{n,stot[m,n]}],{1,0},Background->White]
}
],
PlotRange->All,Frame->True,FrameLabel->{"(Subscript[N, 1])","S/kB"},AspectRatio->1
],
{{m,100,"M"},10,1000},{{n,40,"N"},3,m/2},
{{fstir,False,"Improved Stirling"},{True,False}},
SaveDefinitions->True
]

independent spin s

ei=FullSimplify[n Log[n]-m Log[m]-(n-m)Log[n-m]/.m->(n+ee)/2]

uu=FullSimplify[D[ei,ee]]

Solve[uu==beta,ee]

M[μ_,β_,B_]:=μ(Exp[β ϵ]-Exp[-β ϵ])/(Exp[-β ϵ]+Exp[β ϵ])/.ϵ->(μ B)

pp=Plot[Evaluate[{M[0.055,a/10,x],M[0.055,a/4.5,x],M[0.055,a/2.7,x]}/.a->1.5],{x,-50,50},PlotStyle->{RGBColor[0,0,0.8],RGBColor[0,0.5,0],RGBColor[0.8,0,0]},Frame->True,FrameLabel->{"H","M/at"}]

<<"MCAddOns/NewFrame.m"

xx=xkcdFramedGraphics[pp,FrameLabel->{"H","M/at"},PlotRange->All,ImageResolution->150]

xk=xkcdConvert[xx,ImageResolution->200]

Export["/tmp/paramagnetic-sm.png",xk,"PNG"]

Plot[(Exp[β ϵ]-Exp[-β ϵ])/(Exp[-β ϵ]+Exp[β ϵ])/.ϵ->1,{β,0,5},PlotRange->All]

lnw[n_,ee_]:=n Log[n]+(1/2) (n Log[4]+(ee-n) Log[-ee+n]-(ee+n) Log[ee+n])

beta[n_,ee_]:=(1/2) (Log[-ee+n]-Log[ee+n])

scaled=FullSimplify[lnw[n,ep n]/(n ),Assumptions->n>1]

bscaled=FullSimplify[beta[n, ep n],Assumptions->n>1]

ff=xkcdFramedGraphics[Plot[{scaled,bscaled,0},{ep,-1,1},PlotRange->{-0.5,1},PlotStyle->{Directive[Red,Thick],Directive[Blue,Dotted,Thick],Directive[Thin,Black]}],
FrameLabel->{"E/Nϵ",Row[{Style["ln(Ω)/N  ",Red],Style[" β",Blue]}]},ImageSize->4*72]

xx=xkcdConvert[ff,ImageResolution->250]

Export["/tmp/paramagnetic-lnw.png",xx,"PNG"]

Plot[{lnw[100,ee 100]},{ee,-1,1}]

Manipulate[
Show[
al=If[fstir,1,0];
Graphics[{Red,Dashed,Thick,Line[{{0,stot[m,n]},{n,stot[m,n]}}],
Dashing[{}],Blue,Line[{{0,sstir[m,n]},{n,sstir[m,n]}}],
}],
ListPlot[Table[{n1,s12[m,n,n1]},{n1,0,n}],PlotStyle->Red],
Plot[s12stir[m,n,n1],{n1,0,n},PlotStyle->Blue],
Graphics[{RGBColor[0,0,0.8],
Text["S(N,M) [Stirl.]",Offset[{5,0},{0,sstir[m,n]}],{-1,0},Background->White],
RGBColor[0.8,0,0],
Text["S(N,M) [Ex.]",Offset[{-5,0},{n,stot[m,n]}],{1,0},Background->White]
}
],
PlotRange->All,Frame->True,FrameLabel->{"(Subscript[N, 1])","S/kB"},AspectRatio->1
],
{{m,100,"N"},10,1000},{{n,40,"N"},3,m/2},
{{fstir,False,"Improved Stirling"},{True,False}},
SaveDefinitions->True
]