Lambert.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/Lambert.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/Lambert.f90"
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! Modified      2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are 
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list 
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this 
!        list of conditions and the following disclaimer in the documentation and/or 
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names 
!        of its contributors may be used to endorse or promote products derived from 
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
 
!--------------------------------------------------------------------------
!
!--------------------------------------------------------------------------
module lambert
  use prec
  use math
 
  implicit none
  private
  
  real(preal), parameter :: &
    spi  = sqrt(pi), &
    pref = sqrt(6.0_preal/pi), &
    a    = pi**(5.0_preal/6.0_preal)/6.0_preal**(1.0_preal/6.0_preal), &
    ap   = pi**(2.0_preal/3.0_preal), &
    sc   = a/ap, &
    beta = a/2.0_preal, &
    r1   = (3.0_preal*pi/4.0_preal)**(1.0_preal/3.0_preal), &
    r2   = sqrt(2.0_preal), &
    pi12 = pi/12.0_preal, &
    prek = r1 * 2.0_preal**(1.0_preal/4.0_preal)/beta
 
  public :: &
    lambert_cubetoball, &
    lambert_balltocube
 
contains
 
 
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
pure function lambert_cubetoball(cube) result(ball)
 
  real(preal), intent(in), dimension(3) :: cube
  real(preal),             dimension(3) :: ball, lamxyz, xyz
  real(preal),             dimension(2) :: t
  real(preal)                           :: c, s, q
  real(preal), parameter                :: eps = 1.0e-8_preal
  integer,                 dimension(3) :: p
  integer,                 dimension(2) :: order
  
  if (maxval(abs(cube)) > ap/2.0+eps) then
    ball = ieee_value(cube,ieee_positive_inf)
    return
  end if
  
  ! transform to the sphere grid via the curved square, and intercept the zero point
  center: if (all(deq0(cube))) then
    ball  = 0.0_preal
  else center
    ! get pyramide and scale by grid parameter ratio
    p = getpyramidorder(cube)
    xyz = cube(p) * sc
  
    ! intercept all the points along the z-axis
    special: if (all(deq0(xyz(1:2)))) then
      lamxyz = [ 0.0_preal, 0.0_preal, pref * xyz(3) ]
    else special
      order = merge( [2,1], [1,2], abs(xyz(2)) <= abs(xyz(1)))                                      ! order of absolute values of XYZ
      q = pi12 *  xyz(order(1))/xyz(order(2))                                                       ! smaller by larger
      c = cos(q)
      s = sin(q)
      q = prek * xyz(order(2))/ sqrt(r2-c)
      t = [ (r2*c - 1.0), r2 * s] * q
  
      ! transform to sphere grid (inverse Lambert)
      ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
      c = sum(t**2)
      s = pi * c/(24.0*xyz(3)**2)
      c = spi * c / sqrt(24.0_preal) / xyz(3)
      q = sqrt( 1.0 - s )
      lamxyz = [ t(order(2)) * q, t(order(1)) * q, pref * xyz(3) - c ]
    endif special
  
    ! reverse the coordinates back to order according to the original pyramid number
    ball = lamxyz(p)
 
  endif center
 
end function lambert_cubetoball
 
 
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
pure function lambert_balltocube(xyz) result(cube)
 
  real(preal), intent(in), dimension(3) :: xyz
  real(preal),             dimension(3) :: cube, xyz1, xyz3
  real(preal),             dimension(2) :: tinv, xyz2
  real(preal)                           :: rs, qxy, q2, sq2, q, tt
  integer,                 dimension(3) :: p
  
  rs = norm2(xyz)
  if (rs > r1) then 
    cube = ieee_value(cube,ieee_positive_inf)
    return
  endif
  
  center: if (all(deq0(xyz))) then 
    cube = 0.0_preal
  else center
    p = getpyramidorder(xyz)
    xyz3 = xyz(p)
    
    ! inverse M_3
    xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
    
    ! inverse M_2
    qxy = sum(xyz2**2)
    
    special: if (deq0(qxy)) then 
      tinv = 0.0_preal
    else special
      q2 = qxy + maxval(abs(xyz2))**2
      sq2 = sqrt(q2)
      q = (beta/r2/r1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
      tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/r2/qxy 
      tinv = q * sign(1.0_preal,xyz2) * merge([ 1.0_preal, acos(math_clip(tt,-1.0_preal,1.0_preal))/pi12], &
                                              [ acos(math_clip(tt,-1.0_preal,1.0_preal))/pi12, 1.0_preal], &
                                               abs(xyz2(2)) <= abs(xyz2(1)))
    endif special
    
    ! inverse M_1
    xyz1 = [ tinv(1), tinv(2), sign(1.0_preal,xyz3(3)) * rs / pref ] /sc
    
    ! reverse the coordinates back to order according to the original pyramid number
    cube = xyz1(p)
 
  endif center
 
end function lambert_balltocube
 
 
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
pure function getpyramidorder(xyz)
 
  real(preal),intent(in),dimension(3) :: xyz
  integer,               dimension(3) :: getpyramidorder
 
  if      (((abs(xyz(1)) <=  xyz(3)).and.(abs(xyz(2)) <=  xyz(3))) .or. &
           ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
    getpyramidorder = [1,2,3]
  else if (((abs(xyz(3)) <=  xyz(1)).and.(abs(xyz(2)) <=  xyz(1))) .or. &
           ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
    getpyramidorder = [2,3,1]
  else if (((abs(xyz(1)) <=  xyz(2)).and.(abs(xyz(3)) <=  xyz(2))) .or. &
           ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
    getpyramidorder = [3,1,2]
  else
    getpyramidorder = -1                                                                            ! should be impossible, but might simplify debugging
  end if
 
end function getpyramidorder
 
end module lambert