#!/bin/bash

T=310 # Temperarure
dx=1 # bin size

paste ../2-fe/density.txt ../0-colvars/cvs.txt | \
awk -v dx=$dx -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
}
{
    Pr=$3
    for (i=1;i<=12;i++) {
	x=round((360/acos(-1))*acos($(3+4*i))/dx)
	p0[i,x]++
	p[i,x]+=Pr
	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++) printf "%s ",-KBT*log(p[i,x])
	for(i=1;i<=12;i++) printf "%s ",-KBT*log(p0[i,x])
        printf "\n"
    }
}' > pmf.txt

