Commit bc7c1988 authored by aarograh's avatar aarograh
Browse files

Add more to Geom_Hex 

bound_angle, hex_side_length_from_pitch,
hex_array_element_centroid, rotate_2d_vector,
and orientation enumerations.

Also makes these things public in the Geom module
parent 11d6c80d

src/Geom.f90

+12 −0
Original line number Diff line number Diff line
@@ -20,6 +20,7 @@ USE Geom_Plane 
USE Geom_CircCyl

USE Geom_Box

USE Geom_Poly

USE Geom_Hex

USE ParameterLists



IMPLICIT NONE
@@ -48,6 +49,17 @@ PUBLIC :: OPERATOR(==) 
PUBLIC :: OPERATOR(/=)

PUBLIC :: OPERATOR(.APPROXEQA.)

PUBLIC :: ASSIGNMENT(=)

PUBLIC :: HEX_FLAT_TOP

PUBLIC :: HEX_POINTY_TOP

PUBLIC :: hex_rows

PUBLIC :: hex_columns

PUBLIC :: hex_mid

PUBLIC :: hex_max_columns

PUBLIC :: hex_elements_from_rings

PUBLIC :: bound_angle

PUBLIC :: hex_side_length_from_pitch

PUBLIC :: hex_array_element_centroid

PUBLIC :: rotate_2d_vector



INTERFACE newGeom

  !> @copybrief Geom::newGeom_line

src/Geom_Hex.f90

+140 −0
Original line number Diff line number Diff line
@@ -10,7 +10,10 @@ 
!>

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

MODULE Geom_Hex

#include "Futility_DBC.h"

USE Futility_DBC

USE IntrType

USE Constants_Conversion



IMPLICIT NONE

PRIVATE
@@ -20,6 +23,15 @@ PUBLIC :: hex_columns 
PUBLIC :: hex_mid

PUBLIC :: hex_max_columns

PUBLIC :: hex_elements_from_rings

PUBLIC :: bound_angle

PUBLIC :: hex_side_length_from_pitch

PUBLIC :: hex_array_element_centroid

PUBLIC :: rotate_2d_vector



!> Enumerations for hexagonal unit orientations

INTEGER(SIK),PARAMETER,PUBLIC :: HEX_POINTY_TOP=0

INTEGER(SIK),PARAMETER,PUBLIC :: HEX_FLAT_TOP=1

INTEGER(SIK),PARAMETER,PUBLIC :: HEX_ORIENTATIONS(2)=[HEX_POINTY_TOP, HEX_FLAT_TOP]



!

!===============================================================================
@@ -92,4 +104,132 @@ ELEMENTAL FUNCTION hex_elements_from_rings(rings) RESULT(elements) 


ENDFUNCTION hex_elements_from_rings

!

!-------------------------------------------------------------------------------

!> @brief Bounds an angle to a given range by adding or subtracting 360.0

!> @param a the ange to bound

!> @param lower the lower end of the range

!> @param upper the upper end of the range

!> @param rad indicates if radians should be used; defaults to false

!>

!> The range is inclusive on the lower end and exclusive on the upper end, e.g.,

!> @c [lower,upper).  If the range specified is narrower than 360.0, no checks

!> are performed to ensure the final result is in the range.  The client code

!> is responsible for ensuring this.

!>

ELEMENTAL SUBROUTINE bound_angle(a, lower, upper, rad)

  REAL(SRK),INTENT(INOUT) :: a

  REAL(SRK),INTENT(IN) :: lower

  REAL(SRK),INTENT(IN) :: upper

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: rad

  !

  REAL(SRK) :: increment



  increment = 360.0

  IF(PRESENT(rad)) THEN

    IF(rad) increment = TWOPI

  ENDIF

  DO WHILE(a > upper)

    a = a - increment

  ENDDO

  DO WHILE(a <= lower)

    a = a + increment

  ENDDO



ENDSUBROUTINE bound_angle

!

!-------------------------------------------------------------------------------

