From: hackbard Date: Wed, 27 Jan 2010 17:08:30 +0000 (+0100) Subject: added angle calc (acos problem still existing) X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=4b51fdc9537cc3bf389d7b596bc8dbd973ad157e;p=physik%2Fposic.git added angle calc (acos problem still existing) --- diff --git a/vasp_tools/angle_calc b/vasp_tools/angle_calc new file mode 100755 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 +