Databases
StackOverflow2013
Postsdbo
POerror creating composite rotation matrices
Body
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POerror creating composite rotation matrices
    text
    copied!<p>I had first posted this question on Math.Stackexchange, but the guys there suggested this was the right place for the problem.</p> <p>I'm trying to rotate a set of points on a sphere along more than one axis, to get a specific orientation. I wrote code for the same in Fortran and the transformation seems to work well if I rotate the points in any one direction. But as soon as I specify more than one axis of rotation, the solution goes awry. Most of the points, save maybe the extreme ones, go out of the sphere's volume.</p> <p>My code (pertinent to the question) goes as follows:</p> <pre><code>DO iSlice = 1, nSlices IF(iSlice&lt;10) THEN WRITE(string1,'("slice",i1.1,".dat")'), iSlice ELSE WRITE(string1,'("slice",i2.2,".dat")'), iSlice END IF !PRINT*,string1 OPEN (1, file = string1) DO j = 1,4 READ(1,*) z(j),y(j),t(j) END DO READ(1,*) temp x(1:4) = -1.*temp CLOSE (1) !Transform the nodal positions as desired CALL init_Trans CALL form_translate(-1.*CG) IF (iRotz) THEN ang = (the3*pi/180.0)*(nSlices-iSlice)/(nSlices-1) r_axis = '-rz' PRINT*,'angz = ',ang*180./pi CALL form_rotate(r_axis,ang) END IF IF (iRoty) THEN ang = (the2*pi/180.0)*(nSlices-iSlice)/(nSlices-1) r_axis = '-ry' PRINT*,'angy = ',ang*180./pi CALL form_rotate(r_axis,ang) END IF IF (iRotx) THEN ang = (the1*pi/180.0)*(nSlices-iSlice)/(nSlices-1) r_axis = '-rx' PRINT*,'angx = ',ang*180./pi CALL form_rotate(r_axis,ang) END IF CALL form_translate(CG) CALL transform WRITE(fout,'(A,I2)')'SLICE: ',iSlice DO i = 1, 4 WRITE(fout,71) xnew(i), ynew(i), znew(i), t(i) END DO WRITE(fout,*)'' DO i = 1, 4 WRITE(fout_tec,'(3(D25.17,1X))') xnew(i), ynew(i), znew(i) END DO END DO </code></pre> <p>Where the subroutine definitions are as follows:</p> <pre><code> SUBROUTINE init_Trans USE global_param, ONLY: Trans IMPLICIT NONE INTEGER :: i,j ! This subroutine initializes the final Transformation matrix as a 4x4 identity matrix !--------------------------------------------------------------------------------------! Trans (:,:) = 0 DO i = 1, 4 Trans(i,i) = 1 END DO END SUBROUTINE init_Trans ! -------------------------------------------------------------------------------------------- ! SUBROUTINE form_translate (D_MOVE) USE global_param, ONLY: Trans IMPLICIT NONE INTERFACE SUBROUTINE mat_prod(M1, M2, M3) REAL, DIMENSION(:,:), INTENT(IN) :: M1, M2 REAL, DIMENSION(:,:), INTENT(OUT) :: M3 END SUBROUTINE END INTERFACE REAL, DIMENSION(3), INTENT(IN) :: D_MOVE REAL, DIMENSION(4,4) :: T, Temp INTEGER :: i, j ! Initialize Translation Matrix !-------------------------------! ! WRITE(*,*)'FORMULATING TRANSLATION MATRIX...' DO i = 1, 3 DO j = 1, 3 IF (i .eq. j) THEN T(i,j) = 1.0; ELSE T(i,j) = 0.0; END IF END DO T(i,4) = D_MOVE(i) END DO T(4,4) = 1.0; CALL mat_prod(T,Trans,Temp) Trans = Temp END SUBROUTINE form_translate ! -------------------------------------------------------------------------------------------- ! SUBROUTINE form_rotate(r_axis, theta) USE global_param, ONLY: Trans IMPLICIT NONE INTERFACE SUBROUTINE mat_prod(M1, M2, M3) REAL, DIMENSION(:,:), INTENT(IN) :: M1, M2 REAL, DIMENSION(:,:), INTENT(OUT) :: M3 END SUBROUTINE END INTERFACE REAL, INTENT(IN) :: theta REAL, DIMENSION(4,4) :: R, Temp CHARACTER(LEN = 3), INTENT(IN) :: r_axis INTEGER :: i,j SELECT CASE (r_axis) CASE('-rx') ! WRITE(*,*)'FORMULATING X-ROTATION MATRIX...' ! Initialize x-Rotation matrix !----------------------------! DO i = 2,3 DO j = 2, 3 IF (i .eq. j) THEN R(i,j) = cos(theta) ELSE R(i,j) = ((-1)**(i-1))*sin(theta) END IF END DO END DO DO i = 1, 4, 3 DO j = 1, 4, 3 IF (i .eq. j) R(i,j) = 1.0 END DO END DO CASE('-ry') ! WRITE(*,*)'FORMULATING Y-ROTATION MATRIX...' ! Initialize y-Rotation matrix !----------------------------! DO i = 1,3,2 R(i,i) = cos(theta) IF (mod(i+1,4) .eq. 0) THEN R(i,4-i) = -sin(theta) ELSE R(i,4-i) = sin(theta) END IF END DO DO i = 2, 4, 2 R(i,i) = 1.0 END DO CASE('-rz') ! WRITE(*,*)'FORMULATING Z-ROTATION MATRIX...' ! Initialize z-Rotation matrix !----------------------------! DO i = 1,2 DO j = 1, 2 IF (i .eq. j) THEN R(i,j) = cos(theta) ELSE R(i,j) = ((-1)**i)*sin(theta) END IF END DO END DO DO i = 3, 4 DO j = 3, 4 IF (i .eq. j) R(i,j) = 1.0 END DO END DO END SELECT CALL mat_prod(R,Trans,Temp) Trans = Temp DO i = 1, 4 PRINT*,Trans(i,:) END DO END SUBROUTINE form_rotate ! -------------------------------------------------------------------------------------------- ! SUBROUTINE transform ! Transform coordinates of each point using p_new = Tran*p !----------------------------------------------------------! USE global_param IMPLICIT NONE INTERFACE SUBROUTINE mat_prod(M1, M2, M3) REAL, DIMENSION(:,:), INTENT(IN) :: M1, M2 REAL, DIMENSION(:,:), INTENT(OUT) :: M3 END SUBROUTINE END INTERFACE REAL, DIMENSION(4,1) :: p, p_new INTEGER :: i,j p(4,1) = 1.0 DO i = 1, n_nodes p(1,1) = x(i); p(2,1) = y(i); p(3,1) = z(i) CALL mat_prod(Trans, p, p_new) xnew(i) = p_new(1,1); ynew(i) = p_new(2,1); znew(i) = p_new(3,1) END DO END SUBROUTINE transform </code></pre> <p>I have tried creating a composite rotation matrix so that the resulting rotation, along both the axes, may be handled at once. And I guess this is where I'm going wrong, as the code works well for rotating points about individual axes. Could anyone help me figure out the problem?</p> <p>Thanks a ton!</p>
 

Querying!

 
Guidance
An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload