added angle calc (acos problem still existing)

diff --git a/vasp_tools/angle_calc b/vasp_tools/angle_calc
new file mode 100755 (executable)
index 0000000..da01dfb
--- /dev/null
+++ b/vasp_tools/angle_calc
@@ -0,0 +1,99 @@
+#!/bin/bash
+
+file=$1
+atom=$2
+btom=$3
+ctom=$4
+
+((offset=8))
+
+scale=`sed -n 2p $file`
+X1=`sed -n 3p $file | awk '{ print $1 }'`
+X2=`sed -n 3p $file | awk '{ print $2 }'`
+X3=`sed -n 3p $file | awk '{ print $3 }'`
+Y1=`sed -n 4p $file | awk '{ print $1 }'`
+Y2=`sed -n 4p $file | awk '{ print $2 }'`
+Y3=`sed -n 4p $file | awk '{ print $3 }'`
+Z1=`sed -n 5p $file | awk '{ print $1 }'`
+Z2=`sed -n 5p $file | awk '{ print $2 }'`
+Z3=`sed -n 5p $file | awk '{ print $3 }'`
+export X1 X2 X3
+export Y1 Y2 Y3
+export Z1 Z2 Z3
+
+((line1=atom+offset))
+((line2=btom+offset))
+((line3=ctom+offset))
+
+temp="`sed -n ${line1}p $1`"
+xa=`echo $temp | awk '{ print $1 }'`
+ya=`echo $temp | awk '{ print $2 }'`
+za=`echo $temp | awk '{ print $3 }'`
+
+temp="`sed -n ${line2}p $1`"
+xb=`echo $temp | awk '{ print $1 }'`
+yb=`echo $temp | awk '{ print $2 }'`
+zb=`echo $temp | awk '{ print $3 }'`
+
+temp="`sed -n ${line3}p $1`"
+xc=`echo $temp | awk '{ print $1 }'`
+yc=`echo $temp | awk '{ print $2 }'`
+zc=`echo $temp | awk '{ print $3 }'`
+
+echo -en "angle: "
+echo "$xa $ya $za $xb $yb $zb $xc $yc $zc $scale" | \
+       awk ' \
+       BEGIN {
+               X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"]
+               Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"]
+               Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"]
+               pi=3.14159265
+       }       
+       {
+               X=sqrt(X1^2+X2^2+X3^2)
+               Y=sqrt(Y1^2+Y2^2+Y3^2)
+               Z=sqrt(Z1^2+Z2^2+Z3^2)
+               dx=$1-$4
+               dy=$2-$5
+               dz=$3-$6
+               if(dx>1/2)
+                       dx-=1
+               if(dx<-1/2)
+                       dx+=1
+               if(dy>1/2)
+                       dy-=1
+               if(dy<-1/2)
+                       dy+=1
+               if(dz>1/2)
+                       dz-=1
+               if(dz<-1/2)
+                       dz+=1
+               Dx=$1-$7
+               Dy=$2-$8
+               Dz=$3-$9
+               if(Dx>1/2)
+                       Dx-=1
+               if(Dx<-1/2)
+                       Dx+=1
+               if(Dy>1/2)
+                       Dy-=1
+               if(Dy<-1/2)
+                       Dy+=1
+               if(Dz>1/2)
+                       Dz-=1
+               if(Dz<-1/2)
+                       Dz+=1
+               dxt=dx*X1+dy*Y1+dz*Z1
+               dyt=dx*X2+dy*Y2+dz*Z2
+               dzt=dx*X3+dy*Y3+dz*Z3
+               Dxt=Dx*X1+Dy*Y1+Dz*Z1
+               Dyt=Dx*X2+Dy*Y2+Dz*Z2
+               Dzt=Dx*X3+Dy*Y3+Dz*Z3
+               sp=dxt*Dxt+dyt*Dyt+dzt*Dzt
+               d=sqrt(dxt^2+dyt^2+dzt^2)
+               D=sqrt(Dxt^2+Dyt^2+Dzt^2)
+               print 180.0/pi*acos2(sp/(d*D))
+       }'
+
+echo
+