# gnuplot-file
#
# s=1/2 Heisenberg model for the square lattice
#   finite size + periodic boundary conditions (PBC)
# see also
#   http://deposit.ddb.de/cgi-bin/dokserv?idn=967430151
#   http://diglib.uni-magdeburg.de/Dissertationen/2003/joeschulenburg.htm
#   http://www-e.uni-magdeburg.de/jschulen/diss/diss.pdf
#
#set term dumb (grep -v ps=0 for eps-files)
# E     = A0      + A3*L^-3
# GAP   = D4*L^-2 + D5*L^-3
# N*GAP = D4      + D5*L^-1
# Mag   = M0      + M1*L^-1
# gute Uebersicht in PRB45_7229_runge92
# http://www.comp-phys.org:16080/lugano04/Talks/CoupledCluster.pdf
#
ps=1
# outcommend next line for eps-output
ps=0
#
set key left
set pointsize 1.0
d4=12.1936
a4=d4/4
b4=a4+d4
e0(x)=a0 + a3*x+ a4*x**(4/3.)
e1(x)=b0 + b3*x+ b4*x**(4/3.)
g(x)=            d4          + d5*x
g2(x)=           d42         + d52*x
m(x)=            m0          + m1*x 
m2(x)=           m02         + m12*x + m22*x**2 + m23*x**3
fx(x)=x**(-3/2.)
# fit Okt2004 36er?
a0= -0.667504 # +/- 0.0007583    (0.1136%)
a3= -3.18117  # +/- 0.3775       (11.87%)
a4= 4.91518   # +/- 1.386        (28.2%)
#fit [0:fx(17)] e0(x) "square.dat" i 2 u (fx($1)):($6/$1) via a0,a3,a4
fit [0:fx(29)] e0(x) "square.dat" i 2 u (fx($1)):($6/$1) via a0,a3,a4
b0= -0.668239 # +/- 0.0006648    (0.09948%)
b3= -2.15287  # +/- 0.331        (15.37%)
b4= 9.9464    # +/- 1.215        (12.22%)
#fit [0:fx(17)] e1(x) "square.dat" i 2 u (fx($1)):($7/$1) via b0,b3,b4
fit [0:fx(29)] e1(x) "square.dat" i 2 u (fx($1)):($7/$1) via b0,b3,b4
#fit e0(x) "square.dat" i 1 u (fx($1)):($6/$1):($6/$1-e0(fx($1))) via a0,a3,a4

g31= 0.0145842 # +/- 0.0101       (69.24%)
g32= 11.9428   # +/- 0.8022       (6.717%)
g33= -12.9514  # +/- 2.674        (20.65%)
g3(x)= g31 + g32*x + g33*x**(3./2)
#fit [:1./17] g3(x) "square.dat" i 1 u (1./$1):($7-$6) via g31,g32,g33

fx(x)=x**(-1/2.)
d42= 12.9348  # +/- 0.04216      (0.3259%)
d52= -16.0839 # +/- 0.2139       (1.33%)
#fit g2(x) "square.dat" i 2 u (fx($1)):(($7-$6)*$1) via d42,d52
m02= 0.0489373 # +/- 0.006436     (13.15%)
m12= 3.22294   # +/- 0.07945      (2.465%)
m22= -7.00379  # +/- 0.3175       (4.534%)
m23= 6.77145   # +/- 0.4111       (6.071%)
#fit m2(x) "square.dat" i 3 u (fx($1)):(($8)**.5) via m02,m12,m22,m23
m0= 0.317271 # +/- 0.001053     (0.332%)
m1= 0.851629 # +/- 0.005343     (0.6274%)
#fit m(x)  "square.dat" i 2 u (fx($1)):(($8)**.5) via m0,m1
# fit taken by betts
# i0,1,2,3... = all,squared,highestE,(1,3),...
d4 =  14.3135        
d5 = -24.3338       
#fit g(x) "-" u (fx($1)):(($7-$6)*$1) via d4,d5
fit r0+r1*x "-" u ($1**(-1/2.)):($8**.5) via r0,r1
#32   2  6    0 16   -21.76397307 -21.44913413 0.21926565
#34   1 13    0 34   -23.09745020 -22.79798839 0.21467480
36   3  5    0 12   -24.42883226 -24.143307   0.21071424 # ll=36,37,45
#38   1  7    0 38   -25.7607925  -25.4880015  0.20751801
40   2  6    0 20   -27.09485025 -26.83320    0.20361937 # ??? nicht maxE
e

