pol_angle.f 1.91 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
      MODULE m_polangle
c***********************************************************************
c calculates the polar angle theta and phi of a vector with components
c vx, vy and vz.
c Philipp Kurz 2000-02-08
c***********************************************************************
      CONTAINS
      SUBROUTINE pol_angle(
     >                     vx,vy,vz,
     <                     theta,phi)

      USE m_constants, ONLY : pimach
      IMPLICIT NONE

C     .. Scalar Arguments ..
      REAL, INTENT    (IN) :: vx,vy,vz
      REAL, INTENT   (OUT) :: theta,phi
C     ..
C     .. Local Scalars ..
      REAL eps,vxyz,vxy,pi
C     ..

      pi = pimach()
      eps = 1.0e-8

      vxy  = sqrt(vx**2 + vy**2)
      vxyz = sqrt(vx**2 + vy**2 + vz**2)

      IF ( (vxyz.LT.eps) .OR. (vxy.LT.eps) ) THEN
         theta = 0.0
         phi   = 0.0
      ELSE
c--->    due to rounding errors vxy/vxyz can become >1, if vz is very
c--->    small. therefore, make sure that vz is not to small.
         IF (abs(vz).LT.eps) THEN
            theta = pi/2
         ELSE
            theta = asin(vxy/vxyz)
         ENDIF
         IF ( vz.LT.0 ) theta = pi - theta

c--->    due to rounding errors vy/vxy can become >1, if vx is very
c--->    small. therefore, make sure that vx is not to small.
         IF (abs(vx).LT.eps) THEN
            phi = pi/2
         ELSE
            phi   = asin(abs(vy)/vxy)
         ENDIF
         IF ( vx.LT.0 ) phi = pi - phi
         IF ( vy.LT.0 ) phi = -phi
      ENDIF

      END SUBROUTINE pol_angle
      END MODULE m_polangle