BASH,四点二面角
BASH, Dihedral angle with four 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
到目前为止,我已经完成了:
#!/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 }' LINEA` # x1
y1=`awk '{print }' LINEA` # y1
z1=`awk '{print }' LINEA` # z1
x2=`awk '{print }' LINEB` # x2
y2=`awk '{print }' LINEB` # y2
z2=`awk '{print }' LINEB` # z2
x3=`awk '{print }' LINEC` # x3
y3=`awk '{print }' LINEC` # y3
z3=`awk '{print }' LINEC` # z3
x4=`awk '{print }' LINED` # x4
y4=`awk '{print }' LINED` # y4
z4=`awk '{print }' 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'
我找到了这些相关链接 1, 2 and 3。
参考值为58.7 $^{o}$
我得到的值是:70.62525933704842342761 $^{o}$
和64.23010091217222985704 $^{o}$
有人知道正确获取它的最佳算法是什么吗?
编辑: 如果您在互联网上搜索 torsion.awk
已将您带到这里,请跳到上面已接受的答案,因为它使用 O.P.s改进算法来计算扭转,但仍然演示将 shell 代码转换为 awk
.
以前的读者,还要注意在下面的第 2 个 edit 中使用此代码的改进。
我刚刚注意到最后的 "properly" 资格 ;-/
这是您转换为一个 awk 进程的代码。
我对这个水平的数学没有经验,所以不能说它真的是在计算你需要的结果。
此外,awk 程序中经常存在关于精度的问题,这实际上与编译的基础 c
语言库有关。
无论如何,考虑到所有的注意事项,这里是您的代码的基本转换。
猫torsion_docd.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
# x1=`awk '{print }' LINEA` # x1
# y1=`awk '{print }' LINEA` # y1
# z1=`awk '{print }' LINEA` # z1
# x2=`awk '{print }' LINEB` # x2
# y2=`awk '{print }' LINEB` # y2
# z2=`awk '{print }' LINEB` # z2
# x3=`awk '{print }' LINEC` # x3
# y3=`awk '{print }' LINEC` # y3
# z3=`awk '{print }' LINEC` # z3
# x4=`awk '{print }' LINED` # x4
# y4=`awk '{print }' LINED` # y4
# z4=`awk '{print }' LINED` # z4
NR==1 {x1=; y=; z1=}
NR==2 {x2=; y=; z2=}
NR==3 {x3=; y=; z3=}
NR==4 {
x4=; y=; z4=
# all of this code below is only executed when you read in the 4th line
# becuase then you have all the data
# 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'`
v1x=x1-x2 ; v1y=y1-y2 ; v1z=z1-z2
v2x=x4-x3 ; v2y=y4-y3 ; v2z=z4-z3
# 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'`
v1mag=sqrt((v1x)**2+(v1y)**2+(v1z)**2)
v2mag=sqrt((v2x)**2+(v2y)**2+(v2z)**2)
# calc "acos((($v1x)/($v1mag))*(($v2x)/($v2mag))+(($v1y)/($v1mag))*(($v2y)/($v2mag))+(($v1z)/($v1mag))*(($v2z)/($v2mag)))*180/3.141
592653589793" | 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.14159
2653589793" | sed 's/^\t//g' | sed 's/^~//g'
print acos(((v1x)/(v1mag))*((v2x)/(v2mag))+((v1y)/(v1mag))*((v2y)/(v2mag))+((v1z)/(v1mag))*((v2z)/(v2mag)))*180/3.141592653589793
print acos(((x1)*(x4)+(y1)*(y4)+(z1)*(z4))/(sqrt((x1)**2+(y1)**2+(z1)**2)*sqrt((x4)**2+(y4)**2+(z4)**2)))*180/3.141592653589793
}
如果没有转换文档,它看起来像
猫torsion.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
NR==1 {x1=; y=; z1=}
NR==2 {x2=; y=; z2=}
NR==3 {x3=; y=; z3=}
NR==4 {
x4=; y=; z4=
# all of this code below is only executed when you read in the 4th line
# because then you have all the data
v1x=x1-x2 ; v1y=y1-y2 ; v1z=z1-z2
v2x=x4-x3 ; v2y=y4-y3 ; v2z=z4-z3
v1mag=sqrt((v1x)**2+(v1y)**2+(v1z)**2)
v2mag=sqrt((v2x)**2+(v2y)**2+(v2z)**2)
print acos(((v1x)/(v1mag))*((v2x)/(v2mag))+((v1y)/(v1mag))*((v2y)/(v2mag))+((v1z)/(v1mag))*((v2z)/(v2mag)))*180/3.141592653589793
print acos(((x1)*(x4)+(y1)*(y4)+(z1)*(z4))/(sqrt((x1)**2+(y1)**2+(z1)**2)*sqrt((x4)**2+(y4)**2+(z4)**2)))*180/3.141592653589793
}
请注意,我在你的最后两行前面添加了打印语句 acos
。
在我的机器上,我运行它是
awk -f torsion.awk data.txt
EDIT :我已将 #!/bin/awk
固定在脚本的顶部。因此,您需要使用
将脚本标记为可执行文件
chmod +x ./torsion.awk
然后你可以运行它就像
`./torsion.awk data.txt
您的系统可能需要与顶部的 she-bang 行 (#!/bin/awk
) 不同的 awk
路径。键入 which awk
,然后在 #!
之后使用该值。此外,遗留 Unix 实现通常安装了其他版本的 awk
,因此如果这是您的操作环境,请进行一些研究以找出哪个是您系统上最好的 awk
(通常是 gawk
).
# -------------- end edit --------------------
输出
87.6318
131.872
但鉴于您同意 -58.7
是您想要的输出,我将把如何使用 2 acos
计算留给您。
无论如何,希望您能看到使用 awk
进行此类计算要简单得多。
当然,希望真正的数学家参与进来(大笑之后)并帮助纠正这个问题(或提供他们自己的想法)。
IHTH
经过漫长的一夜,我找到了解决方案:
awk -v var=$((x+2)) 'NR==var' $FILE > LINEaa
awk -v var=$((y+2)) 'NR==var' $FILE > LINEbb
awk -v var=$((z+2)) 'NR==var' $FILE > LINEcc
awk -v var=$((w+2)) 'NR==var' $FILE > LINEd
x1=`awk '{print }' LINEaa` # x1
y1=`awk '{print }' LINEaa` # y1
z1=`awk '{print }' LINEaa` # z1
x2=`awk '{print }' LINEbb` # x2
y2=`awk '{print }' LINEbb` # y2
z2=`awk '{print }' LINEbb` # z2
x3=`awk '{print }' LINEcc` # x3
y3=`awk '{print }' LINEcc` # y3
z3=`awk '{print }' LINEcc` # z3
x4=`awk '{print }' LINEd` # x4
y4=`awk '{print }' LINEd` # y4
z4=`awk '{print }' LINEd` # z4
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
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
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
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
calc "acos(($vn1x)*($vn2x)+($vn1y)*($vn2y)+($vn1z)*($vn2z))*180/3.141592653589793" | sed 's/^\t//g' | sed 's/^~//g'
根据您在本线程其他地方精炼的 shell 代码,我也将其转录为 awk
解决方案。由于人们似乎已经找到了 _docd
版本的用法,所以我将在最后包含它。我还包括一个调试版本(在回复中间)。
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=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y4=; 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
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
}
保存文件后,必须将脚本标记为可执行文件:
chmod +x ./torsion2.awk
然后您可以 运行 使用提供的示例数据:
./torsion2.awk data.txt
输出为
58.6892
这是完整的调试版本。我需要它,因为我有编辑错误,比如将 y2=
更改为 y=
! (这些事情发生了 ;-/ )
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=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y4=; z4=
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
}
这里是转录 shell 到 awk 版本
我强烈推荐 Grymoire's Awk Tutorial 来帮助您理解 awk
编程范式及其内置变量,例如 NR
(记录数)。
cat torsion2_docd.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
# x1=`awk '{print }' LINEA` # x1
# y1=`awk '{print }' LINEA` # y1
# z1=`awk '{print }' LINEA` # z1
# x2=`awk '{print }' LINEB` # x2
# y2=`awk '{print }' LINEB` # y2
# z2=`awk '{print }' LINEB` # z2
# x3=`awk '{print }' LINEC` # x3
# y3=`awk '{print }' LINEC` # y3
# z3=`awk '{print }' LINEC` # z3
# x4=`awk '{print }' LINED` # x4
# y4=`awk '{print }' LINED` # y4
# z4=`awk '{print }' LINED` # z4
NR==1 {x1=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y=; z4=
# 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
}
积分:
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
到目前为止,我已经完成了:
#!/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 }' LINEA` # x1
y1=`awk '{print }' LINEA` # y1
z1=`awk '{print }' LINEA` # z1
x2=`awk '{print }' LINEB` # x2
y2=`awk '{print }' LINEB` # y2
z2=`awk '{print }' LINEB` # z2
x3=`awk '{print }' LINEC` # x3
y3=`awk '{print }' LINEC` # y3
z3=`awk '{print }' LINEC` # z3
x4=`awk '{print }' LINED` # x4
y4=`awk '{print }' LINED` # y4
z4=`awk '{print }' 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'
我找到了这些相关链接 1, 2 and 3。
参考值为58.7 $^{o}$
我得到的值是:70.62525933704842342761 $^{o}$
和64.23010091217222985704 $^{o}$
有人知道正确获取它的最佳算法是什么吗?
编辑: 如果您在互联网上搜索 torsion.awk
已将您带到这里,请跳到上面已接受的答案,因为它使用 O.P.s改进算法来计算扭转,但仍然演示将 shell 代码转换为 awk
.
以前的读者,还要注意在下面的第 2 个 edit 中使用此代码的改进。
我刚刚注意到最后的 "properly" 资格 ;-/
这是您转换为一个 awk 进程的代码。
我对这个水平的数学没有经验,所以不能说它真的是在计算你需要的结果。
此外,awk 程序中经常存在关于精度的问题,这实际上与编译的基础 c
语言库有关。
无论如何,考虑到所有的注意事项,这里是您的代码的基本转换。
猫torsion_docd.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
# x1=`awk '{print }' LINEA` # x1
# y1=`awk '{print }' LINEA` # y1
# z1=`awk '{print }' LINEA` # z1
# x2=`awk '{print }' LINEB` # x2
# y2=`awk '{print }' LINEB` # y2
# z2=`awk '{print }' LINEB` # z2
# x3=`awk '{print }' LINEC` # x3
# y3=`awk '{print }' LINEC` # y3
# z3=`awk '{print }' LINEC` # z3
# x4=`awk '{print }' LINED` # x4
# y4=`awk '{print }' LINED` # y4
# z4=`awk '{print }' LINED` # z4
NR==1 {x1=; y=; z1=}
NR==2 {x2=; y=; z2=}
NR==3 {x3=; y=; z3=}
NR==4 {
x4=; y=; z4=
# all of this code below is only executed when you read in the 4th line
# becuase then you have all the data
# 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'`
v1x=x1-x2 ; v1y=y1-y2 ; v1z=z1-z2
v2x=x4-x3 ; v2y=y4-y3 ; v2z=z4-z3
# 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'`
v1mag=sqrt((v1x)**2+(v1y)**2+(v1z)**2)
v2mag=sqrt((v2x)**2+(v2y)**2+(v2z)**2)
# calc "acos((($v1x)/($v1mag))*(($v2x)/($v2mag))+(($v1y)/($v1mag))*(($v2y)/($v2mag))+(($v1z)/($v1mag))*(($v2z)/($v2mag)))*180/3.141
592653589793" | 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.14159
2653589793" | sed 's/^\t//g' | sed 's/^~//g'
print acos(((v1x)/(v1mag))*((v2x)/(v2mag))+((v1y)/(v1mag))*((v2y)/(v2mag))+((v1z)/(v1mag))*((v2z)/(v2mag)))*180/3.141592653589793
print acos(((x1)*(x4)+(y1)*(y4)+(z1)*(z4))/(sqrt((x1)**2+(y1)**2+(z1)**2)*sqrt((x4)**2+(y4)**2+(z4)**2)))*180/3.141592653589793
}
如果没有转换文档,它看起来像
猫torsion.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
NR==1 {x1=; y=; z1=}
NR==2 {x2=; y=; z2=}
NR==3 {x3=; y=; z3=}
NR==4 {
x4=; y=; z4=
# all of this code below is only executed when you read in the 4th line
# because then you have all the data
v1x=x1-x2 ; v1y=y1-y2 ; v1z=z1-z2
v2x=x4-x3 ; v2y=y4-y3 ; v2z=z4-z3
v1mag=sqrt((v1x)**2+(v1y)**2+(v1z)**2)
v2mag=sqrt((v2x)**2+(v2y)**2+(v2z)**2)
print acos(((v1x)/(v1mag))*((v2x)/(v2mag))+((v1y)/(v1mag))*((v2y)/(v2mag))+((v1z)/(v1mag))*((v2z)/(v2mag)))*180/3.141592653589793
print acos(((x1)*(x4)+(y1)*(y4)+(z1)*(z4))/(sqrt((x1)**2+(y1)**2+(z1)**2)*sqrt((x4)**2+(y4)**2+(z4)**2)))*180/3.141592653589793
}
请注意,我在你的最后两行前面添加了打印语句 acos
。
在我的机器上,我运行它是
awk -f torsion.awk data.txt
EDIT :我已将 #!/bin/awk
固定在脚本的顶部。因此,您需要使用
chmod +x ./torsion.awk
然后你可以运行它就像
`./torsion.awk data.txt
您的系统可能需要与顶部的 she-bang 行 (#!/bin/awk
) 不同的 awk
路径。键入 which awk
,然后在 #!
之后使用该值。此外,遗留 Unix 实现通常安装了其他版本的 awk
,因此如果这是您的操作环境,请进行一些研究以找出哪个是您系统上最好的 awk
(通常是 gawk
).
# -------------- end edit --------------------
输出
87.6318
131.872
但鉴于您同意 -58.7
是您想要的输出,我将把如何使用 2 acos
计算留给您。
无论如何,希望您能看到使用 awk
进行此类计算要简单得多。
当然,希望真正的数学家参与进来(大笑之后)并帮助纠正这个问题(或提供他们自己的想法)。
IHTH
经过漫长的一夜,我找到了解决方案:
awk -v var=$((x+2)) 'NR==var' $FILE > LINEaa
awk -v var=$((y+2)) 'NR==var' $FILE > LINEbb
awk -v var=$((z+2)) 'NR==var' $FILE > LINEcc
awk -v var=$((w+2)) 'NR==var' $FILE > LINEd
x1=`awk '{print }' LINEaa` # x1
y1=`awk '{print }' LINEaa` # y1
z1=`awk '{print }' LINEaa` # z1
x2=`awk '{print }' LINEbb` # x2
y2=`awk '{print }' LINEbb` # y2
z2=`awk '{print }' LINEbb` # z2
x3=`awk '{print }' LINEcc` # x3
y3=`awk '{print }' LINEcc` # y3
z3=`awk '{print }' LINEcc` # z3
x4=`awk '{print }' LINEd` # x4
y4=`awk '{print }' LINEd` # y4
z4=`awk '{print }' LINEd` # z4
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
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
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
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
calc "acos(($vn1x)*($vn2x)+($vn1y)*($vn2y)+($vn1z)*($vn2z))*180/3.141592653589793" | sed 's/^\t//g' | sed 's/^~//g'
根据您在本线程其他地方精炼的 shell 代码,我也将其转录为 awk
解决方案。由于人们似乎已经找到了 _docd
版本的用法,所以我将在最后包含它。我还包括一个调试版本(在回复中间)。
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=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y4=; 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
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
}
保存文件后,必须将脚本标记为可执行文件:
chmod +x ./torsion2.awk
然后您可以 运行 使用提供的示例数据:
./torsion2.awk data.txt
输出为
58.6892
这是完整的调试版本。我需要它,因为我有编辑错误,比如将 y2=
更改为 y=
! (这些事情发生了 ;-/ )
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=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y4=; z4=
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
}
这里是转录 shell 到 awk 版本
我强烈推荐 Grymoire's Awk Tutorial 来帮助您理解 awk
编程范式及其内置变量,例如 NR
(记录数)。
cat torsion2_docd.awk
#!/bin/awk -f
function acos(x) { return atan2((1.-x^2)^0.5,x) }
# x1=`awk '{print }' LINEA` # x1
# y1=`awk '{print }' LINEA` # y1
# z1=`awk '{print }' LINEA` # z1
# x2=`awk '{print }' LINEB` # x2
# y2=`awk '{print }' LINEB` # y2
# z2=`awk '{print }' LINEB` # z2
# x3=`awk '{print }' LINEC` # x3
# y3=`awk '{print }' LINEC` # y3
# z3=`awk '{print }' LINEC` # z3
# x4=`awk '{print }' LINED` # x4
# y4=`awk '{print }' LINED` # y4
# z4=`awk '{print }' LINED` # z4
NR==1 {x1=; y1=; z1=}
NR==2 {x2=; y2=; z2=}
NR==3 {x3=; y3=; z3=}
NR==4 {
x4=; y=; z4=
# 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
}