From 9153648f1b773a1c1186d0a0ac29e78b659cd1f8 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Tue, 9 Feb 2021 09:13:12 +0100
Subject: [PATCH] avoid atan2(0.0, 0.0) in Dihedral::SetAngleICS
---
modules/mol/base/src/impl/dihedral.cc | 41 ++++++++++++++-------------
1 file changed, 22 insertions(+), 19 deletions(-)
diff --git a/modules/mol/base/src/impl/dihedral.cc b/modules/mol/base/src/impl/dihedral.cc
index ee26230be..2c8c003c4 100644
--- a/modules/mol/base/src/impl/dihedral.cc
+++ b/modules/mol/base/src/impl/dihedral.cc
@@ -125,26 +125,29 @@ void Dihedral::SetAngleICS(Real angle, bool update_other) {
// not all connectors are living in the same coordinate system. We first
// have to bring all of them into the local coordinate system of a3.
if (a2->GetPrimaryConnector() && a2->GetPrimaryConnector()->GetFirst()==a1) {
- v1=geom::Vec3(0, 0,-1)*conn2->GetLocalRot();
- v2=geom::Vec3(0, 0, 1);
- ConnectorImplP c=GetConnector(a4, a3);
- v3=c->GetDir();
- Real phi1=atan2(v1[1], v1[0]);
- geom::Vec3 n3=vec_for_angle(v3, angle+phi1);
- c->SetDir(n3);
- // check if we have to update the other connectors
- const ConnectorImplList& sc=a3->GetSecondaryConnectors();
- if (update_other && sc.size()>1) {
- Real phi2=(angle-phi1)-atan2(v3[1], v3[0]);
- for (ConnectorImplList::const_iterator i=sc.begin(); i!=sc.end(); ++i) {
- if (*i==c)
- continue;
- geom::Vec3 v=(*i)->GetDir();
- Real phi4=phi2+atan2(v[1], v[0]);
- geom::Vec3 n=vec_for_angle(v, phi4);
- (*i)->SetDir(n);
- }
+ // equivalent to geom::Transpose(conn2->GetLocalRot()) * geom::Vec3(0,0,-1)
+ v1=geom::Vec3(0, 0,-1)*conn2->GetLocalRot();
+ ConnectorImplP c=GetConnector(a4, a3);
+ v2=c->GetDir();
+ Real phi1 = 0.0;
+ if(std::abs(v1[1]) + std::abs(v1[0]) > Real(1e-6)) phi1=atan2(v1[1], v1[0]);
+ geom::Vec3 n2=vec_for_angle(v2, angle+phi1);
+ c->SetDir(n2);
+ // check if we have to update the other connectors
+ const ConnectorImplList& sc=a3->GetSecondaryConnectors();
+ if (update_other && sc.size()>1) {
+ Real phi2=(angle-phi1);
+ if(std::abs(v2[1]) + std::abs(v2[0]) > Real(1e-6)) phi2-=atan2(v2[1], v2[0]);
+ for (ConnectorImplList::const_iterator i=sc.begin(); i!=sc.end(); ++i) {
+ if (*i==c)
+ continue;
+ v3=(*i)->GetDir();
+ Real phi3=phi2;
+ if(std::abs(v3[1]) + std::abs(v3[0]) > Real(1e-6)) phi3+=atan2(v3[1], v3[0]);
+ geom::Vec3 n=vec_for_angle(v3, phi3);
+ (*i)->SetDir(n);
}
+ }
}
}
--
GitLab