BASH,具有四个点的二面角 [英] BASH, Dihedral angle with four points

查看:98
本文介绍了BASH,具有四个点的二面角的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

点:

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'

我找到了以下相关链接 1 2 3 .

I have found these related links 1, 2 and 3.

参考值为58.7 $^{o}$

我得到的值是:70.62525933704842342761 $^{o}$64.23010091217222985704 $^{o}$

有人知道正确地获得最佳算法是什么?

推荐答案

基于您在此线程中其他地方提炼的shell代码,我也将其转录为awk解决方案.由于人们似乎已经找到了_docd版本的用法,因此我将在末尾添加它.我还提供了一个调试版本(在回复的中间).

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

输出为

58.6892


这里是完整的调试版本.我需要它是因为我遇到了编辑错误,例如将y2=$3更改为y=$3! (这些事情发生了;-/)


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
}


这是转录到awk版本的shell


And here is the transcribed shell to awk version

我强烈推荐Grymoire的Awk教程为您提供帮助了解awk编程范例及其内置变量,例如NR((记录数)).

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
}

这篇关于BASH,具有四个点的二面角的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