Points:
A -2.08576 1.76533 -0.46417
B -0.95929 0.87554 0.03365
C 0.28069 1.66193 0.42640
D 0.62407 2.22927 -0.44649
So far, I have done:
#!/bin/bash
awk 'NR==1' $FILE > LINEA
awk 'NR==1' $FILE > LINEB
awk 'NR==1' $FILE > LINEC
awk 'NR==1' $FILE > LINED
x1=`awk '{print $2}' LINEA` # x1
y1=`awk '{print $3}' LINEA` # y1
z1=`awk '{print $4}' LINEA` # z1
x2=`awk '{print $2}' LINEB` # x2
y2=`awk '{print $3}' LINEB` # y2
z2=`awk '{print $4}' LINEB` # z2
x3=`awk '{print $2}' LINEC` # x3
y3=`awk '{print $3}' LINEC` # y3
z3=`awk '{print $4}' LINEC` # z3
x4=`awk '{print $2}' LINED` # x4
y4=`awk '{print $3}' LINED` # y4
z4=`awk '{print $4}' LINED` # z4
v1x=`calc "($x1)-($x2)" | sed 's/^\t//g'`
v1y=`calc "($y1)-($y2)" | sed 's/^\t//g'`
v1z=`calc "($z1)-($z2)" | sed 's/^\t//g'`
v2x=`calc "($x4)-($x3)" | sed 's/^\t//g'`
v2y=`calc "($y4)-($y3)" | sed 's/^\t//g'`
v2z=`calc "($z4)-($z3)" | sed 's/^\t//g'`
v1mag=`calc "sqrt(($v1x)**2+($v1y)**2+($v1z)**2)" | sed 's/^\t//g'`
v2mag=`calc "sqrt(($v2x)**2+($v2y)**2+($v2z)**2)" | sed 's/^\t//g'`
calc "acos((($v1x)/($v1mag))*(($v2x)/($v2mag))+(($v1y)/($v1mag))*(($v2y)/($v2mag))+(($v1z)/($v1mag))*(($v2z)/($v2mag)))*180/3.141592653589793" | sed 's/^\t//g' | sed 's/^~//g'
calc "acos((($x1)*($x4)+($y1)*($y4)+($z1)*($z4))/(sqrt(($x1)**2+($y1)**2+($z1)**2)*sqrt(($x4)**2+($y4)**2+($z4)**2)))*180/3.141592653589793" | sed 's/^\t//g' | sed 's/^~//g'
I have found these related links 1, 2 and 3.
The referenced value is 58.7 $^{o}$
The value that I get is: 70.62525933704842342761 $^{o}$ and 64.23010091217222985704 $^{o}$
Someone knows what would the best algorithm be to obtain it properly?
Based on your refined shell code elsewhere in this thread, I've transcribed that into an awk solution as well. As people seem to have found the _docd version of use, I will include that at the end. I'm also including a debug version (in the middle of the reply).
cat torsion2.awk
-
#!/bin/awk -f
BEGIN {
# dbg=0 # turns off dbg output
# see below for debug version of this script
}
function acos(x) { return atan2((1.-x^2)^0.5,x) }
NR==1 {x1=$2; y1=$3; z1=$4}
NR==2 {x2=$2; y2=$3; z2=$4}
NR==3 {x3=$2; y3=$3; z3=$4}
NR==4 {
x4=$2; y4=$3; z4=$4
# all of this code below is only executed when you read in the 4th line
# because then you have all the data
#
v1x=x2-x1 ; v1y=y2-y1 ; v1z=z2-z1 #plane1
v2x=x3-x2 ; v2y=y3-y2 ; v2z=z3-z2 #plane1
v3x=x2-x3 ; v3y=y2-y3 ; v3z=z2-z3 #plane2
v4x=x3-x4 ; v4y=y3-y4 ; v4z=z3-z4 #plane2
plane1_x=(v1y*v2z)-(v1z*v2y) # normal vector 1
plane1_y=(v2x*v1z)-(v2z*v1x) # normal vector 1
plane1_z=(v1x*v2y)-(v1y*v2x) # normal vector 1
plane2_x=(v3y*v4z)-(v3z*v4y) # normal vector 2
plane2_y=(v4x*v3z)-(v4z*v3x) # normal vector 2
plane2_z=(v3x*v4y)-(v3y*v4x) # normal vector 2
v1mag=sqrt(((plane1_x)**2)+((plane1_y)**2)+((plane1_z)**2)) # magnitude normal vector 1
v2mag=sqrt(((plane2_x)**2)+((plane2_y)**2)+((plane2_z)**2)) # magnitude normal vector 2
vn1x=(plane1_x)/(v1mag) ; vn1y=(plane1_y)/(v1mag) ; vn1z=(plane1_z)/(v1mag) # normalization normal vector 1
vn2x=(plane2_x)/(v2mag) ; vn2y=(plane2_y)/(v2mag) ; vn2z=(plane2_z)/(v2mag) # normalization normal vector 2
print acos((vn1x*vn2x)+(vn1y*vn2y)+(vn1z*vn2z))*180/3.141592653589793
}
Once the file is saved, you must mark the script as executable:
chmod +x ./torsion2.awk
Then you can run it with the sample data supplied:
./torsion2.awk data.txt
The output is
58.6892
Here is the full debug version. I needed it because I had editing errors like changing y2=$3 to just y=$3! (These things happen ;-/ )
cat torsion2_debug.awk
#!/bin/awk -f
BEGIN {
dbg=1 # turns on dbg output
# dbg=0 # turns off dbg output
}
function acos(x) { return atan2((1.-x^2)^0.5,x) }
NR==1 {x1=$2; y1=$3; z1=$4}
NR==2 {x2=$2; y2=$3; z2=$4}
NR==3 {x3=$2; y3=$3; z3=$4}
NR==4 {
x4=$2; y4=$3; z4=$4
if (dbg) {
print "x1="x1 "\ty1="y1 "\tz1=" z1
print "x2="x2 "\ty2="y2 "\tz2=" z2
print "x3="x3 "\ty3="y3 "\tz3=" z3
print "x4="x4 "\ty4="y4 "\tz4=" z4
}
# all of this code below is only executed when you read in the 4th line
# because then you have all the data
#
v1x=x2-x1 ; v1y=y2-y1 ; v1z=z2-z1 #plane1
v2x=x3-x2 ; v2y=y3-y2 ; v2z=z3-z2 #plane1
v3x=x2-x3 ; v3y=y2-y3 ; v3z=z2-z3 #plane2
v4x=x3-x4 ; v4y=y3-y4 ; v4z=z3-z4 #plane2
if (dbg) {
print "#dbg: v1x="v1x "\tv1y=" v1y "\tv1z="v1z
print "#dbg: v2x="v2x "\tv2y=" v2y "\tv2z="v2z
print "#dbg: v3x="v3x "\tv3y=" v3y "\tv3z="v3z
print "#dbg: v4x="v4x "\tv4y=" v4y "\tv4z="v4z
}
plane1_x=(v1y*v2z)-(v1z*v2y) # normal vector 1
plane1_y=(v2x*v1z)-(v2z*v1x) # normal vector 1
plane1_z=(v1x*v2y)-(v1y*v2x) # normal vector 1
plane2_x=(v3y*v4z)-(v3z*v4y) # normal vector 2
plane2_y=(v4x*v3z)-(v4z*v3x) # normal vector 2
plane2_z=(v3x*v4y)-(v3y*v4x) # normal vector 2
if (dbg) {
print "#dbg: plane1_x=" plane1_x "\tplane1_y=" plane1_y "\tplane1_z=" plane1_z
print "#dbg: plane2_x=" plane2_x "\tplane2_y=" plane2_y "\tplane2_z=" plane2_z
}
v1mag=sqrt(((plane1_x)**2)+((plane1_y)**2)+((plane1_z)**2)) # magnitude normal vector 1
v2mag=sqrt(((plane2_x)**2)+((plane2_y)**2)+((plane2_z)**2)) # magnitude normal vector 2
if (dbg) {
print "#dbg: v1mag=" v1mag "\tv2mag="v2mag
}
vn1x=(plane1_x)/(v1mag) ; vn1y=(plane1_y)/(v1mag) ; vn1z=(plane1_z)/(v1mag) # normalization normal vector 1
vn2x=(plane2_x)/(v2mag) ; vn2y=(plane2_y)/(v2mag) ; vn2z=(plane2_z)/(v2mag) # normalization normal vector 2
if (dbg) {
print "#dbg: " (vn1x*vn2x) " " (vn1y*vn2y) " " ((vn1z*vn2z)*180/3.141592653589793)
}
print acos((vn1x*vn2x)+(vn1y*vn2y)+(vn1z*vn2z))*180/3.141592653589793
}
And here is the transcribed shell to awk version
I highly recommend the Grymoire's Awk Tutorial to help you understand the awk programming paradigm and its built in variables like NR (Number (of) Record).
cat torsion2_docd.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
# x1=`awk '{print $2}' LINEA` # x1
# y1=`awk '{print $3}' LINEA` # y1
# z1=`awk '{print $4}' LINEA` # z1
# x2=`awk '{print $2}' LINEB` # x2
# y2=`awk '{print $3}' LINEB` # y2
# z2=`awk '{print $4}' LINEB` # z2
# x3=`awk '{print $2}' LINEC` # x3
# y3=`awk '{print $3}' LINEC` # y3
# z3=`awk '{print $4}' LINEC` # z3
# x4=`awk '{print $2}' LINED` # x4
# y4=`awk '{print $3}' LINED` # y4
# z4=`awk '{print $4}' LINED` # z4
NR==1 {x1=$2; y1=$3; z1=$4}
NR==2 {x2=$2; y2=$3; z2=$4}
NR==3 {x3=$2; y3=$3; z3=$4}
NR==4 {
x4=$2; y=$3; z4=$4
# all of this code below is only executed when you read in the 4th line
# because then you have all the data
#
# v1x=`calc "($x2)-($x1)" | sed 's/^\t//g'` #plane1
# v1y=`calc "($y2)-($y1)" | sed 's/^\t//g'` #plane1
# v1z=`calc "($z2)-($z1)" | sed 's/^\t//g'` #plane1
# v2x=`calc "($x3)-($x2)" | sed 's/^\t//g'` #plane1
# v2y=`calc "($y3)-($y2)" | sed 's/^\t//g'` #plane1
# v2z=`calc "($z3)-($z2)" | sed 's/^\t//g'` #plane1
# v3x=`calc "($x2)-($x3)" | sed 's/^\t//g'` #plane2
# v3y=`calc "($y2)-($y3)" | sed 's/^\t//g'` #plane2
# v3z=`calc "($z2)-($z3)" | sed 's/^\t//g'` #plane2
# v4x=`calc "($x3)-($x4)" | sed 's/^\t//g'` #plane2
# v4y=`calc "($y3)-($y4)" | sed 's/^\t//g'` #plane2
# v4z=`calc "($z3)-($z4)" | sed 's/^\t//g'` #plane2
v1x=x2-x1 ; v1y=y2-y1 ; v1z=z2-z1 #plane1
v2x=x3-x2 ; v2y=y3-y2 ; v2z=z3-z2 #plane1
v3x=x2-x3 ; v3y=y2-y3 ; v3z=z2-z3 #plane2
v1x=x2-x1 ; v1y=y2-y1 ; v1z=z2-z1 #plane1
v2x=x3-x2 ; v2y=y3-y2 ; v2z=z3-z2 #plane1
v3x=x2-x3 ; v3y=y2-y3 ; v3z=z2-z3 #plane2
v4x=x3-x4 ; v4y=y3-y4 ; v4z=z3-z4 #plane2
# plane1_x=`calc "($v1y)*($v2z)-($v1z)*($v2y)" | sed 's/^\t//g'` # normal vector 1
# plane1_y=`calc "($v2x)*($v1z)-($v2z)*($v1x)" | sed 's/^\t//g'` # normal vector 1
# plane1_z=`calc "($v1x)*($v2y)-($v1y)*($v2x)" | sed 's/^\t//g'` # normal vector 1
# plane2_x=`calc "($v3y)*($v4z)-($v3z)*($v4y)" | sed 's/^\t//g'` # normal vector 2
# plane2_y=`calc "($v4x)*($v3z)-($v4z)*($v3x)" | sed 's/^\t//g'` # normal vector 2
# plane2_z=`calc "($v3x)*($v4y)-($v3y)*($v4x)" | sed 's/^\t//g'` # normal vector 2
plane1_x=(v1y*v2z)-(v1z*v2y) # normal vector 1
plane1_y=(v2x*v1z)-(v2z*v1x) # normal vector 1
plane1_z=(v1x*v2y)-(v1y*v2x) # normal vector 1
plane2_x=(v3y*v4z)-(v3z*v4y) # normal vector 2
plane2_y=(v4x*v3z)-(v4z*v3x) # normal vector 2
plane2_z=(v3x*v4y)-(v3y*v4x) # normal vector 2
# v1mag=`calc "sqrt(($plane1_x)**2+($plane1_y)**2+($plane1_z)**2)" | sed 's/^\t//g'` # magnitude normal vector 1
# v2mag=`calc "sqrt(($plane2_x)**2+($plane2_y)**2+($plane2_z)**2)" | sed 's/^\t//g'` # magnitude normal vector 2
v1mag=sqrt((plane1_x)**2+(plane1_y)**2+(plane1_z)**2) # magnitude normal vector 1
v2mag=sqrt((plane2_x)**2+(plane2_y)**2+(plane2_z)**2) # magnitude normal vector 2
# vn1x=`calc "($plane1_x)/($v1mag)" | sed 's/^\t//g'` # normalization normal vector 1
# vn1y=`calc "($plane1_y)/($v1mag)" | sed 's/^\t//g'` # normalization normal vector 1
# vn1z=`calc "($plane1_z)/($v1mag)" | sed 's/^\t//g'` # normalization normal vector 1
# vn2x=`calc "($plane2_x)/($v2mag)" | sed 's/^\t//g'` # normalization normal vector 2
# vn2y=`calc "($plane2_y)/($v2mag)" | sed 's/^\t//g'` # normalization normal vector 2
# vn2z=`calc "($plane2_z)/($v2mag)" | sed 's/^\t//g'` # normalization normal vector 2
vn1x=(plane1_x)/(v1mag) ; vn1y=(plane1_y)/(v1mag) ; vn1z=(plane1_z)/(v1mag) # normalization normal vector 1
vn2x=(plane2_x)/(v2mag) ; vn2y=(plane2_y)/(v2mag) ; vn2z=(plane2_z)/(v2mag) # normalization normal vector 2
# calc "acos(($vn1x)*($vn2x)+($vn1y)*($vn2y)+($vn1z)*($vn2z))*180/3.141592653589793" | sed 's/^\t//g' | sed 's/^~//g'
print acos((vn1x*vn2x)+(vn1y*vn2y)+(vn1z*vn2z))*180/3.141592653589793
}