!> @brief Calculates the hexagon side length from the pitch

!> @param pitch the pitch of the hexagon

!> @returns length the length of one side of the hexagon

!>

FUNCTION hex_side_length_from_pitch(pitch) RESULT(length)

  REAL(SRK),INTENT(IN) :: pitch

  REAL(SRK) :: length



  length = pitch / SQRT(3.0)



ENDFUNCTION hex_side_length_from_pitch

!

!-------------------------------------------------------------------------------

!> @brief calculates the centroid of a hexagonal array element

!> @param pitch the pitch of the hexagonal array

!> @param row the row index in the array (1 is the top row)

!> @param col the column index in the array (1 is the first element in the row)

!> @param rings the number of rings in the array

!> @param orientation whether the elements are flat or pointy topped

!> @returns coords the coordinates of the centroid

!>

!> If the orientation is flat-topped, then the first "row" is assumed to be at

!> a 30-degree angle relative to the positive x-axis.  In either orientation,

!> the origin (0.0, 0.0) is assumed to be at the center of the hexagon at the

!> center of the array.

!>

FUNCTION hex_array_element_centroid(pitch, row, col, rings, orientation) RESULT(coords)

  REAL(SRK),INTENT(IN) :: pitch

  INTEGER(SIK),INTENT(IN) :: row

  INTEGER(SIK),INTENT(IN) :: col

  INTEGER(SIK),INTENT(IN) :: rings

  INTEGER(SIK),INTENT(IN) :: orientation

  REAL(SRK) :: coords(2)

  !

  INTEGER(SIK) :: mid

  REAL(SRK) :: vec(2)



  REQUIRE(pitch > 0.0)

  REQUIRE(rings > 0)

  REQUIRE(row > 0)

  REQUIRE(row <= hex_rows(rings))

  REQUIRE(col > 0)

  REQUIRE(col <= hex_columns(rings, row))

  REQUIRE(ANY(orientation == HEX_ORIENTATIONS))



  mid = hex_mid(rings)



  ! Calculate x,y for array of flat-topped hexagons first

  vec(1) = (-(mid - 1) + ABS(mid - row) * 0.5 + col) * pitch

  vec(2) = (mid - row) * pitch * COS(SIXTHPI)



  ! Now rotate 30 degrees for pointy-topped

  IF(orientation == HEX_POINTY_TOP) THEN

    coords = rotate_2d_vector(vec, SIXTHPI, rad=.TRUE.)

  ELSE

    coords = vec

  ENDIF



ENDFUNCTION hex_array_element_centroid

!

!-------------------------------------------------------------------------------

!> @brief rotates a 2D vector by some angle

!> @param vec the vector to rotate

!> @param angle the angel by which to rotate the vector

!> @param rad whether to use radians or not; defaults to false

!> @returns rot_vec the rotated vector

!>

!> The @c vec and @c angle must both be degrees (@c rad=false) or both be radians

!> (@c rad=true).

!>

FUNCTION rotate_2d_vector(vec, angle, rad) RESULT(rot_vec)

  REAL(SRK),INTENT(IN) :: vec(2)

  REAL(SRK),INTENT(IN) :: angle

  LOGICAL(SBK),INTENT(IN),OPTIONAL :: rad

  REAL(SRK) :: rot_vec(2)

  !

  LOGICAL(SBK) :: lrad

  REAL(SRK) :: a, lower,upper



  lrad = .FALSE.

  IF(PRESENT(rad)) lrad = rad

  a = angle

  lower = ZERO

  IF(lrad) THEN

    upper = TWOPI

  ELSE

    upper = 360.0

  ENDIF

  CALL bound_angle(a, lower, upper, RAD=lrad)

  rot_vec(1) = COS(a)*vec(1) - SIN(a)*vec(2)

  rot_vec(2) = SIN(a)*vec(1) + COS(a)*vec(2)



ENDFUNCTION rotate_2d_vector

!

ENDMODULE Geom_Hex
 
 No newline at end of file