#!/bin/ksh

cat << EOF > script.gnuplot
set terminal postscript
set size square
set output "test.ps"
EOF

echo
for npts in 1000
do
  for rad in 20 40 60 80
  do
    suff=${npts}_${rad}
    cat hardspheres.f90 | 	sed "s/n     = 10000/n     = $npts/g" | \
				sed "s/rad   = 0.25d0\/sqrt(dble(n))/rad   = 0.${rad}d0\/sqrt(dble(n))/g" | \
				sed "s/nstop = 100000/nstop = 1000000/g" \
				> dummycodexflr6.f90
    gfortran dummycodexflr6.f90
    echo $npts points, rad=0.${rad}d0\/sqrt\(n\)
    ./a.out > test${suff}.dat
    rm dummycodexflr6.f90 a.out
ndone=$(tail -1 test${suff}.dat | awk '{print $4;}')
cat <<EOF >> script.gnuplot
set title "$ndone pts (asked ${npts}...), dmin 2*$rad (arb.)  \n rejects=f(n)"
set xrange [0:*]; set yrange [0:*]
plot 'test${suff}.dat' u 0:5 w impulses notitle

set title "$ndone pts (asked ${npts}...), dmin = 2*rad of circles ($rad) \n [0:1]x[0:1] coverage"
set xrange [0:1]; set yrange [0:1]
plot 'test${suff}.dat' u 1:2:3 w circles notitle

set title "$ndone pts (asked ${npts}...), dmin = 2*rad of circles ($rad) \n [0:1]x[0:1] distribution"
set xrange [0:1]; set yrange [0:1]
plot 'test${suff}.dat' u 1:2 w p pt 7 ps 1 notitle
EOF
  done
done

gnuplot script.gnuplot
ps2pdf test.ps

#clean
for fxflr6 in test*.dat test.ps script.gnuplot; do rm $fxflr6; done

evince test.pdf