#! /bin/sh
hstep=$1
facc=$2
nconnect=$3
nburn=$4
wait=$5
plotno=$6
date=`date '+%D %H:%M'`
sig1=2
#sig2=0.5
sig2=0.1
alpha=30

sm -s << FIN
if ($wait) {xterm_ss} else {laser_sp mc2d.$plotno.ps\n lweight 2}
define alpha ($alpha*pi/180)
define cosa (cos(\$alpha))
define sina (sin(\$alpha))
define sigx2 ($sig1**2*\$cosa*\$cosa+$sig2**2*\$sina*\$sina)
define sigy2 ($sig1**2*\$sina*\$sina+$sig2**2*\$cosa*\$cosa)
expand 1.1
limits -8 8 -8 8
box

da mc2d.out
read {x 1 y 2}
set i=0,dimen(x)-1
ctype 3
if ($nconnect>0) {
    connect x y if (i<$nconnect)
}
set x=x if (i>=$nburn)
set y=y if (i>=$nburn)
ptype 1 1
points x y
set xh=-6,6,0.2
set yh=histogram(x:xh)
vecminmax yh min xmax
set yh = yh/\$xmax
limits -8 8 0 4
histogram xh yh
set yh=histogram(y:xh)
vecminmax yh min ymax
set yh = yh/\$ymax
limits 0 4 -8 8
histogram angle yh xh

# plot the original 2-d Gaussian
ctype 0
limits -8 8 -8 8
da gauss2d.out2
if ($sig2==0.1)  {
    da gauss2d.out2b
}
read {x 1 y 2}
ptype 1 1
points x y
set xh=-6,6,0.2
set yh=histogram(x:xh)
vecminmax yh min xmax
set yh = yh/\$xmax
limits -8 8 0 4
histogram xh yh
set yh=histogram(y:xh)
vecminmax yh min ymax
set yh = yh/\$ymax
limits 0 4 -8 8
histogram angle yh xh
ctype 0

# plot the starting point and first Nconnect points last so not overwritten
limits -8 8 -8 8
da mc2d.out
read {x 1 y 2}
set i=0,dimen(x)-1
ctype 3
if ($nconnect>0) {
    connect x y if (i<$nconnect)
}
expand 1.5
ptype 20 3
rel (x[0]) (y[0])
dot
define xi (x[0])
define yi (y[0])

expand 1.3
ctype 0
limits 0 10 0 10
relocate 5 -1.1
putlabel 5 x
relocate -1.1 5
putlabel 5 y
relocate 0.5 9.5
putlabel 6 Black: bivariate Gaussian
relocate 0.5 9
putlabel 6 Red: MCMC emulation
relocate 0.5 8.5
putlabel 6 x_i=\$xi, y_i=\$yi, h = $hstep
relocate 0.5 8.0
putlabel 6 N_{burn}=$nburn, f_{accept}=$facc

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