if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_erg.eps"
set key left bottom
set pointsize 1.5
set key 0.0285,-0.67 Right spacing 1.1 
fx(x)=x**(-3/2.)   # x-trafo
fi(x)=x**(-2/3.)   # invers trafo
yy=-0.648
set label 5 "N:" at fx(56),yy center
set label 1 "12" at fx(12),yy center
set label 2 "16" at fx(16),yy center
set label 3 "24" at fx(24),yy center
set label 4 "36" at fx(36),yy center
set xlabel "L^{-3}"
set ylabel "E/N"
# -0.66806-2.4884*x+1.754*x**(4/3.),
# "square.dat" i 5 u ($1**(-3/2.)):($6/$1) t "Eodd",\
# x==1/L^3 SWT 3th oder PRB47weihong93 (2.11)
#  pades-koeff-abschaetzung (serie?) e=-0.6693
esw(x)=-0.6693-(4.06656*.5+0.321151)*x
# lin01 ED+QMC
elin01(x)=-0.66944-2.272*x+1.64*x**(4./3)
#
# L=(1,3),(0,l) => J1-J3-chain
# L=(2,2),(0,l) => S=1-chain E=a0+a2/L^2 L=N (!) PRB55_2721qin97
#  weihong93 (2.6) S=1 chain
#  es1(L,s=1)      =-1.396391525  -2.474925987/L^2
#  es1(L,s=1/2+1/2)=-1.396391525/2-2.474925987/L^2
s=1
pi=3.141592654
# Weihong93: 1D-chain SWT E(L) es1(1/0)=-1.39639  =>  (es1/2)=-0.6982
es1(x)=(-s*s+s*(2/pi-1)-1./4*(2/pi-1)**2 - pi/3*(x**(-2.))*(2*s-2/pi+1) )/2
#  es1(fi(x))  Faktoren???
#    irgendwie scheint die Potenz der FS-Korrekturen nie zu stimmen ???
# "< echo 0 -0.700742019485" t "Haldane E/2N" w p 1 12,
#
# plot only N>8, its better looking
plot [0:fx(9.8)] [-0.732:-0.644]  \
e0(x) t "" w l 1 1,e1(x) t "" w l 1 1,\
 "square.dat" i 0 u (fx($1)):($6/$1) t "E_0/N" w p 1 1,\
 "square.dat" i 1 u (fx($1)):($6/$1) t "" w p 1 6,\
 "square.dat" i 0 u (fx($1)):($7/$1) t "E_1/N" w p 1 2,\
 "square.dat" i 1 u (fx($1)):($7/$1) t "" w p 1 6,\
 -0.6907 -3.6 *fi(x)**(-2.) t "(1,3),(0,l)" w l 3 3,\
 -0.6919 +1.45*fi(x)**(-3/2.) t "" w l 3 3,\
 -0.700742-30 *fi(x)**(-3.) t "(2,2),(0,l)" w l 4 4,\
 -0.697  +1.65*fi(x)**(-3/2.) t "" w l 4 4,\
elin01(x) t "Lin01" w l 5 5,\
esw(x) t "Weihong93" w l 2 2,\
-0.66960*((x>0.001)?1/0:1) t "CCM extrapolated" w l 7 7,\
0 t ""
# (2,2) exact S=1-chain! scaling? E0/N==E1/N !!!???
#"square.dat" i 5 u ($1**(-3/2.)):($6/$1) t "E0/N" w p 8 8,
if( ps==0 ) pause -1

