Commit 3525df6c by Alexander Neukirchen

### Shifted angle calculation to math directory; modernised.

parent 3abb683b
 ... ... @@ -31,4 +31,5 @@ math/grule.f90 math/lh_tofrom_lm.f90 math/gradYlm.f90 math/analysistests.f90 math/pol_angle.f90 )
 !-------------------------------------------------------------------------------- ! 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. !-------------------------------------------------------------------------------- MODULE m_polangle !----------------------------------------------------------------------------- ! Calculates the polar angles theta and phi of a vector with components ! vx, vy and vz. ! Philipp Kurz 2000-02-08 ! Modernised A.N. 2020 !----------------------------------------------------------------------------- CONTAINS SUBROUTINE pol_angle(vx, vy, vz, theta, phi) USE m_constants, ONLY : pimach IMPLICIT NONE REAL, INTENT(IN) :: vx, vy, vz REAL, INTENT(OUT) :: theta, phi REAL :: eps, r, rho, pi pi = pimach() eps = 1.0e-9 rho = SQRT(vx**2 + vy**2) r = SQRT(vx**2 + vy**2 + vz**2) ! Zero vector: IF ( (r.LT.eps) .AND. (rho.LT.eps) ) THEN theta = 0.0 phi = 0.0 ! v on positive z-axis: ELSE IF ( (rho.LT.eps) .AND. vz.GT.0) THEN theta = 0.0 phi = 0.0 ! v on negative z-axis: ELSE IF ( (rho.LT.eps) .AND. vz.LT.0) THEN theta = pi phi = 0.0 ELSE ! v in xy-plane: IF (ABS(vz).LT.eps) THEN theta = pi/2 ELSE theta = ASIN(rho/r) END IF IF ( vz.LT.0 ) THEN theta = pi - theta END IF ! v in yz-plane IF (ABS(vx).LT.eps) THEN phi = pi/2 ELSE phi = ASIN(ABS(vy)/rho) END IF IF ( vx.LT.0 ) phi = pi - phi IF ( vy.LT.0 ) phi = -phi END IF IF(phi<0) phi=phi+2*pi IF(theta<0) theta=theta + 2*pi IF(theta>pi) THEN theta=theta-pi phi=phi+pi END IF phi=mod(phi,2*pi) END SUBROUTINE pol_angle SUBROUTINE sphericaltocart(r, theta, phi, x, y, z) IMPLICIT NONE REAL, INTENT(IN) :: r, theta, phi REAL, INTENT(OUT) :: x, y, z x=r*SIN(theta)*COS(phi) y=r*SIN(theta)*SIN(phi) z=r*COS(theta) END SUBROUTINE sphericaltocart END MODULE m_polangle
 ... ... @@ -9,7 +9,6 @@ vgen/od_mkgz.f vgen/vacp5_0.f vgen/vacp5_z.f vgen/visp5_0.f vgen/pol_angle.f vgen/visp5_z.f ) set(fleur_F90 \${fleur_F90} ... ...
 !-------------------------------------------------------------------------------- ! 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. !-------------------------------------------------------------------------------- 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-9 vxy = sqrt(vx**2 + vy**2) vxyz = sqrt(vx**2 + vy**2 + vz**2) IF ( (vxyz.LT.eps) .AND. (vxy.LT.eps) ) THEN theta = 0.0 phi = 0.0 ELSE IF ( (vxy.LT.eps) .AND. vz.LT.0) THEN theta = pi phi = 0.0 ELSE IF ( (vxy.LT.eps) .AND. vz.GT.0) 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 IF(phi<0) phi=phi+2*pi IF(theta<0) theta=theta + 2*pi IF(theta>pi) THEN theta=theta-pi phi=phi+pi END IF phi=mod(phi,2*pi) END SUBROUTINE pol_angle SUBROUTINE sphericaltocart(r,theta,phi,x,y,z) IMPLICIT NONE REAL, INTENT(IN) :: r,theta,phi REAL, INTENT(OUT):: x,y,z x=r*SIN(theta)*COS(phi) y=r*SIN(theta)*SIN(phi) z=r*COS(theta) END SUBROUTINE sphericaltocart END MODULE m_polangle
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!