Shape Object

 TYPE ShapeDef
      INTEGER :: Type = 1      
      REAL(KIND=qPREC) :: Position(3)=(/0d0,0d0,0d0/)
      REAL(KIND=qPREC) :: psi=0   !Initial rotation about z axis                                                                                              
      REAL(KIND=qPREC) :: theta=0 !Angle rotated around y (from z towards x)                                                                                  
      REAL(KIND=qPREC) :: phi=0   !Angle rotated around z (from x towards y)                                                                                  
      REAL(KIND=qPREC) :: size_param(3) = (/1d0,1d0,1d0/)
      REAL(KIND=qPREC) :: RotationMatrix(3,3)=RESHAPE((/1d0,0d0,0d0,0d0,1d0,0d0,0d0,0d0,1d0/),(/3,3/))    !Rotates lab frame into object frame                
      REAL(KIND=qPREC) :: xBounds(3,2)=RESHAPE((/-1d0,-1d0,-1d0,1d0,1d0,1d0/),(/3,2/))
   END TYPE ShapeDef

Where type is one of

• SPHERE
• ELLIPSOID
• CYLINDER
• ELLIPTICAL_PRISM
• RECTANGULAR_PRISM

To initialize a shape object you need to do the following in the order below

 CALL SetShapeType(Shape, type, d(3), a, b, c)
 CALL SetShapeOrientation(Shape, psi, theta, phi) !creates rotation matrix
 Shape%Position=position
 CALL SetShapeBounds(Shape) !Sets shape physical bounds

The call to SetShapeType has optional parameters depending on the type. You can always pass an array 'd(' that has the radius or 'half-width' along the x, y, and z directions. Or if the shape type is sphere you can pass one parameter - just the radius. If the shape type is cylinder you can pass the radius and half-height. For all other types you have to pass all three extents in an array or as three separate parameters.

Once a shape is fully specified, you can get the physical bounds of a box that encompasses the shape in Shape%xBounds. You can also check to see if a cell is within a shape and get the rotated coordinates by calling

IsInShape(Shape, pos, coords)

You can also transform a vector in the Shape's frame back to the lab frame by calling

RotateVectorFromShape(Shape, vec)