if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_erg2.eps"
yy=-0.666
set label 5 "N:" at fx(64),yy center
set label 1 "40" at fx(40),yy center
set label 2 "36" at fx(36),yy center
set label 3 "32" at fx(32),yy center
set label 4 "28" at fx(28),yy center
set key left bottom Right spacing 1.1 
plot [0:fx(23.8)] [-0.685:-0.665]  \
e0(x) t "" w l 1 1,e1(x) t "" w l 1 1,\
 "square.dat" i 0 u (fx($1)):($6/$1) t "E_0/N" w p 1 1,\
 "square.dat" i 1 u (fx($1)):($6/$1) t "" w p 1 6,\
 "square.dat" i 0 u (fx($1)):($7/$1) t "E_1/N" w p 1 2,\
 "square.dat" i 1 u (fx($1)):($7/$1) t "" w p 1 6,\
 -0.6907 -3.6 *fi(x)**(-2.) t "(1,3),(0,l)" w l 3 3,\
 -0.6919 +1.45*fi(x)**(-3/2.) t "" w l 3 3,\
 -0.700742-30 *fi(x)**(-3.) t "(2,2),(0,l)" w l 4 4,\
 -0.697  +1.65*fi(x)**(-3/2.) t "" w l 4 4,\
elin01(x) t "Lin01" w l 5 5,\
esw(x) t "Weihong93" w l 2 2,\
-0.66960*((x>0.001)?1/0:1) t "CCM extrapolated" w l 7 7,\
0 t ""
if( ps==0 ) pause -1
if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_gap2.eps"
set key bottom right
set pointsize 2
yy=0.16
fx(x)=x**(-1.)
set label 5 "N:" at fx(56),yy center
set label 1 "40" at fx(40),yy center
set label 2 "20" at fx(20),yy center
set label 3 "26" at fx(26),yy center
set label 4 "36" at fx(36),yy center
set xlabel "N^{-1}"
set ylabel "{/Symbol D}"
set ylabel "E-E_0"
set xtics ("0.1" 0.1, "0.05" 0.05, "0.025" 0.025, "0.0125" 0.0125)
set ytics ("4" 4, "2" 2, "1" 1, "0.5" 0.5, "0.25" 0.25, "0.125" 0.125)
set logscale xy
### fit r0*x**r1 "square.dat" i 1 u (1./$1):($10-$6) via r0,r1
# Singlett-Tripplet-GAP \sim N**-1    5.64(10)*x**0.830(7)   q=(pi,pi)-
# lowest-Magnon-Gap     \sim N**-1/2  7.26(68)*x**0.400(30)  q=(x,x)-
# Singlett-Singlett-GAP \sim N**-1/4  5.16(6)*x**0.198(4)    q=(0,pi)+ N%4
plot [0.11*.15:.37**2] [0.11:4*1.1]\
 5.64*x**0.830 t "", 7.26*x**0.400 t "", 5.16*x**0.198 t "", "square.dat"\
   i 1 u (1./$1):($7-$6)  t "E_0(S=1)-E_0" w p 1 6,\
"" i 1 u (1./$1):($9-$6)  t "E_1(S=1)-E_0" w p 2 7,\
"" i 1 u (1./$1):($10-$6) t "E_1(S=0)-E_0" w p 2 8
if ( ps==0 ) pause -1
set nologscale xy
set xtics auto
set ytics auto

if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_gap1.eps"
set key top left
set pointsize 2
yy=0.05
fx(x)=1./x
set label 5 "N:" at fx(40),yy center
set label 1 "10" at fx(10),yy center
set label 2 "16" at fx(16),yy center
set label 3 "24" at fx(24),yy center
set label 4 "36" at fx(36),yy center
set xlabel "L^{-2}"
set ylabel "{/Symbol D}"
plot [0*.15:fx(7.8)] [0:1.1] g2(x**(.5))*x t "", g3(x) t "",\
 "square.dat" i 5 u (fx($1)):(2./3*($7-$6)) t "2/3 (E(S=3/2)-E(S=1/2))",\
 "square.dat" i 3 u (fx($1)):($7-$6) smooth csplines t "(1,3),(0,l)" w l 3 3,\
 "square.dat" i 4 u (fx($1)):($7-$6) smooth csplines t "(2,2),(0,l)" w l 4 4,\
 "square.dat" i 0 u (fx($1)):($7-$6) t "{/Symbol D}" w p 1 1,\
 "square.dat" i 1 u (fx($1)):($7-$6) t "" w p 1 6
if ( ps==0 ) pause -1


