Initialization and parameters

Import packages and define directories

these are custom packages that should be present in ~.Mathematica/Applications/

<<"MCAddOns/NewFrame.m"
<<"MCAddOns/FixPolygons.m"
Needs["ErrorBarPlots`"]

DROOT=NotebookDirectory[]
DDIR=DROOT<>"data/"
FDIR=DROOT<>"figures/"

Define plot style options

size of graphics for output (journal column size)

onecol=3.4*72
twocol=6.7*72

style for text

tstyle={FontFamily->"Helvetica",FontSize->10,SingleLetterItalics->False};
txtFormat[text_,style_:tstyle]:=Style[text,style]

style for drawing frames

framestyle=Sequence[LabelOffset->-3,PlotRangeOffset->2,FrameStyle->Directive[AbsoluteThickness[0.5],Black],LabelStyle->tstyle];

style for lines in plots

linestyles[xopts_:{}]:={
Directive[Flatten[{RGBColor[0.8,0,0],AbsoluteThickness[1],xopts}]],
Directive[Flatten[{RGBColor[0.1,0.2,0.9],AbsoluteThickness[1],xopts}]],
Directive[Flatten[{RGBColor[0.0,0.5,0.1],AbsoluteThickness[1],xopts}]],
Directive[Flatten[{RGBColor[0.5,0.5,0.5],AbsoluteThickness[1],xopts}]],
Directive[Flatten[{Black,AbsoluteThickness[1],xopts}]]
}
linedashing=Directive[AbsoluteDashing[{2,4,2,4,2,4}]];
fillstyles[xopts_:{},mix_:0.75]:={
Directive[Flatten[{Blend[{RGBColor[0.8,0,0],White},mix],xopts}]],
Directive[Flatten[{Blend[{RGBColor[0.1,0.2,0.9],White},mix],xopts}]],
Directive[Flatten[{Blend[{RGBColor[0.0,0.5,0.1],White},mix],xopts}]],
Directive[Flatten[{Blend[{RGBColor[0.5,0.5,0.5],White},mix],xopts}]]
}

style for color scales

Module[{f},
cbarHot[alpha_,range_:{0,1},opfun_:((1.0)&)]:=Directive[f=Min[1,Max[0,(alpha-range[[1]])/(range[[2]]-range[[1]])]];
Blend[{
Blend[{ColorData["DeepSeaColors"][(1-f)^1.1],White},(1-f)^4],
RGBColor[1.0,0.5,0]},(f)^10
]
,Opacity[opfun[(alpha-range[[1]])/(range[[2]]-range[[1]])]]
]
];

Module[{f},
cbarBWR[alpha_,range_:{0,1},opfun_:((1.0)&)]:=Directive[f=Min[1,Max[0,(alpha-range[[1]])/(range[[2]]-range[[1]])]];
Blend[{RGBColor[1,0.3,0],RGBColor[0,0.4,1],White,Black},{Sin[f*Pi/2],Cos[f*Pi/2],Sin[f*Pi]^4*2+Sin[f*Pi]^16*4+Sin[f*Pi]^64*8+Sin[f*Pi]^256*32,Cos[f*Pi]^16}] ,Opacity[opfun[(alpha-range[[1]])/(range[[2]]-range[[1]])]]
]
];

ContourPlot[x,{x,0,1},{y,0,1},ColorFunction->cbarHot,Contours->40]

periodic color scale

a2p[r_,a_,o_:{0,0}]:=o+r*{Cos[a],Sin[a]}
wedge[r_,dr_,a_,da_,n_:10,o_:{0,0}]:=
Polygon[Join[{a2p[r-dr/2,a-da/2,o],a2p[r+dr/2,a-da/2,o]},
Table[a2p[r+dr/2,a-da/2+da*i/n,o],{i,0,n}],
{a2p[r+dr/2,a+da/2,o],a2p[r-dr/2,a+da/2,o]},
Table[a2p[r-dr/2,a+da/2-da*i/n,o],{i,0,n}]
]]
wheel[o_,r_,dr_,opts__:{ColorFunction->cbarPHI}]:=Graphics[Table[{(ColorFunction/.opts)[ t/(2Pi)],wedge[r,r*dr,t,2 Pi/180,3,o]},{t,0,2 Pi,2 Pi/180}]]

Remove[wheel]

wheel[{0,0},1,0.5,ColorFunction->cbarPHI]

Module[{f},
cbarPHI[alpha_,range_:{0,1},opfun_:((1.0)&)]:=
Directive[f=Min[1,Max[0,(alpha-range[[1]])/(range[[2]]-range[[1]])]];
Blend[{Black,White,RGBColor[0,0.4,0.9],RGBColor[1.0,0.3,0]},{0.5 Cos[f*Pi+Pi/4+Pi/4]^4,1Sin[ f Pi+Pi/4+Pi/4]^8,Cos[f Pi+Pi/4]^4,Cos[f Pi+Pi/2+Pi/4]^4}],Opacity[opfun[(alpha-range[[1]])/(range[[2]]-range[[1]])]]
]
];

