#! /bin/sh
nconnect=0
nburn=0
wait=$1
plotno=$2
date=`date '+%D %H:%M'`

sig1=	# up to you to compute these values and insert them
sig2=

xmin=-1
xmax=5
ymin=1.6
ymax=2.4

sm -s << FIN
if ($wait) {xterm_p} else {laser_pp linemc.$plotno.ps\n lweight 2}

# plot the MCMC points for the first random number seed
sixwin 1
expand 0.9
limits $xmin $xmax $ymin $ymax
box
da line.n20.s12.mcmc
read {x 1 y 2 p 3}
set i=0,dimen(x)-1
ctype 3
if ($nconnect>0) {
    connect x y if (i<$nconnect)
}
ctype 0
set x=x if (i>=$nburn)
set y=y if (i>=$nburn)
ptype 1 1
points x y

# plot marginal histograms along the axes
# and compare to Gaussians with predicted width
set xh=-2,6,0.2
set yh=histogram(x:xh)
vecminmax yh min xmax
set yh = yh/\$xmax
limits $xmin $xmax 0 4
histogram xh yh
ctype 4
connect xh (exp(-(xh-2.)**2/(2*($sig1)**2)))
ctype 0
set xh=1.7,2.3,0.01
set yh=histogram(y:xh)
vecminmax yh min ymax
set yh = yh/\$ymax
limits 0 4 $ymin $ymax 
histogram angle yh xh
ctype 4
connect (exp(-(xh-2.)**2/(2*($sig2)**2))) xh
ctype 0

limits 0 10 0 10
relocate 5 -1.8
putlabel 5 \theta_1
relocate -1.8 5
putlabel 5 \theta_2

# plot the data
sixwin 2
da line.n20.s12.dat
limits 0 22 0 44
box
read {xx 1 yy 2 err 3}
ptype 4 0
points xx yy
errorbar xx yy err 2
errorbar xx yy err 4

# plot five random lines from the MCMC chain
set rr=random(5)
set x0=0,21,2
lweight 1
ltype 1
do i=1,5 {
    define k (int(dimen(x)*rr[\$i-1]))
    define cc (\$i+2)
    ctype \$cc
    set y0 = x[\$k]+y[\$k]*x0
    connect x0 y0
}
ltype 0
lweight 2

# find the best-fit parameters and plot this line
ctype 7
sort {p x y}
set y0 = x[dimen(x)-1]+y[dimen(x)-1]*x0
connect x0 y0
ctype 0
if ($wait) {lweight 2}
limits 0 10 0 10
relocate 5 -1.8
putlabel 5 x
relocate -2.0 5
putlabel 5 y

# repeat the above for seed 17
sixwin 3
expand 0.9
limits $xmin $xmax $ymin $ymax
box
da line.n20.s17.mcmc
read {x 1 y 2 p 3}
set i=0,dimen(x)-1
ctype 3
if ($nconnect>0) {
    connect x y if (i<$nconnect)
}
ctype 0
set x=x if (i>=$nburn)
set y=y if (i>=$nburn)
ptype 1 1
points x y
set xh=-2,6,0.2
set yh=histogram(x:xh)
vecminmax yh min xmax
set yh = yh/\$xmax
limits $xmin $xmax 0 4
histogram xh yh
ctype 4
connect xh (exp(-(xh-2.)**2/(2*($sig1)**2)))
ctype 0
set xh=1.7,2.3,0.01
set yh=histogram(y:xh)
vecminmax yh min ymax
set yh = yh/\$ymax
limits 0 4 $ymin $ymax 
histogram angle yh xh
ctype 4
connect (exp(-(xh-2.)**2/(2*($sig2)**2))) xh
ctype 0

limits 0 10 0 10
relocate 5 -1.8
putlabel 5 \theta_1
relocate -1.8 5
putlabel 5 \theta_2

# plot the data
sixwin 4
da line.n20.s17.dat
limits 0 22 0 44
box
read {xx 1 yy 2 err 3}
ptype 4 0
points xx yy
errorbar xx yy err 2
errorbar xx yy err 4

# plot five random lines from the MCMC chain
set rr=random(5)
set x0=0,21,2
lweight 1
ltype 1
do i=1,5 {
    define k (int(dimen(x)*rr[\$i-1]))
    define cc (\$i+2)
    ctype \$cc
    set y0 = x[\$k]+y[\$k]*x0
    connect x0 y0
}
ltype 0
lweight 2

# find the best-fit parameters and plot this line
ctype 7
sort {p x y}
set y0 = x[dimen(x)-1]+y[dimen(x)-1]*x0
connect x0 y0
ctype 0
if ($wait) {lweight 2}
limits 0 10 0 10
relocate 5 -1.8
putlabel 5 x
relocate -2.0 5
putlabel 5 y

# repeat the above for "seed 0" (points on line)
sixwin 5
expand 0.9
limits $xmin $xmax $ymin $ymax
box
da line.n20.s0.mcmc
read {x 1 y 2 p 3}
set i=0,dimen(x)-1
ctype 3
if ($nconnect>0) {
    connect x y if (i<$nconnect)
}
ctype 0
set x=x if (i>=$nburn)
set y=y if (i>=$nburn)
ptype 1 1
points x y
set xh=-2,6,0.2
set yh=histogram(x:xh)
vecminmax yh min xmax
set yh = yh/\$xmax
limits $xmin $xmax 0 4
histogram xh yh
ctype 4
connect xh (exp(-(xh-2.)**2/(2*($sig1)**2)))
ctype 0
set xh=1.7,2.3,0.01
set yh=histogram(y:xh)
vecminmax yh min ymax
set yh = yh/\$ymax
limits 0 4 $ymin $ymax 
histogram angle yh xh
ctype 4
connect (exp(-(xh-2.)**2/(2*($sig2)**2))) xh
ctype 0

limits 0 10 0 10
relocate 5 -1.8
putlabel 5 \theta_1
relocate -1.8 5
putlabel 5 \theta_2

# plot the data
sixwin 6
da line.n20.s0.dat
limits 0 22 0 44
box
read {xx 1 yy 2 err 3}
ptype 4 0
points xx yy
errorbar xx yy err 2
errorbar xx yy err 4

# plot five random lines from the MCMC chain
set rr=random(5)
set x0=0,21,2
lweight 1
ltype 1
do i=1,5 {
    define k (int(dimen(x)*rr[\$i-1]))
    define cc (\$i+2)
    ctype \$cc
    set y0 = x[\$k]+y[\$k]*x0
    connect x0 y0
}
ltype 0
lweight 2

# find the best-fit parameters and plot this line
ctype 7
sort {p x y}
set y0 = x[dimen(x)-1]+y[dimen(x)-1]*x0
connect x0 y0
ctype 0
if ($wait) {lweight 2}
limits 0 10 0 10
relocate 5 -1.8
putlabel 5 x
relocate -2.0 5
putlabel 5 y

#putdate "$date"
if ($wait) {!sleep $wait} else {hardcopy}
FIN