if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_gap.eps"
set key right
set pointsize 2
yy=7.8
fx(x)=x**(-1/2.)
set label 5 "N:" at fx(64),yy center
set label 1 "10" at fx(10),yy center
set label 2 "16" at fx(16),yy center
set label 3 "24" at fx(24),yy center
set label 4 "36" at fx(36),yy center
# D=j(j+1)/(\chi L^2) mit j=1
set arrow from 0.03,14.404 to 0,14.404
set label "Weihong93" at 0.031,14.2
set xlabel "L^{-1}"
set ylabel "N{/Symbol D}"
#
plot [0*.15:fx(7.8)] [7.5:14.5] g2(x) t "fit 18...36",g(x) t "fit 38...40",\
 "square.dat" i 3 u (fx($1)):(($7-$6)*$1) smooth csplines t "(1,3),(0,l)" w l 3 3,\
 "square.dat" i 4 u (fx($1)):(($7-$6)*$1) smooth csplines t "(2,2),(0,l)" w l 4 4,\
 "square.dat" i 0 u (fx($1)):(($7-$6)*$1) t "" w p 1 1,\
 "square.dat" i 1 u (fx($1)):(($7-$6)*$1) t "" w p 1 6
if ( ps==0 ) pause -1

if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_gap3.eps"
set key right
set pointsize 2
yy=7.8
fx(x)=x**(-1/2.)
set label 5 "N:" at fx(64),yy center
set label 1 " 9" at fx( 9),yy center
set label 2 "15" at fx(15),yy center
set label 3 "21" at fx(21),yy center
set label 4 "27" at fx(27),yy center
set xlabel "L^{-1}"
set ylabel "N{/Symbol D}"
plot [0*.15:fx(7.8)] [7.1:14.5] \
14.4-22.2*x t "",\
"square.dat" i 5 u (fx($1)):(2./3*($7-$6)*$1) t "2/3 N {/Symbol D}",\
"<awk '{if($1%2==1 && $2==1 && $3==3)print}' square.dat" u (fx($1)):(2./3*($7-$6)*$1) smooth csplines t "(1,3),(0,l)" w l,\
"<awk '{if($1%2==1 && $2==1 && $3==4)print}' square.dat" u (fx($1)):(2./3*($7-$6)*$1) smooth csplines t "(1,4),(0,l)" w l,\
"<awk '{if($1%2==1 && $2==1 && $3==5)print}' square.dat" u (fx($1)):(2./3*($7-$6)*$1) smooth csplines t "(1,5),(0,l)" w l,\
"<awk '{if($1%2==1 && $2==1 && $3==7)print}' square.dat" u (fx($1)):(2./3*($7-$6)*$1) smooth csplines t "(1,7),(0,l)" w l
# "square.dat" i 4 u (fx($1)):(($7-$6)*$1) smooth csplines t "(2,2),(0,l)" w l 4 4,\
#
if ( ps==0 ) pause -1
set nolabel
set noarrow

if( ps==1 ) set term postscript eps enhanced mono dashed "Helvetica" 20; set out "tmp/sq_mag.eps"
set key left
yy=0.32
fx(x)=x**(-1/2.)
set label 5 "N:" at fx(64),yy center
set label 1 "10" at fx(10),yy center
set label 2 "16" at fx(16),yy center
set label 3 "24" at fx(24),yy center
set label 4 "36" at fx(36),yy center
set xlabel "L^{-1}"
set ylabel "m^+"
# weihong93 x=1/L (2.12)
msw(x)=0.3037+0.8778680*x
# lin01 x=1/L (8) mz^2   m=(3mz)**.5=0.30530313
mLin01z2(x)=0.03107+0.2009*x+0.196*x**2+0.007*x**3
plot [0*.14:.37] [0.99*0.3:0.65] m(x) t "", m2(x) t "(1,3),(0,l)" w l 3 3,\
 "square.dat" i 4 u (fx($1)):($8**.5) t "(2,2),(0,l)" w l 4 4,\
 "square.dat" i 0 u (fx($1)):($8**.5) t "m^+" w p 1 1,\
 "square.dat" i 1 u (fx($1)):($8**.5) t "" w p 1 6,\
 msw(x) t "Weihong93" w l 2 2, (3*mLin01z2(x))**.5 t "Lin01" w l 5 5,\
 0.610/2*((x>0.01)?1/0:1) t "CCM extrapolated" w l 7 7
if ( ps==0 ) pause -1