#!/bin/bash

T=310 # Temperature
IT=100 # number of bootstraps

rm cvs.txt; for i in `seq $IT`; do cat ../0-colvars/cvs.txt >> cvs.txt; done

paste ../2-fe/bootstrap.txt cvs.txt | \
awk -v dx=1 -v T=$T '
function acos(o) {
    if(o==1||o==-1) {
	return 2*(1-o)*atan2(0.5,0.5)
    } else {
	return atan2(-o,sqrt(1-o*o))+2*atan2(0.5,0.5)
    }
}function round(x) {
    ival = int(x)
    if (ival == x) {
	return x
    } else {
	if (x < 0) {
	    aval = -x
	    ival = int(aval)
	    fraction = aval - ival
	    if (fraction >= 0.5)
		return int(x)-1
	    else
		return ((int(x)==-0)?0:int(x))
	} else {
	    fraction = x - ival
	    if (fraction >= .5)
		return ival + 1
	    else
		return ival
	}
    }
}BEGIN{
    x0=1e8
    x1=-1e8
}{
    t=$1; IT=t
    for (i=1;i<=12;i++) {
	x=round((360/acos(-1))*acos($(5+4*i))/dx)
	p0[t,i,x]++
	p[t,i,x]+=$5
	if(x<x0)x0=x
	if(x>x1)x1=x
    }
}END{
    KBT=0.0019872041*T
    for(x=x0;x<=x1;x++) {
	printf "%s ",x*dx
	for(i=1;i<=12;i++) {f=0; f2=0; n=0; for(t=1;t<=IT;t++) if (p[t,i,x]>0) {f+=log(p[t,i,x]); f2+=log(p[t,i,x])*log(p[t,i,x]); n++} printf "%s %s ",-KBT*f/n,KBT*sqrt(n*f2-f*f)/n}
	for(i=1;i<=12;i++) {f=0; f2=0; n=0; for(t=1;t<=IT;t++) if (p0[t,i,x]>0) {f+=log(p0[t,i,x]); f2+=log(p0[t,i,x])*log(p0[t,i,x]); n++} printf "%s %s ",-KBT*f/n,KBT*sqrt(n*f2-f*f)/n}
        printf "\n"
    }
}' > pmf.err

