#!/bin/bash

##collect the generalised forces for each value of lambda
for system in C7eq alphaL C7ax alphaR 
do

  rm -f $system/gF_all.dat

  for lambda in 0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
  do  
    f=$(ls $system/emil.$system.$lambda.log)
    wc -l $f
     
    ##get the total generalised force at the given timepoint and lambda
    grep "nstep:" $f |\
     awk '{print '$lambda',$4+$6+$8}' > $system/gF$lambda.dat

    ##take the time-average generalised force also
    awk 'NF==2{t+=$2;c++}END{print '$lambda',t/c}' \
          $system/gF$lambda.dat >> $system/gF_all.dat

  done
done

ref=C7eq
for system in alphaL C7ax alphaR 
do

rfiles=$(ls -rt $ref/gF[0-9]*.dat)
sfiles=$(ls -rt $system/gF[0-9]*.dat)
files="$rfiles $sfiles"

echo "#Forces from files:"
for f in $files
do
  echo "# "$f
done

##Do the integration timepoint by timepoint
##using just a simple Euler method for this example.
paste $files | \
 awk 'NF==36{lx=$1;ly=($20-$2);t=0.0;\
             for(i=1;i<=8;i++){\
                               x=$(1+2*i);h=x-lx;lx=x;\
                               y=$(20+2*i)-$(2+2*i);\
                               t+=h*0.5*(y+ly);ly=y;}\
             print t}' > ddA_$system.dat

done

##make a plot and find the means
gnuplot --persist <<EOF
f(x)=A
g(x)=B
h(x)=C
fit f(x) './ddA_alphaR.dat' via A
fit g(x) './ddA_C7ax.dat'   via B
fit h(x) './ddA_alphaL.dat' via C

plot './ddA_alphaR.dat'  w l,\
     './ddA_C7ax.dat'    w l,\
     './ddA_alphaL.dat'  w l, A, B, C
print A, B, C
EOF