subroutine aptcsyv (n1, n2, iunit, u, v, w, au, av, aw, np,
& tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCSYV
c
c call aptcsyv (n1, n2, iunit, u, v, w, au, av, aw, np,
c & tol, nerr)
c
c Version: aptcsyv Updated 1989 November 27 14:00.
c aptcsyv 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) and
c the associated bound vectors a = (au, av, aw) from
c coordinate system n1 to coordinate system n2. Allowed
c coordinate systems are cartesian, cylindrical, spherical.
c Note: if (u, v, w) is at the origin, then au and av are
c independent of the coordinate system, and aw changes sign
c between the spherical coordinate system and the other two.
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, au, av, aw, np, tol.
c
c Output: u, v, w, au, av, aw, nerr.
c
c Glossary:
c
c au,av,aw In/Out The u, v, w components of a bound vector at (u, v, w).
c Size np.
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), and
c number of bound vectors (au, av, aw).
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 Any spatial coordinate or vector component with
c a final value less than the estimated error in its
c calculation, based on tol, will be truncated to
c zero.
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---- Component u of bound vector.
dimension au (1)
c---- Component v of bound vector.
dimension av (1)
c---- Component w of bound vector.
dimension aw (1)
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 /laptcsyv/ angfin
c---- Output angle conversion factor.
common /laptcsyv/ angfout
c---- Estimated truncation error in au.
common /laptcsyv/ auerr
c---- Initial value of au.
common /laptcsyv/ aus
c---- Estimated truncation error in av.
common /laptcsyv/ averr
c---- Initial value of av.
common /laptcsyv/ avs
c---- Estimated truncation error in aw.
common /laptcsyv/ awerr
c---- Initial value of aw.
common /laptcsyv/ aws
c---- Cosine (v) (cyl or spherical).
common /laptcsyv/ costh
c---- Cosine (w) (spherical).
common /laptcsyv/ cosph
c---- A very small number.
common /laptcsyv/ fuz
c---- Index in u, v, w, au, av, aw.
common /laptcsyv/ n
c---- Mathematical constant, pi.
common /laptcsyv/ pi
c---- Distance from z axis.
common /laptcsyv/ rcyl
c---- Sine (v) (cyl or spherical).
common /laptcsyv/ sinth
c---- Sine (w) (spherical).
common /laptcsyv/ sinph
c---- Initial value of u.
common /laptcsyv/ us
c---- Initial value of v.
common /laptcsyv/ vs
c---- Initial value of w.
common /laptcsyv/ ws
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcsyv 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',12x,'v',12x,'w',12x,'au',11x,'av',11x,'aw' /
cbug & (i2,1p6e13.5))
cbug write ( 3, 9902) (n, u(n), v(n), w(n), au(n), av(n), aw(n),
cbug & 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
aus = au(n)
avs = av(n)
au(n) = aus * costh + avs * sinth
auerr = tol * (abs (aus * costh) + abs (avs * sinth))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
av(n) = -aus * sinth + avs * costh
averr = tol * (abs (aus * sinth) + abs (avs * costh))
if (abs (av(n)) .lt. averr) then
av(n) = 0.0
endif
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
aus = au(n)
avs = av(n)
aws = aw(n)
au(n) = aus * costh * sinph + avs * sinth * sinph +
& aws * cosph
auerr = tol * (abs (aus * costh * sinph) +
& abs (avs * sinth * sinph) +
& abs (aws * cosph))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
av(n) = -aus * sinth + avs * costh
averr = tol * (abs (aus * sinth) + abs (avs * costh))
if (abs (av(n)) .lt. averr) then
av(n) = 0.0
endif
aw(n) = aus * costh * cosph + avs * sinth * cosph -
& aws * sinph
awerr = tol * (abs (aus * costh * cosph) +
& abs (avs * sinth * cosph) +
& abs (aws * sinph))
if (abs (aw(n)) .lt. awerr) then
aw(n) = 0.0
endif
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
aus = au(n)
avs = av(n)
au(n) = aus * costh - avs * sinth
auerr = tol * (abs (aus * costh) + abs (avs * sinth))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
av(n) = aus * sinth + avs * costh
averr = tol * (abs (aus * sinth) + abs (avs * costh))
if (abs (av(n)) .lt. averr) then
av(n) = 0.0
endif
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
aus = au(n)
avs = av(n)
aws = aw(n)
au(n) = aus * sinph + aws * cosph
auerr = tol * (abs (aus * sinph) + abs (aws * cosph))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
aw(n) = aus * cosph - aws * sinph
awerr = tol * (abs (aus * cosph) + abs (aws * sinph))
if (abs (aw(n)) .lt. awerr) then
aw(n) = 0.0
endif
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
aus = au(n)
avs = av(n)
aws = aw(n)
au(n) = aus * costh * sinph - avs * sinth +
& aws * costh * cosph
auerr = tol * (abs (aus * costh * sinph) +
& abs (avs * sinth) +
& abs (aws * costh * cosph))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
av(n) = aus * sinth * sinph + avs * costh +
& aws * sinth * cosph
averr = tol * (abs (aus * sinth * sinph) +
& abs (avs * costh) +
& abs (aws * sinth * cosph))
if (abs (av(n)) .lt. averr) then
av(n) = 0.0
endif
aw(n) = aus * cosph - aws * sinph
awerr = tol * (abs (aus * cosph) + abs (aws * sinph))
if (abs (aw(n)) .lt. awerr) then
aw(n) = 0.0
endif
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
aus = au(n)
aws = aw(n)
au(n) = aus * sinph + aws * cosph
auerr = tol * (abs (aus * sinph) + abs (aws * cosph))
if (abs (au(n)) .lt. auerr) then
au(n) = 0.0
endif
aw(n) = aus * cosph - aws * sinph
awerr = tol * (abs (aus * cosph) + abs (aws * sinph))
if (abs (aw(n)) .lt. awerr) then
aw(n) = 0.0
endif
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',12x,'v',12x,'w',12x,'au',11x,'av',11x,'aw' /
cbug & (i2,1p6e13.5))
cbug write ( 3, 9903) (n, u(n), v(n), w(n), au(n), av(n), aw(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptcsyv. (+1 line.)
end
UCRL-WEB-209832