p.75: N_nu= 8 pi/c^3 nu^2/(exp(x)-1) where x=vu/vu*= h nu/kT, nu*=kT/h
p.74: F_nu= 2 pi/c^2 h nu^3/(exp(x)-1) 

N_lambda=8 pi  lambda^-4/(exp(x)-1) where x=lambda*/lambda, lambda*=hc/kT
F_lambda=2 pi c^2 h lambda^-5/(exp(x)-1) 

N_nu ~ x^2/(exp(x)-1)  max @ x=1.59362, average: <x>=2.70118
F_nu ~ x^3/(exp(x)-1)  max @ x=2.82144
let   y=1/x=lambda kT/hc
N_lambda ~ y^-4/(exp(1/y)-1)  max @ y=0.255057, average: <y>=0.684216, <1/y>=2.70118
F_lambda ~ y^-5/(exp(1/y)-1)  max @ y=0.201405
median x=2.35677, y=0.42431 (these must be inverses)

p. 74: F nu_max = 2.821 kT/h
p. 75: <energy>=2.701 kT
Wien: lambda_max = 0.201405 hc/kT  = 2.898 mm K/T
total energy flux=sigma T^4
---
math:
Integrate[x^3/(Exp[x]-1),{x,0,Infinity}]/Integrate[x^2/(Exp[x]-1),{x,0,Infinity}]

             4
           Pi
Out[1]= ----------
        30 Zeta[3]

In[2]:= N[%]                                                                    

Out[2]= 2.70118

In[1]:= FindRoot[D[x^2/(Exp[x]-1),x]==0,{x,1}]                                  

Out[1]= {x -> 1.59362}
* ? kb/h*1.59362
    33205641974.20406    

In[15]:= FindRoot[D[x^3/(Exp[x]-1),x]==0,{x,2}]                                 

Out[15]= {x -> 2.82144}


In[10]:=  FindRoot[D[y^-4/(Exp[1/y]-1),y]==0,{y,.3}]                            

Out[10]= {y -> 0.255057}
* ? h*c/kb*.255057
    .3669702355352369E-02

In[17]:= FindRoot[D[y^-5/(Exp[1/y]-1),y]==0,{y,.2}]                             

Out[17]= {y -> 0.201405}


Integrate[x^-3/(Exp[1/x]-1),{x,0,Infinity}]/Integrate[x^-4/(Exp[1/x]-1),
{x,0,Infinity}]                                                                 

             2
           Pi
Out[7]= ----------
        12 Zeta[3]

In[8]:= N[%]                                                                    

Out[8]= 0.684216

Integrate[x^-5/(Exp[1/x]-1),{x,0,Infinity}]/Integrate[x^-4/(Exp[1/x]-1)
,{x,0,Infinity}]                                                                

              4
            Pi
Out[12]= ----------
         30 Zeta[3]

In[13]:= N[%]                                                                   

Out[13]= 2.70118

---nu cff---
Integrate[x^2/(Exp[x]-1),x]
          3
        -x     2          x                    x                  x
Out[1]= --- + x  Log[1 - E ] + 2 x PolyLog[2, E ] - 2 PolyLog[3, E ]
         3
In[16]:= Integrate[x^2/(Exp[x]-1),{x,0,Infinity}]                               

Out[16]= 2 Zeta[3]

x->0: - 2 Zeta[3]= -2.40411; x->Infinity: 0
In[7]:= FindRoot[%1 == -2.40411/2, {x,2}]                                       

                                  -15
Out[7]= {x -> 2.35677 - 6.11274 10    I}

so median x=2.35677

total photon count:  
Integrate[x^2/(Exp[x]-1),{x,0,Infinity}]                               

Out[10]= 2 Zeta[3]
N[2 Zeta[3]]                                                           

Out[19]= 2.40411


so total number of photons is  8 pi/c^3 vu*^3 2 Zeta[3]

* ? 8*pi*(kb/h/c)^3*2.40411
    20286793.76616503    /m^3
p. 75: N=20.2 T^3 photons/cm^3
    
In[13]:= Integrate[x^3/(Exp[x]-1),{x,0,Infinity}]                               

           4
         Pi
Out[13]= ---
         15

median y from ccf
Integrate[y^-4/(Exp[1/y]-1),{y,0,Infinity}]
Out[27]= 2 Zeta[3]
FindRoot[NIntegrate[y^-4/(Exp[1/y]-1),{y,0,q}]==Zeta[3],{q,.42}]
NIntegrate::nlim: y = q is not a valid limit of integration.

Out[29]= {q -> 0.424311}

dev.new()
postscript("file.eps")
dev.off()

D=read.csv("1000_lambdaC.dat",sep="",header= FALSE,skip=1)
Dv=as.vector(unlist(D))
> mean(1/Dv)
[1] 2.770042

