A Simple Fortran Transformation Module

Was cleaning my home laptop and found this fortran transformation module that I wrote long time ago.

Has all the transformation that are generally used. Just posting it in case, I need it again.

module transformation 
implicit none
!

! Code: Sukhbinder
! 31st March 2010
!

real (kind=8) :: angle,agrad
real (kind=8),parameter :: PI = 3.141592653589793
real (kind=8),parameter :: torad = PI/180.0d0

real (kind=8) :: rotz(4,4),rotx(4,4),roty(4,4)
real (kind=8) :: shearxy(4,4),shearxz(4,4),shearyz(4,4)
real (kind=8) :: trans(4,4),scale(4,4)


contains
subroutine setshearxy(a,b)
real (kind=8) :: a,b

  shearxy(1:4,1:4) = 0.0d0
  shearxy(1,3) = a
  shearxy(2,3) = b
  shearxy(1,1) =1.0d0
  shearxy(2,2) =1.0d0
  shearxy(3,3) =1.0d0
  shearxy(4,4) =1.0d0
  
end subroutine

subroutine setshearxz(a,b)
real (kind=8) :: a,b

  shearxz(1:4,1:4) = 0.0d0
  shearxz(1,2) = a
  shearxz(3,2) = b
  shearxz(1,1) =1.0d0
  shearxz(2,2) =1.0d0
  shearxz(3,3) =1.0d0
  shearxz(4,4) =1.0d0
  
end subroutine


subroutine setshearyz(a,b)
real (kind=8) :: a,b

  shearyz(1:4,1:4) = 0.0d0
  shearyz(2,1) = a
  shearyz(3,1) = b
  shearyz(1,1) =1.0d0
  shearyz(2,2) =1.0d0
  shearyz(3,3) =1.0d0
  shearyz(4,4) =1.0d0
  
end subroutine

subroutine setscale(p,q,r)
real (kind=8) :: p,q,r

  scale(1:4,1:4) = 0.0d0
  scale(1,1) = p
  scale(2,2) = q
  scale(3,3) = r
  scale(4,4) = 1.0d0
end subroutine

subroutine settrans(p,q,r)
real (kind=8) :: p,q,r

  trans(1:4,1:4) = 0.0d0
  trans(1,4) = p
  trans(2,4) = q
  trans(3,4) = r
  trans(1,1) =1.0d0
  trans(2,2) =1.0d0
  trans(3,3) =1.0d0
  trans(4,4) =1.0d0
  
end subroutine

subroutine setrotz(ang)
real(kind=8) :: ang

rotz(1:4,1:4)=0.0d0

rotz(1,1) = cos(ang*torad)
rotz(1,2) = -1*sin(ang*torad)
rotz(2,1) = sin(ang*torad)
rotz(2,2) = cos(ang*torad)
rotz(3,3) = 1.0d0
rotz(4,4) = 1.0d0

end subroutine


subroutine setrotx(ang)
real(kind=8) :: ang

rotx(1:4,1:4)=0.0d0

rotx(1,1) = 1.0d0
rotx(2,2) = cos(ang*torad) 
rotx(3,2) = sin(ang*torad)
rotx(2,3) = -1*sin(ang*torad)
rotx(3,3) = cos(ang*torad)
rotx(4,4) = 1.0d0

end subroutine

subroutine setroty(ang)
real(kind=8) :: ang

roty(1:4,1:4)=0.0d0

roty(1,1) = cos(ang*torad) 
roty(3,1) = -1*sin(ang*torad) 
roty(2,2) = 1.0d0
roty(1,3) = sin(ang*torad) 
roty(3,3) = cos(ang*torad)
roty(4,4) = 1.0d0

end subroutine


subroutine getnewcoords(oldcords,rotmat)
real (kind=8)              :: oldcords(:,:),rotmat(4,4)
real (kind=8), allocatable :: newcoords(:,:)
real (kind=8) :: oldone(4),newone(4)

integer :: isize,i

isize=size(oldcords,2)

allocate(newcoords(3,isize))


do i=1,isize
 oldone(1:4)=1.0d0
 oldone(1:3) = oldcords(1:3,i)
 newone(1:4)=1.0d0
 newone=matmul(rotmat,oldone)
 newcoords(1:3,i)=newone(1:3)
end do

oldcords = newcoords
deallocate(newcoords)


end subroutine

end module

! !!!!!!!!!!!!!!!!! example use !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program test
use transformation 

!open(20,file='bladeps.txt',status='old')
!open(40,file='30degree.txt',status='new')
real(kind=8) :: oldcords(4),xyz(3,10),temp(3,10)

call random_number(oldcords)
call random_number(xyz)
oldcords=oldcords*10.0d0
oldcords(4)=1.0d0

xyz=xyz*10.0d0
write(*, '(3f10.2)') ((xyz(i,j),i=1,3),j=1,10)
write(*,*)
write(*, '(4f10.2)') (oldcords(i), i=1,4)
write(*,*)
call setrotz(30.0d0)
write(*, '(4f10.2)') ((rotz(i,j),j=1,4),i=1,4)
temp=xyz
call getnewcoords(temp,rotz)
write(*,*)
write(*, '(3f10.2)') ((temp(i,j),i=1,3),j=1,10)

!close(30)
!close(20)

end program test