tic
temperatura=300;  % u K
vibraciona_frek=300;  % u cm-1
rotaciona_konstanta=10;  % u cm-1
broj_nivoa=20;
nivoi=0:broj_nivoa;

teta_vibr=1.44*vibraciona_frek;
teta_rot=1.44*rotaciona_konstanta;

vibraciona=exp(-nivoi*teta_vibr/temperatura);
vibraciona=vibraciona/sum(vibraciona);
rotaciona=(2*nivoi+1).*exp(-nivoi.*(nivoi+1)*teta_rot/temperatura);
rotaciona=rotaciona/sum(rotaciona);

figure
plot(nivoi,vibraciona)
xlabel('Broj nivoa');
ylabel('Verovatnoca');
legend ('Bolcmanova raspodela za vibracione nivoe')
figure
plot(nivoi,rotaciona)
xlabel('Broj nivoa');
ylabel('Verovatnoca');
legend ('Bolcmanova raspodela za rotacione nivoe')
toc