mean energy= <1/lambda>hc = 2.70118 kT
* ? 2.770042e6*h*c/(2.70118*kb)
    14754.56540753169    

mean(Dv)
[1] 0.644351

mean <lambda kT/hc> = 0.684216
* ? 0.684216*h*c/(kb*0.644351e-6)
    15277.92283905181    
    
median(Dv)
[1] 0.422
median y = kT lambda/hc = 0.424311
* ? h*c*0.424311/(.422e-6*kb)
    14466.56535023682    



> library("moments")
> moment(Dv,order=-1)
[1] 2.770042
> moment(Dv,order=1)
[1] 0.644351

hist(Dv)
hist(log(Dv))
plot(density(Dv))
qq=density(Dv,from=.1,to=1.2,n=15)
qq2=data.frame(qq$x,qq$y)
colnames(qq2)=c("lambda","dens")
 qq2
      lambda      dens
1  0.1000000 0.5915750
2  0.1785714 1.3078130
3  0.2571429 1.7154025
4  0.3357143 1.6294342
5  0.4142857 1.3511596
6  0.4928571 1.0627381
7  0.5714286 0.8604672
8  0.6500000 0.7120465
9  0.7285714 0.5754148
10 0.8071429 0.4491532
11 0.8857143 0.3171017
12 0.9642857 0.2293483
13 1.0428571 0.1703364
14 1.1214286 0.1304505
15 1.2000000 0.1089141

out2=nls(dens~ k1*lambda^-4/(exp(k2/lambda)-1),start=list(k1=1,k2=1.2),data=qq2)
> summary(out2)

Formula: dens ~ k1 * lambda^-4/(exp(k2/lambda) - 1)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
k1  0.42859    0.05145    8.33 1.43e-06 ***
k2  1.02875    0.03325   30.94 1.46e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1372 on 13 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 2.04e-06

lambda*=hc/kT=k2=1.02875e-6
* ? h*c/(kb*1.02875e-6)
    13985.68509188549    

library(ggplot2)
lambda=(seq(from = .05, to =1.5, length.out=256))
predf2 <-  data.frame(lambda)
predf2$dens=predict(out2, predf2)
ggplot(qq2, aes(lambda, dens))+  geom_point()+
  geom_line(data= predf2, aes(lambda, dens  ), col=2)


q=hist(log(3e14/Dv))
 str(q)
List of 6
 $ breaks  : num [1:11] 31 31.5 32 32.5 33 33.5 34 34.5 35 35.5 ...
 $ counts  : int [1:10] 2 8 26 57 83 219 273 234 83 15
 $ density : num [1:10] 0.004 0.016 0.052 0.114 0.166 0.438 0.546 0.468 0.166 0.03
 $ mids    : num [1:10] 31.2 31.8 32.2 32.8 33.2 ...
 $ xname   : chr "log(3e+14/Dv)"
 $ equidist: logi TRUE
 - attr(*, "class")= chr "histogram"

qq=data.frame(exp(q$mids),q$counts/(exp(q$mids)*sinh(.25)))
colnames(qq)=c("freq","dens")
out=nls(dens~ k1*freq^2/(exp(freq/k2)-1),start=list(k1=1e-30,k2=5e14),data=qq)

> summary(out)

Formula: dens ~ k1 * freq^2/(exp(freq/k2) - 1)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
k1 3.210e-41  4.864e-42   6.599 0.000169 ***
k2 2.916e+14  2.056e+13  14.184 5.94e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.774e-13 on 8 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 5.231e-06

nu*=kT/h
* ? 2.916e+14*h/kb
    13994.59743500830    

freq=exp(seq(from = min(log(qq$freq))-1, to = max(log(qq$freq)), length.out=256))
predf <-  data.frame(freq)
predf$dens=predict(out, predf)
ggplot(qq, aes(freq, dens))+  geom_point()+
  geom_line(data= predf, aes(freq, dens  ), col=2)
  

  math:
data=Import["1000_lambdaC.dat","Data",HeaderLines -> 1]
Dv=Flatten[data]

Mean[Dv]
Out[33]= 0.644351

Median[Dv]
Out[34]= 0.422

Dvi=1/Dv;
Mean[Dvi]
Out[36]= 2.77004

Moment[Dv,-1]
Out[37]= 2.77004

Dh=HistogramDistribution[Dv]
Plot[CDF[Dh,x],{x,0,10},PlotRange->Full] 
Plot[PDF[Dh,x],{x,0,1},PlotRange->Full]

Ds=SmoothKernelDistribution[Dv]
Plot[PDF[Ds,x],{x,0,1},PlotRange->Full]


In[34]:= FindMaximum[PDF[Ds,x],{x,.3}]                                          

FindMaximum::lstol: 
   The line search decreased the step size to within the tolerance specified by
     AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase
     in the function. You may need more than MachinePrecision
     digits of working precision to meet these tolerances.

