Commit 2609e6f1 authored by Alexander Neukirchen's avatar Alexander Neukirchen

Appearantly found the reason for exploding potentials.

parent b00dcb9c
......@@ -100,7 +100,7 @@ CONTAINS
TYPE(t_coreSpecInput) :: coreSpecInput
TYPE(t_wann) :: wann
TYPE(t_potden) :: vTot, vx, vCoul, vTemp, vxcForPlotting
TYPE(t_potden) :: inDen, outDen, EnergyDen
TYPE(t_potden) :: inDen, outDen, EnergyDen, dummyDen
TYPE(t_potden), DIMENSION(3) :: testDen
CLASS(t_xcpot), ALLOCATABLE :: xcpot
......@@ -517,8 +517,10 @@ CONTAINS
! Test: Build a field, for which the theoretical divergence etc. are known and
! compare with the result of the routine.
CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeBxc(stars,atoms,sphhar,vacuum,input,noco,vTot,testDen)
CALL dummyDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL matrixsplit(stars, atoms, sphhar, vacuum, input, noco, 1.0, inDen, dummyDen, testDen(1), testDen(2), testDen(3))
CALL checkplotinp()
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', testDen(1), testDen(1), testDen(2), testDen(3))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDeny ', testDen(2))
......
......@@ -200,16 +200,16 @@ CONTAINS
! Calculate the charge and magnetization densities in the muffin tins.
DO ityp = 1,atoms%ntype
theta = noco%beta(ityp)
phi = noco%alph(ityp)
DO ilh = 0,sphhar%nlh(atoms%ntypsy(ityp))
DO iri = 1,atoms%jri(ityp)
IF (SIZE(denmat%mt,4).LE.2) THEN
cdnup = rho(iri,ilh,ityp,1)
cdndown = rho(iri,ilh,ityp,2)
theta = noco%beta(ityp)
phi = noco%alph(ityp)
chden = cdnup + cdndown
mgden = cdnup - cdndown
rho(iri,ilh,ityp,1) = chden
rho(iri,ilh,ityp,2) = mgden*COS(phi)*SIN(theta)
rho(iri,ilh,ityp,3) = mgden*SIN(phi)*SIN(theta)
......@@ -805,8 +805,353 @@ CONTAINS
END SUBROUTINE savxsf
SUBROUTINE savcsv
SUBROUTINE savcsv(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, score, potnorm, denName, denf, denA1, denA2, denA3)
USE m_outcdn
USE m_xsf_io
! Takes one/several t_potden variable(s), i.e. scalar fields in MT-sphere/
! plane wave representation and makes it/them into plottable .xsf file(s)
! according to a scheme given in plot_inp.
!
! This is an adaptation of the old plotdop.f90 by P. Kurz.
!
! Naming convention:
! The output .xsf files will have names composed of the plotted density
! (c.f. identifier) and an optional tag. For a scalar density, only said
! density is plotted. For a spin-polarized density with an up and down
! component, f denotes the sum of both densities and g denotes their
! difference u-d possibly with additional modifying factors as explained
! above in vector-/matrixsplit. f and g are to be understood as scalar
! fields f(r_vec) and g(r_vec). For a density matrix, f is still the
! scalar part. Aditionally, three vector components A_[1-3] arise.
!
! E.g.: scalar density rho (n(r_vec))
! ---> rho.xsf (n(r_vec))
! spin-polarized density rho (n_up(r_vec),n_down(r_vec))
! ---> rho_f.xsf, rho_g.xsf [n(r_vec),m(r_vec)]
! matrix density rho ((n_11(r_vec),n_12(r_vec)),
! (n_21(r_vec),n_22(r_vec)))
! ---> rho_f.xsf, rho_A1.xsf, rho_A2.xsf, rho_A3.xsf
! [n(r_vec),m_vec(r_vec)]
IMPLICIT NONE
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_input), INTENT(IN) :: input
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
LOGICAL, INTENT(IN) :: score, potnorm
CHARACTER(len=20), INTENT(IN) :: denName
TYPE(t_potden), INTENT(IN) :: denf
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA1
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA2
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA3
REAL :: tec, qint, phi0, angss
INTEGER :: i, j, ix, iy, iz, na, nplo, iv, iflag, nfile
INTEGER :: nplot, nt, jm, jspin, numInDen, numOutFiles
LOGICAL :: twodim, cartesian, xsf, unwind, polar
TYPE(t_potden), ALLOCATABLE :: den(:)
REAL, ALLOCATABLE :: xdnout(:)
REAL :: pt(3), vec1(3), vec2(3), vec3(3), &
zero(3), help(3), qssc(3), point(3)
INTEGER :: grid(3)
REAL :: rhocc(atoms%jmtd)
CHARACTER (len=20), ALLOCATABLE :: outFilenames(:)
CHARACTER (len=30) :: filename
CHARACTER (len=7) :: textline
REAL, PARAMETER :: eps = 1.0e-15
NAMELIST /plot/twodim,cartesian,unwind,vec1,vec2,vec3,grid,zero,phi0,filename
nfile = 120
IF (PRESENT(denA2)) THEN
ALLOCATE(den(4))
den(1) = denf
den(2) = denA1
den(3) = denA2
den(4) = denA3
numInDen = 4
numOutFiles = 4
ELSE IF (PRESENT(denA1)) THEN
ALLOCATE(den(2))
den(1) = denf
den(2) = denA1
numInDen = 2
numOutFiles = 2
ELSE
ALLOCATE(den(1))
den(1) = denf
numInDen = 1
numOutFiles = 1
END IF
DO i = 1, numInDen
IF ((numInDen.NE.4).AND.(score)) THEN
OPEN (17,file='cdnc',form='unformatted',status='old')
REWIND 17
DO jspin = 1, input%jspins
DO nt = 1, atoms%ntype
jm = atoms%jri(nt)
READ (17) (rhocc(j),j=1,jm)
DO j = 1, jm
den(i)%mt(j,0,nt,jspin) = den(i)%mt(j,0,nt,jspin) - rhocc(j)/2.0/SQRT(pi_const)
END DO
READ (17) tec
END DO
READ (17) qint
den(i)%pw(1,jspin) = den(i)%pw(1,jspin) - qint/cell%volint
END DO
CLOSE (17)
ELSE IF (score) THEN
CALL juDFT_error('Subtracting core charge in noco calculations not supported', calledby = 'plot')
END IF
END DO
IF (noco%l_ss) THEN
qssc = MATMUL(TRANSPOSE(cell%bmat),noco%qss)
END IF
! Open the plot_inp file for input
OPEN (18,file='plot_inp')
READ(18,'(i2,5x,l1,1x,a)') nplot,xsf,textline
polar = .FALSE.
IF ((noco%l_noco).AND.(numInDen.EQ.4)) THEN
polar = (textline(1:7)=='polar=T').OR.(textline(1:7)=='polar=t')
IF (polar) THEN
numOutFiles = 7
END IF
END IF
ALLOCATE(outFilenames(numOutFiles))
ALLOCATE(xdnout(numOutFiles))
IF (numOutFiles.EQ.1) THEN
outFilenames(1) = TRIM(denName)
ELSE IF (numOutFiles.EQ.2) THEN
outFilenames(1) = TRIM(denName) // '_f'
outFilenames(2) = TRIM(denName) // '_g'
ELSE
outFilenames(1) = TRIM(denName) // '_f'
outFilenames(2) = TRIM(denName) // '_A1'
outFilenames(3) = TRIM(denName) // '_A2'
outFilenames(4) = TRIM(denName) // '_A3'
IF (polar) THEN
outFilenames(5) = TRIM(denName) // '_Aabs'
outFilenames(6) = TRIM(denName) // '_Atheta'
outFilenames(7) = TRIM(denName) // '_Aphi'
END IF
END IF
! If xsf is specified we create input files for xcrysden
IF (.NOT.xsf) THEN
DO i = 1, numOutFiles
OPEN(nfile+i,file=TRIM(ADJUSTL(outFilenames(i)))//'.xsf',form='formatted')
CALL xsf_WRITE_atoms(nfile+i,atoms,input%film,oneD%odi%d1,cell%amat)
END DO
END IF
! Loop over all plots
DO nplo = 1, nplot
! the defaults
twodim = .TRUE.
cartesian = .TRUE.
grid = (/100,100,100/)
vec1 = (/0.,0.,0./)
vec2 = (/0.,0.,0./)
vec3 = (/0.,0.,0./)
zero = (/0.,0.,0./)
filename = "default"
READ(18,plot)
IF (twodim.AND.ANY(grid(1:2)<1)) &
CALL juDFT_error("Illegal grid size in plot",calledby="plot")
IF (.NOT.twodim.AND.ANY(grid<1)) &
CALL juDFT_error("Illegal grid size in plot",calledby="plot")
IF (twodim) grid(3) = 1
!calculate cartesian coordinates if needed
IF (.NOT.cartesian) THEN
vec1=matmul(cell%amat,vec1)
vec2=matmul(cell%amat,vec2)
vec3=matmul(cell%amat,vec3)
zero=matmul(cell%amat,zero)
END IF
!Open the file
IF (filename =="default") WRITE(filename,'(a,i2)') "plot",nplo
DO i = 1, numOutFiles
IF (xsf) THEN
CALL xsf_WRITE_header(nfile+i,twodim,filename,vec1,vec2,vec3,zero,grid)
ELSE
IF (numOutFiles.NE.1) THEN
OPEN (nfile+i,file = filename//outFilenames(i),form='formatted')
ELSE
OPEN (nfile+i,file = filename,form='formatted')
END IF
END IF
END DO
!loop over all points
DO iz = 0, grid(3)-1
DO iy = 0, grid(2)-1
DO ix = 0, grid(1)-1
point = zero + vec1*REAL(ix)/(grid(1)-1) +&
vec2*REAL(iy)/(grid(2)-1)
IF (.NOT.twodim) point = point + vec3*REAL(iz)/(grid(3)-1)
! Set region specific parameters for point
! Get MT sphere for point if point is in MT sphere
CALL getMTSphere(input,cell,atoms,oneD,point,nt,na,pt)
IF (na.NE.0) THEN
! In MT sphere
iv = 0
iflag = 1
ELSE IF (input%film.AND..NOT.oneD%odi%d1.AND.ABS(point(3))>=cell%z1) THEN
! In vacuum in 2D system
iv = 1
iflag = 0
pt(:) = point(:)
ELSE IF ((oneD%odi%d1).AND.(SQRT((point(1))**2 + (point(2))**2)>=cell%z1)) THEN
! In vacuum in 1D system
iv = 1
iflag = 0
pt(:) = point(:)
ELSE
! In interstitial region
iv = 0
iflag = 2
pt(:) = point(:)
END IF
DO i = 1, numInDen
CALL outcdn(pt,nt,na,iv,iflag,1,potnorm,stars,&
vacuum,sphhar,atoms,sym,cell,oneD,&
den(i),xdnout(i))
END DO
IF (na.NE.0) THEN
IF (noco%l_ss) THEN
! rotate magnetization "backward"
angss = DOT_PRODUCT(qssc,pt-atoms%pos(:,na))
help(1) = xdnout(2)
help(2) = xdnout(3)
xdnout(2) = +help(1)*COS(angss)+help(2)*SIN(angss)
xdnout(3) = -help(1)*SIN(angss)+help(2)*COS(angss)
! xdnout(2)=0. ; xdnout(3)=0. ; xdnout(4)=0.
END IF
END IF
IF (noco%l_ss .AND. (.NOT. unwind)) THEN
! rotate magnetization
angss = DOT_PRODUCT(qssc,point)
help(1) = xdnout(2)
help(2) = xdnout(3)
xdnout(2) = +help(1)*COS(angss) -help(2)*SIN(angss)
xdnout(3) = +help(1)*SIN(angss) +help(2)*COS(angss)
END IF
IF (polar) THEN
xdnout(5) = SQRT(ABS(xdnout(2)**2+xdnout(3)**2+xdnout(4)**2))
IF (xdnout(5)<eps) THEN
xdnout(5)= 0.0
!xdnout(6)= -tpi_const
!xdnout(7)= -tpi_const
xdnout(6)= 0.0
xdnout(7)= 0.0
ELSE
DO j = 1, 3
help(j) = xdnout(1+j)/xdnout(5)
END DO
IF (help(3)<0.5) THEN
xdnout(6)= ACOS(help(3))
ELSE
xdnout(6)= pi_const/2.0-ASIN(help(3))
END IF
IF (SQRT(ABS(help(1)**2+help(2)**2)) < eps) THEN
!xdnout(7)= -tpi_const
xdnout(7)= 0.0
ELSE
IF ( ABS(help(1)) > ABS(help(2)) ) THEN
xdnout(7)= ABS(ATAN(help(2)/help(1)))
ELSE
xdnout(7)= pi_const/2.0-ABS(ATAN(help(1)/help(2)))
END IF
IF (help(2)<0.0) THEN
xdnout(7)= -xdnout(7)
END IF
IF (help(1)<0.0) THEN
xdnout(7)= pi_const-xdnout(7)
END IF
!phi0=0
DO WHILE (xdnout(7) > +pi_const)
xdnout(7)= xdnout(7)-tpi_const
END DO
DO WHILE (xdnout(7) < -pi_const)
xdnout(7)= xdnout(7)+tpi_const
END DO
!IF (ABS(xdnout(2)-xdnout(3))<eps) THEN
! IF (xdnout(2)>0) THEN
! xdnout(7)=pi_const/4.0
! ELSE
! xdnout(7)=-3*pi_const/4.0
! END IF
!IF (xdnout(2)>eps) THEN
! xdnout(7)=ATAN(xdnout(3)/xdnout(2))
!ELSE IF (ABS(xdnout(2))<eps) THEN
! xdnout(7)=SIGN(1.0, xdnout(3))*pi_const/2.0
!ELSE IF ((xdnout(2)<-eps).AND.(xdnout(3)>=0.0)) THEN
! xdnout(7)=ATAN(xdnout(3)/xdnout(2))+pi_const
!ELSE
! xdnout(7)=ATAN(xdnout(3)/xdnout(2))-pi_const
!END IF
END IF
END IF
xdnout(6)= xdnout(6)/pi_const
xdnout(7)= xdnout(7)/pi_const
END IF ! (polar)
DO i = 1, numOutFiles
IF (xsf) THEN
WRITE(nfile+i,*) xdnout(i)
ELSE
WRITE(nfile+i,'(4e15.7)') point ,xdnout(i)
END IF
END DO
END DO
END DO
END DO !z-loop
DO i = 1, numOutFiles
IF (xsf) THEN
CALL xsf_WRITE_endblock(nfile+i,twodim)
ELSE
CLOSE(nfile+i)
END IF
END DO
END DO !nplot
CLOSE(18)
IF (xsf) THEN
DO i = 1, numOutFiles
CLOSE(nfile+i)
END DO
END IF
DEALLOCATE(xdnout, outFilenames)
END SUBROUTINE savcsv
SUBROUTINE vectorplot(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
......@@ -1035,7 +1380,7 @@ CONTAINS
factor = 2.0
denName = 'vTot'
score = .FALSE.
potnorm = .FALSE.
potnorm = .TRUE.
IF (input%jspins.EQ.2) THEN
IF (noco%l_noco) THEN
CALL matrixplot(stars, atoms, sphhar, vacuum, input, oneD, sym, &
......@@ -1055,7 +1400,7 @@ CONTAINS
END SUBROUTINE procplot
SUBROUTINE plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, div, phiPot, divGrx, divGry, divGrz)!, &
noco, div, phiPot)!, divGrx, divGry, divGrz, &
!xcBmodx, xcBmody, xcBmodz)
IMPLICIT NONE
......@@ -1069,7 +1414,7 @@ CONTAINS
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), INTENT(IN) :: div, phiPot, divGrx, divGry, divGrz!, &
TYPE(t_potden), INTENT(IN) :: div, phiPot!, divGrx, divGry, divGrz, &
!xcBmodx, xcBmody, xcBmodz
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
......@@ -1078,14 +1423,14 @@ CONTAINS
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'modPot ', phiPot)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'divPotx ', divGrx)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
! .FALSE., .TRUE., 'divPotx ', divGrx)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'divPoty ', divGry)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
! .FALSE., .TRUE., 'divPoty ', divGry)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'divPotz ', divGrz)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
! .FALSE., .TRUE., 'divPotz ', divGrz)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
! .FALSE., .TRUE., 'xcBmodx ', xcBmodx)
......
......@@ -6,6 +6,7 @@ vgen/modcyli.f
vgen/modcylk.f
vgen/od_mkgxyz3.f
vgen/od_mkgz.f
vgen/sphcoord.f
vgen/vacp5_0.f
vgen/vacp5_z.f
vgen/visp5_0.f
......@@ -49,7 +50,6 @@ vgen/xy_av_den.f90
vgen/VYukawaFilm.f90
vgen/divergence.f90
vgen/xcBfield.f90
vgen/sphcoord.f90
)
#vdW Stuff
set(fleur_F90 ${fleur_F90}
......
......@@ -36,10 +36,12 @@ CONTAINS
REAL, ALLOCATABLE :: thet(:), phi(:), div_temp(:, :)
REAL, ALLOCATABLE :: xcBx_mt(:,:,:), xcBy_mt(:,:,:), xcBz_mt(:,:,:)
REAL :: r, th, ph
INTEGER :: jr, k, nsp, kt, i
REAL :: r, th, ph, eps
INTEGER :: jr, k, nsp, kt, i, lh, lhmax
nsp = atoms%nsp()
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
eps=1.e-10
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,1),grady%gr(3,atoms%jri(n)*nsp,1), &
gradz%gr(3,atoms%jri(n)*nsp,1))
......@@ -60,12 +62,6 @@ CONTAINS
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
!DO jr=1,atoms%jri(n)
! xcBx_mt(jr,0:,:) = xcB(1)%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
! xcBy_mt(jr,0:,:) = xcB(2)%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
! xcBz_mt(jr,0:,:) = xcB(3)%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
!ENDDO
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(1)%mt(:,0:,n,:), n, gradx)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(2)%mt(:,0:,n,:), n, grady)
CALL mt_to_grid(.TRUE., 1, atoms, sphhar, bxc(3)%mt(:,0:,n,:), n, gradz)
......@@ -76,14 +72,33 @@ CONTAINS
DO k = 1, nsp
th = thet(k)
ph = phi(k)
div_temp(kt+k,1) = (SIN(th)*COS(ph)*gradx%gr(1,kt+k,1) + SIN(th)*SIN(ph)*grady%gr(1,kt+k,1) + COS(th)*gradz%gr(1,kt+k,1))*r**2&
+(COS(th)*COS(ph)*gradx%gr(2,kt+k,1) + COS(th)*SIN(ph)*grady%gr(2,kt+k,1) - SIN(th)*gradz%gr(2,kt+k,1))*r&
-(SIN(ph)*gradx%gr(3,kt+k,1) - COS(ph)*grady%gr(3,kt+k,1))*r/SIN(th)
div_temp(kt+k,1) = (SIN(th)*COS(ph)*gradx%gr(1,kt+k,1) + SIN(th)*SIN(ph)*grady%gr(1,kt+k,1) + COS(th)*gradz%gr(1,kt+k,1))&
+(COS(th)*COS(ph)*gradx%gr(2,kt+k,1) + COS(th)*SIN(ph)*grady%gr(2,kt+k,1) - SIN(th)*gradz%gr(2,kt+k,1))/r&
-(SIN(ph)*gradx%gr(3,kt+k,1) - COS(ph)*grady%gr(3,kt+k,1))/(r*SIN(th))
ENDDO ! k
kt = kt+nsp
ENDDO ! jr
CALL mt_from_grid(atoms, sphhar, n, 1, div_temp, div%mt(:,0:,n,:))
DO jr = 1, atoms%jri(n)
DO lh=0, lhmax
IF (ABS(div%mt(jr,lh,n,1))<eps) THEN
div%mt(jr,lh,n,:)=0.0
END IF
END DO
END DO
kt = 0
DO jr = 1, atoms%jri(n)
r =atoms%rmsh(jr, n)
DO k = 1, nsp
div_temp(kt+k,1) = div_temp(kt+k,1)*r**2
ENDDO ! k
kt = kt+nsp
ENDDO ! jr
CALL finish_mt_grid
......@@ -112,7 +127,7 @@ CONTAINS
TYPE(t_potden), INTENT(INOUT) :: div
TYPE(t_gradients) :: gradx,grady,gradz
REAL, ALLOCATABLE :: div_temp(:, :)
INTEGER :: i,ifftxc3
......
......@@ -9,7 +9,7 @@ module m_mpmom
! Coulomb potential. This is used in the preconditioning of
! the SCF iteration for metallic systems.
!
! qlmo(m,l,n) : mult.mom. of the mufftn-tin charge density
! qlmo(m,l,n) : mult.mom. of the muffin-tin charge density
! qlmp(m,l,n) : mult.mom. of the plane-wave charge density
! qlm (m,l,n) : (output) difference of the former quantities
!
......
!--------------------------------------------------------------------------------
! Copyright (c) 2019 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_sphcoord
!--------------------------------------------------------------------------------
! Calculates the polar angles theta and phi of a vector with components
! x, y and z.
!
! Based on pol_angle by P. Kurz.
!--------------------------------------------------------------------------------
CONTAINS
SUBROUTINE sphcoord(x,y,z,theta,phi)
USE m_constants
IMPLICIT NONE
REAL, INTENT(IN) :: x, y, z
REAL, INTENT(OUT) :: theta, phi
REAL, PARAMETER :: eps = 1.0e-15
REAL r, rho
rho = SQRT(x**2+y**2)
r = SQRT(x**2+y**2+z**2)
IF (r.LT.eps) THEN
theta = 0.0
phi = 0.0
ELSE
theta = ACOS(z/r)
IF (rho.LT.eps) THEN
phi = 0.0
ELSE IF ((x/rho).LT.eps) THEN
phi = SIGN(1.0, y)*pi_const/2.0
ELSE IF (x.GT.0.0) THEN
phi = ATAN(y/x)!pi/2 - SIGN(1.0, vx)*pi/2
ELSE IF (y.GE.0.0) THEN
phi = ATAN(y/x)+pi_const
ELSE
phi = ATAN(y/x)-pi_const
ENDIF
ENDIF
END SUBROUTINE sphcoord
END MODULE m_sphcoord
!--------------------------------------------------------------------------------