#//****** Переменные ниже нужно указать самостоятельно ******//
nla = 25.;                    #сколько спинов в системе
g={0:1,2:3,4:5}     # массив ln(g(E)).  Ключем является энергия, значение - ln(g(E))
# Например, если энергия E=5 вырождена 10 раз в системе, то g[5]=ln(10)=2.302585093;
tmin=0.001; #минимальная и максимальная температура, в интервале которых строить
tmax=5;     #термодинамическое усреднение. Для модели Изинга можно оставить как есть
#//****** Переменные выше нужно указать самостоятельно ******//


iener0 = next(iter(g))

for ix in range(int(tmin*1000),int(tmax*1000)):
    x=0.001*ix
    beta=1./x
    gmax=-1000.

    for ie in g:
        if g[ie]-beta*(ie-iener0)>gmax:
            gmax=g[ie]-beta*(ie-iener0)

    a0=0; ae=0; ae2=0; aecp=0; as2=0;

    for ie in g:
        weight=math.exp(g[ie]-beta*(ie-iener0)-gmax)
        a0  += weight
        ae  += weight*ie
        ae2 += weight*ie*ie

    aecp=a

    ae  /= nla; ae2 /= nla*nla;
    ae  /= a0; ae2 /= a0; aecp /= a0;

    as2 = (math.log(a0) + gmax - beta*iener0 + aecp*beta)/nla;
    cv = beta*beta*(ae2-ae*ae)*nla;

    printf("%f %e %e %e \n" % (x,ae,cv,as2));