Out[34]= {1.72553, {x -> 0.279527}}
max y = 0.255057 = lambda* kT/hc
* ? 0.255057*h*c/(0.279527e-6*kb)
    13128.25721791587    


sdata=Table[{x,PDF[Ds,x]},{x,.05,1,.05}]
out=NonlinearModelFit[sdata,k1 x^-4 /(Exp[k2/x]-1),{{k1,1},{k2,1.2}},{x}]
In[45]:= out["ParameterConfidenceIntervalTable"]                                
                                                                                
Out[45]=      Estimate   Standard Error   Confidence Interval

                                          0.339237
         k1   0.429263   0.0428508        0.519289

                                          0.971025
         k2   1.02887    0.0275313        1.08671

Show[ListPlot[sdata], Plot[out[x], {x, .01, 1}]]



fit/plot:use spreadsheet by columns to make file  364t104_1.fit
vanvleck 2002$ fitxy
* set file  364t104_1.fit xc 1 xec 2 yc 3 yec 4
* re
 lambda 	Dlambda			prob			Dprob
 READ Done.
 
* set f(x)=k1*x^-4/(exp(k2/x)-1) k1 .2 k2 1
* re
* fit
Enter list of Ks to vary, e.g. K1-K3,K5 k1-k2
 chi-square=   217.8549    
 
   7 iterations used
 REDUCED chi-squared=  0.4685314      chi-squared=   3.748251    
 K1=   .4175241         K2=   1.005667    
Display covariance matrice? No, Screen, File [N,S,F] 

* ? h*c/(1.005666e-6*kb)
    14306.71170972987    
vanvleck 1998$ plot
set file  364t104_1.fit xc 1 xec 2 yc 3 yec 4
re
 lambda         Dlambda                 prob                    Dprob
 READ Done.
 
set f(x)=k1*x^-4/(exp(k2/x)-1)
set K1=   .4175241         K2=   1.005667
sc
bo
dp
fc
set pmin .01
cl

bo
dp
fc


---
Wed Sep 19 15:18:43 CDT 2018
D=read.csv("1000_lambdaD.txt",sep="",header= FALSE,skip=1)
Dv=as.vector(unlist(D))
> mean(1/Dv)
 0.940181
median(Dv)
[1] 1.22095
mean(Dv)
[1] 1.967681

mode: x=1.59362
median: x=2.35676
average: x=2.70118

* ? h*c/1.22095e-6/kb/2.35676
    5000.119176344443    
* ? h*c*0.940181e6/2.70118/kb
    5007.852609967125    

plot(density(Dv))
qq=density(Dv,from=.2,to=4,n=15)
qq2=data.frame(qq$x,qq$y)
colnames(qq2)=c("lambda","dens")

out2=nls(dens~ k1*lambda^-4/(exp(k2/lambda)-1),start=list(k1=1,k2=.5),data=qq2)
Parameters:
   Estimate Std. Error t value Pr(>|t|)    
k1 11.11416    1.01335   10.97 6.10e-08 ***
k2  3.02744    0.07367   41.10 3.78e-15 ***

* ? h*c/(kb*3.02744E-6)
    4752.455387481570    
    
lambda=(seq(from = min((qq2$lambda))-.1, to = max((qq2$lambda))+.5, length.out=256))
predf <-  data.frame(lambda)
predf$dens=predict(out2, predf)
dev.new()
pdf(file="1000D_R.pdf", width=7,height=7)
ggplot(qq2, aes(lambda, dens))+  geom_point()+
  geom_line(data= predf, aes(lambda, dens  ), col=2)
dev.off()
  
fitxy
set file 1000D_density.dat yc 1 yec 2 xc 3 xec 4
re

set f(x)=k1*x^-4/(exp(k2/x)-1) k1 11000 k2  3.02744  
 fit
Enter list of Ks to vary, e.g. K1-K3,K5 k1
   7 iterations used
 REDUCED chi-squared=  4.2832617E-02  chi-squared=  0.4283262    
 K1=   11512.06    

* fit
Enter list of Ks to vary, e.g. K1-K3,K5 k1,k2
 chi-square=   4.981912    
 
   4 iterations used
 REDUCED chi-squared=  1.5962031E-03  chi-squared=  1.4365828E-02
 K1=   10293.15         K2=   2.916973    

 COVARIANCE MATRIX:
  0.36E+07  0.32E+03
  0.32E+03  0.31E-01
                    BETA      LOWER     UPPER      S.D.  ___ 95% Confidence ___
                                                    BETA         Interval

       1   1.02931489E+04  -1.80+308  1.80+308  7.59E+01   1.01E+04 to  1.05E+04
       2   2.91697314E+00  -1.80+308  1.80+308  7.02E-03   2.90E+00 to  2.93E+00