How Can I Plot Energy Level Occupation Numbers in Mathematica?

[Munin]

So, I've managed to get the distribution in a decent way. Using this code;

hw = 1;
kt = 25;
n = 10000;
dist[b_] := 1/(b*Exp[hw*m/kt] - 1);
normsum[b_] := Sum[dist[b]*(m + 2)*(m + 1)/2, {m, 0, 300}]
q = FindRoot[normsum[b] == n, {b, 0.5}]
occnumber = Table[N[dist[b /. q]*(m + 2)*(m + 1)/2], {m, 0, sumlim}];
Part[occnumber, 1];
ListPlot[occnumber, PlotRange -> {0, 500}]

Now I want to plot the Energy level occupation numbers where

E=Sum[m*hw*(m+1)*(m+2)*0.5/b*Exp[m*hw/kt]-1

But, I've been unable to do so, and then of course I want to get the Heat Capacity as a function of the temperature, but I think I might be able to manage that if I get some help on this.

I've been stuck for hours, please help me!

 

[Bill Simpson]

In your code above you didn't say what value sumlim should have. I'm guessing perhaps 300 to get a reasonable plot.

Are you really saying you want to assign the value of a Sum to E==2.71828?
Or do you want to assign the value of a sum to a variable of your own?

What are you trying to Sum over? I can't see if your problem is with the Sum or the Plot until you give more information. Your Sum appears to be incomplete.

 

[Munin]

Yea, the Sum should be 300. It would be easier if I just showed you what I'm trying to do:
The thing I'm trying to do is explained in:
http://www.eg.bucknell.edu/physics/ligare/ajp_bose1.pdf"

I want to get the Heat Capacity-Plot.

I've done it in Java, now all I want to do it in Mathematica as well as I'm trying to learn the language.

 

[Bill Simpson]

If you have an extraordinarily simple and well done Java implementation that looks like and is organized just like the mathematics in the paper and you wanted to show your code then it should be a pretty simple task of replacing that line by line with Mathematica.

Often it is possible to reorganize the result and replace several lines or even whole sections with a single Mathematica function, but that may be optional.