subroutine aptcsys (n1, n2, iunit, u, v, w, np, tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCSYS
c
c call aptcsys (n1, n2, iunit, u, v, w, np, tol, nerr)
c
c Version: aptcsys Updated 1990 November 27 14:00.
c aptcsys Originated 1989 November 2 14:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To transform the np points (u, v, w) from
c coordinate system n1 to coordinate system n2. Allowed
c coordinate systems are cartesian (0), cylindrical (1), or
c spherical (2).
c Angles may be in degrees (iunit = 0) or radians (iunit = 1).
c Results will be truncated to zero if less than the estimated
c error in their calculation, based on tol. Disallowed input
c values of n1, n2, iunit, or np are indicated by a nonzero
c value of nerr.
c
c Input: n1, n2, iunit, u, v, w, np, tol.
c
c Output: u, v, w, nerr.
c
c Glossary:
c
c iunit Input Indicates unit to be used for angles:
c 0 if angles are in degrees.
c 1 if angles are in radians.
c
c n1 Input Indicates initial coordinate system type:
c 0 for cartesian coordinates. u = x, v = y, w = z.
c 1 for cylindrical coordinates. u = radius from z
c axis, v = angle in xy plane, counterclockwise from
c x axis, w = z.
c 2 for spherical coordinates. u = radius from origin,
c v = angle in xy plane, counterclockwise from
c x axis, w = angle from z axis.
c
c n2 Input Indicates final coordinate system. See n1.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if either n1 or n2 is not 0, 1, or 2.
c 3 if n1 = n2.
c 4 if iunit is not 0 or 1.
c
c np Input Number of points (u, v, w).
c
c tol Input Numerical tolerance limit. Any angle with a sine
c or cosine (absolute value) less than tol will be
c adjusted to make its sine or cosine = 0.
c
c u,v,w In/Out The coordinates of a point. Size np.
c See n1, n2. Will be truncated to zero, if less than
c the estimated error in their calculation,
c base on tol.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- First spatial coordinate.
dimension u (1)
c---- Second spatial coordinate.
dimension v (1)
c---- Third spatial coordinate.
dimension w (1)
c.... Local variables.
c---- Input angle conversion factor.
common /laptcsys/ angfin
c---- Output angle conversion factor.
common /laptcsys/ angfout
c---- Cosine (v) (cyl or spherical).
common /laptcsys/ costh
c---- Cosine (w) (spherical).
common /laptcsys/ cosph
c---- A very small number.
common /laptcsys/ fuz
c---- Index in u, v, w.
common /laptcsys/ n
c---- Mathematical constant, pi.
common /laptcsys/ pi
c---- Distance from z axis.
common /laptcsys/ rcyl
c---- Sine (v) (cyl or spherical).
common /laptcsys/ sinth
c---- Sine (w) (spherical).
common /laptcsys/ sinph
c---- Initial value of u.
common /laptcsys/ us
c---- Initial value of v.
common /laptcsys/ vs
c---- Initial value of w.
common /laptcsys/ ws
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcsys converting units from n1=',i2,' to n2=',i2,
cbug & '. iunit=',i2,'.' /
cbug & ' n1,n2: 0 = x,y,z, 1 = r,theta,z, 2 = rho,theta,phi.' /
cbug & ' iunit = 0 (degrees), 1 (radians).')
cbug write ( 3, 9901) n1, n2, iunit
cbugc***DEBUG ends.
c.... Initialize.
c---- Mathematical constant, pi.
pi = 3.14159265358979323
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((n1 .lt. 0) .or. (n1 .gt. 2)) then
nerr = 2
go to 210
endif
if ((n2 .lt. 0) .or. (n2 .gt. 2)) then
nerr = 2
go to 210
endif
if (n1 .eq. n2) then
nerr = 3
go to 210
endif
if ((iunit .ne. 0) .and. (iunit .ne. 1)) then
nerr = 4
go to 210
endif
cbugc***DEBUG begins.
cbug 9902 format (/ ' initial values.' /
cbug & ' n u',21x,'v',21x,'w' / (i5,1p3e22.14))
cbug write ( 3, 9902) (n, u(n), v(n), w(n), n = 1, np)
cbugc***DEBUG ends.
c.... Find factors to convert angle units.
c---- Angles are in degrees.
if (iunit .eq. 0) then
angfin = pi / 180.0
angfout = 180.0 / pi
c---- Angles are in radians.
else
angfin = 1.0
angfout = 1.0
c---- Tested iunit.
endif
c.... Find type of conversion to do, and do it.
c---- Initial system is cartesian.
if (n1 .eq. 0) then
c---- Cartesian to cylindrical.
if (n2 .eq. 1) then
c---- Loop over points.
do 110 n = 1, np
us = u(n)
vs = v(n)
u(n) = sqrt (us**2 + vs**2)
costh = (us + fuz) / (u(n) + fuz)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = vs / (u(n) + fuz)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
v(n) = atan2 (sinth, costh) * angfout
c---- End of loop over points.
110 continue
c---- Cartesian to spherical.
elseif (n2 .eq. 2) then
c---- Loop over points.
do 120 n = 1, np
us = u(n)
vs = v(n)
ws = w(n)
rcyl = sqrt (us**2 + vs**2)
costh = (us + fuz) / (rcyl + fuz)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = vs / (rcyl + fuz)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
u(n) = sqrt (us**2 + vs**2 + ws**2)
cosph = ws / (u(n) + fuz)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = (rcyl + fuz) / (u(n) + fuz)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
v(n) = atan2 (sinth, costh) * angfout
w(n) = atan2 (sinph, cosph) * angfout
c---- End of loop over points.
120 continue
c---- Tested n2 for this n1.
endif
c---- Initial system is cylindrical.
elseif (n1 .eq. 1) then
c---- Cylindrical to cartesian.
if (n2 .eq. 0) then
c---- Loop over points.
do 130 n = 1, np
us = u(n)
vs = v(n) * angfin
costh = cos (vs)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = sin (vs)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
u(n) = us * costh
v(n) = us * sinth
c---- End of loop over points.
130 continue
c---- Cylindrical to spherical.
elseif (n2 .eq. 2) then
c---- Loop over points.
do 140 n = 1, np
us = u(n)
ws = w(n)
u(n) = sqrt (us**2 + ws**2)
cosph = ws / (u(n) + fuz)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = (us + fuz) / (u(n) + fuz)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
w(n) = atan2 (sinph, cosph) * angfout
c---- End of loop over points.
140 continue
c---- Tested n2 for this n1.
endif
c---- Initial system is spherical.
elseif (n1 .eq. 2) then
c---- Spherical to cartesian.
if (n2 .eq. 0) then
c---- Loop over points.
do 150 n = 1, np
us = u(n)
vs = v(n) * angfin
ws = w(n) * angfin
costh = cos (vs)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = sin (vs)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
cosph = cos (ws)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = sin (ws)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
u(n) = us * costh * sinph
v(n) = us * sinth * sinph
w(n) = us * cosph
c---- End of loop over points.
150 continue
c---- Spherical to cylindrical.
elseif (n2 .eq. 1) then
c---- Loop over points.
do 160 n = 1, np
us = u(n)
ws = w(n) * angfin
cosph = cos (ws)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = sin (ws)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
u(n) = us * sinph
w(n) = us * cosph
c---- End of loop over points.
160 continue
c---- Tested n2 for this n1.
endif
c---- Tested n1.
endif
cbugc***DEBUG begins.
cbug 9903 format (/ ' final values.' /
cbug & ' n u',21x,'v',21x,'w' / (i5,1p3e22.14))
cbug write ( 3, 9903) (n, u(n), v(n), w(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptcsys. (+1 line.)
end
UCRL-WEB-209832