plot markers

plotmark[style_:linestyles[][[1]],stylein_:Directive[White], sz_:2]:=Graphics[{EdgeForm[style],FaceForm[stylein],#}]&/@{
Disk[{0,0},Offset[{sz,sz}]],
Polygon[{Offset[{0,sz}*Sqrt[Pi/2]],Offset[{sz,0}*Sqrt[Pi/2]],Offset[{0,-sz}*Sqrt[Pi/2]],Offset[{-sz,0}*Sqrt[Pi/2]]}]
}

2D symmetrized oscillator

phi[n_,x_,a_:1]:=1/Sqrt[2^n n!] (a/Pi)^(1/4) Exp[-a x^2/2] HermiteH[n,Sqrt[a] x]

phiBOLTZ[n1_,n2_,x1_,x2_]:=phi[n1,x1]phi[n2,x2]
phiBOSE[n1_,n2_,x1_,x2_]:=(phi[n1,x1]phi[n2,x2]+phi[n1,x2]phi[n2,x1])/Sqrt[2]
phiFERMI[n1_,n2_,x1_,x2_]:=(phi[n1,x1]phi[n2,x2]-phi[n1,x2]phi[n2,x1])/Sqrt[2]

f1d=FramedGraphics[Plot[{phi[1,x],phi[2,x],0},{x,-5,5},PlotStyle->{Red,Blue,Gray},PlotRange->{-0.7,0.7}],FrameLabel->{"x","ψ(x)"},framestyle,ImageSize->onecol]

Export[DROOT<>"f12-1d.pdf",f1d,"PDF"]

PWF2D[f_]:=FramedGraphics[RasterPlot[ContourPlot[f[x,y],{x,-5,5},{y,-5,5},PlotRange->{-0.5,0.5},ColorFunctionScaling->False,Contours->10,ColorFunction->(cbarBWR[#,{-1,1}]&)]],FrameLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},framestyle,ImageSize->onecol]

fboltz=PWF2D[phiBOLTZ[1,2,#1,#2]&]

Export[DROOT<>"f12-boltz.pdf",fboltz,"PDF"]

fboltz=PWF2D[phiBOLTZ[2,1,#1,#2]&]

Export[DROOT<>"f12-boltzb.pdf",fboltz,"PDF"]

fbose=PWF2D[phiBOSE[1,2,#1,#2]&]

Export[DROOT<>"f12-bose.pdf",fbose,"PDF"]

ffermi=PWF2D[phiFERMI[1,2,#1,#2]&]

Export[DROOT<>"f12-fermi.pdf",ffermi,"PDF"]

pdx=Integrate[phiBOLTZ[1,2,x,x+dx]^2,{x,-Infinity,Infinity},Assumptions->dx>0]

pdxb=Integrate[phiBOSE[1,2,x,x+dx]^2,{x,-Infinity,Infinity},Assumptions->dx>0]

pdxf=Integrate[phiFERMI[1,2,x,x+dx]^2,{x,-Infinity,Infinity},Assumptions->dx>0]

fcc=FramedGraphics[Plot[{pdx,pdxb,pdxf},{dx,-8,8},PlotStyle->{Black,Blue, Red}],FrameLabel->{"Δ","P(Δ)"},framestyle,ImageSize->onecol]

Export[DROOT<>"f12-delta.pdf",fcc,"PDF"]

Comparison of Distributions

nn=Integrate[1/(Exp[β (ϵ-μ)]-1),{ϵ,0,Infinity},Assumptions->β>0]

nn=Integrate[1/(Exp[β (ϵ-μ)]+1),{ϵ,0,Infinity},Assumptions->β>0]

nn=Integrate[Exp[-β (ϵ-μ)],{ϵ,0,Infinity},Assumptions->β>0]

Solve[nn==NN,μ]

Log[1-Exp[10]]

dnBOSE[β_,ϵ_,n_:1]:=1/(Exp[β (ϵ-μ)]-1)/.μ->Log[1-Exp[-n β]]/β
dnFERMI[β_,ϵ_,n_:1]:=1/(Exp[β (ϵ-μ)]+1)/.μ->Log[Exp[n β]-1]/β
dnBOLTZ[β_,ϵ_,n_:1]:=Exp[-β (ϵ-μ)]/.μ->Log[β n]/β

b=0.1;n=10;
Manipulate[
Plot[{dnBOSE[b,x,n],dnBOLTZ[b,x,n],dnFERMI[b,x,n]},{x,0,50},PlotRange->{{0,n*1.1+3/b},{0,n*b+1}},PlotStyle->{Blue,Black,Red},Frame->True,FrameLabel->{"ϵ","n(ϵ)"}],
{{b,2,"β"},0.02,10},
{{n,5,"N"},1,50},SaveDefinitions->True
]