Rotations from FreeCAD to GDML file #112
-
|
Hello, I am making a simple detector setup with three ElTubes with their axes concentric with the origin. So when I move them to a cylindrical position which is not the origin, I have to rotate each ElTube so that its axis intersects with the origin. When I do the rotation in FreeCAD, I used the standard Euler Equations to get the following rotation: I have tried to hard code my geometry in the gdml format, with the dimensions I used in FreeCAD, but at extreme angles for theta, the ElTubes were no longer aligned. I also tested the FreeCAD model at different values of theta and it does remain aligned. Thanks for the help! |
Beta Was this translation helpful? Give feedback.
Replies: 3 comments 5 replies
-
|
Did you check if the gdml file produces the correct rotation when displaying the Eltubes? In geant the rotation is specified by a sequence of rotations around the x-axis, then y, then z (or vice-versa I have to check the code again). But in FreeCAD the rotation is given by only one angle around an arbitrary axis. One can always rewrite a rotation about one axis as a sequence of rotations around x, y, and z. So what what you see in FreeCAD (one angle, x, y, z components of the axis) is different than what you see in geant (angle around x, angle around y, angle around z), but the result should be the same. |
Beta Was this translation helpful? Give feedback.
-
|
Sorry for the late reply, but are the rotations in GDML done like euler angles, where the x, y, and z rotations are done with respect to the body's frame, or are the rotations done with respect to a fixed frame that does not change as each x, y, and z rotation is done on the body? I have been trying to get from the FreeCAD rotation to the GDML one and I think I am not interpreting the GDML rotation correctly. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Ian. Did you take a look at my comment in the GitHub forum. I believe I
answered this question there. If you are interested, I can post the code
that I use to get the x, y, z rotation angles from the FreeCAD rotation.
Munther
…On Fri, Jun 9, 2023, 12:34 PM Ian L ***@***.***> wrote:
Sorry for the late reply, but are the rotations in GDML done like euler
angles, where the x, y, and z rotations are done with respect to the body's
frame, or are the rotations done with respect to a fixed frame that does
not change as each x, y, and z rotation is done on the body? I have been
trying to get from the FreeCAD rotation to the GDML one and I think I am
not interpreting the GDML rotation correctly.
—
Reply to this email directly, view it on GitHub
<#112 (reply in thread)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AWV3NCCGXHRHZ5BE6T25LWLXKN3FBANCNFSM6AAAAAAY3NVF6Y>
.
You are receiving this because you commented.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
|
Yes that would be helpful, thank you. My main source of confusion is the which coordinate system the rotation is being done under, whether that be the body's local coordinate system or a fixed lab coordinate system. I apologize for the lack of clarity. Seeing the code to get from FreeCAD rotation to GDML rotation would clear up my confusion. Thank you again for your help! |
Beta Was this translation helpful? Give feedback.
-
|
@IanLuebbers-UIUC , here is the code we use (and which I wrote). FreeCAD had its own conversion from quaternions to rotations about x, y, z, but I found that it did not work all the time. I believe my code works in all cases, but I obviously have not tested the universe of possibilities! The rotations are about the fixed X, Y, Z axes, not body axes. I am sorry, but I forget now whether in geant the translation happens before or after the rotation. Easisest thing to do is to do some experiments and check. def quaternion2XYZ(rot):
"""
convert a quaternion rotation to a sequence of rotations around X, Y, Z
Here is my (Munther Hindi) derivation:
First, the rotation matrices for rotations around X, Y, Z axes.
Rx = [ 1 0 0]
[ 0 cos(a) -sin(a)]
[ 0 sin(a) cos(a)]
Ry= [ cos(b) 0 sin(b)]
[ 0 1 0]
[-sin(b) 0 cos(b)]
Rz = [ cos(g) -sin(g) 0]
[ sin(g) cos(g) 0]
[ 0 0 1]
Rederivation from the previous version. Geant processes the rotation from
the gdml as R = Rz Ry Rx, i.e, Rx applied last, not first, so now we have
R = Rx Ry Rz =
[cosb*cosg, -cosb*sing, sinb],
[cosa*sing+cosg*sina*sinb, cosa*cosg-sina*sinb*sing, -cosb*sina],
[sina*sing-cosa*cosg*sinb, cosa*sinb*sing+cosg*sina, cosa*cosb]
To get the angles a(lpha), b(eta) for rotations around x, y axes, transform the unit vector (0,0,1)
[x,y,z] = Q*(0,0,1) = R*(0,0,1) ==>
x = sin(b)
y = -sin(a)cos(b)
z = cos(a)cos(b)
==> a = atan2(-y, x) = atan2(sin(a)*cos(b), cos(a)*cos(b)) = atan2(sin(a), cos(a))
then b = atan2(x*cos(a), z) = atan2(sin(b)*cos(a), cos(b)*cos(a)] = atan2(sin(b), cos(b))
Once a, b are found, g(amma) can be found by transforming (1, 0, 0), or (0,1,0)
Since now a, b are known, one can form the inverses of Rx, Ry:
Rx^-1 = Rx(-a)
Ry^-1 = Ry(-b)
Now R*(1,0,0) = Rx*Ry*Rz(1,0,0) = (x, y, z)
multiply both sides by Ry^-1 Rx^-1:
Ry^-1 Rx^-1 Rx Ry Rz (1,0,0) = Rz (1,0,0) = Ry(-b) Rx(-a) (x, y, z) = (xp, yp, zp)
==>
xp = cos(g)
yp = sin(g)
zp = 0
and g = atan2(yp, zp)
"""
v = rot * Vector(0, 0, 1)
print(v)
# solution 1.
if v.x > 1.0:
v.x = 1.0
if v.x < -1.0:
v.x = -1.0
b = math.asin(v.x)
if math.cos(b) > 0:
a = math.atan2(-v.y, v.z)
else:
a = math.atan2(v.y, -v.z)
# sanity check 1
ysolution = -math.sin(a) * math.cos(b)
zsolution = math.cos(a) * math.cos(b)
if v.y * ysolution < 0 or v.z * zsolution < 0:
print("Trying second solution")
b = math.pi - b
if math.cos(b) > 0:
a = math.atan2(-v.y, v.z)
else:
a = math.atan2(v.y, -v.z)
# sanity check 2
ysolution = -math.sin(a) * math.cos(b)
zsolution = math.cos(a) * math.cos(b)
if v.y * ysolution < 0 or v.z * zsolution < 0:
print("Failed both solutions!")
print(v.y, ysolution)
print(v.z, zsolution)
Ryinv = FreeCAD.Rotation(Vector(0, 1, 0), math.degrees(-b))
Rxinv = FreeCAD.Rotation(Vector(1, 0, 0), math.degrees(-a))
vp = Ryinv * Rxinv * rot * Vector(1, 0, 0)
g = math.atan2(vp.y, vp.x)
return [math.degrees(a), math.degrees(b), math.degrees(g)]
|
Beta Was this translation helpful? Give feedback.
-
|
This is very helpful, Thank you again! |
Beta Was this translation helpful? Give feedback.
-
|
Good! Please let me know if you questions. |
Beta Was this translation helpful? Give feedback.
-
|
Note the comments (left over from the first iteration), say a = atan2(-y, x), but the correct value is atan2(-y, z), which is in the code. |
Beta Was this translation helpful? Give feedback.
@IanLuebbers-UIUC , here is the code we use (and which I wrote). FreeCAD had its own conversion from quaternions to rotations about x, y, z, but I found that it did not work all the time. I believe my code works in all cases, but I obviously have not tested the universe of possibilities! The rotations are about the fixed X, Y, Z axes, not body axes. I am sorry, but I forget now whether in geant the translation happens before or after the rotation. Easisest thing to do is to do some experiments and check.