#!/bin/bash

T=310 # Temperature
IT=100 # number of bootstraps
l=`wc ../0-colvars/cvs.txt | awk '{print $1}'` # number of data points
dx=1 # bin size

cat ../0-colvars/cvs.txt ../2-fe/bootstrap.txt | \
awk -v dx=$dx -v T=$T -v l=$l '
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
}{
    if (NR<=l){
	t=$1
	for (i=1;i<=12;i++) {
	    x=round((360/acos(-1))*acos($(4*i))/dx)
	    X[i,t]=x
	    if(x<x0)x0=x
	    if(x>x1)x1=x
	}
    } else {
	b=$1
	IT=b
	t=$4
	Pr=$5
	for (i=1;i<=12;i++) {
	    x=X[i,t]
	    p0[b,i,x]++
	    p[b,i,x]+=Pr
	}
    }
}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(b=1;b<=IT;b++) if (p[b,i,x]>0) {f+=log(p[b,i,x]); f2+=log(p[b,i,x])*log(p[b,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(b=1;b<=IT;b++) if (p0[b,i,x]>0) {f+=log(p0[b,i,x]); f2+=log(p0[b,i,x])*log(p0[b,i,x]); n++} printf "%s %s ",-KBT*f/n,KBT*sqrt(n*f2-f*f)/n}
        printf "\n"
    }
}' > pmf.err

