subroutine aptqvac (au, av, bu, bv, cu, cv, du, dv, np, tol, nopt,
& arpa, arpb, arpc, arpd, ntype, qu, qv, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQVAC
c
c call aptqvac (au, av, bu, bv, cu, cv, du, dv, np, tol, nopt,
c & arpa, arpb, arpc, arpd, ntype, qu, qv, nerr)
c
c Version: aptqvac Updated 1990 November 29 15:30.
c aptqvac Originated 1990 February 8 14:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the shape and size of np quadrangles in the uv plane,
c with vertices a = (au, av), b = (bu, bv), c = (cu, cv) and
c d = (du, dv) counterclockwise around the quadrangle.
c Optionally (nopt = 0 or 2), the vertex parallelogram areas arpa,
c arpb, arpc and arpd will be returned, and (nopt = 1 or 2) the
c number and location of any concave (negative area) vertices
c will be returned, encoded in ntype, and if the quadrangle is a
c bowtie, the intersection of the sides q = (qu, qv), will be
c returned. The net area of the quadrangle is
c 0.25 * (arpa + arpb + arpc + arpd).
c Flag nerr indicates any input error.
c
c Input: au, av, bu, bv, cu, cv, du, dv, np, tol, nopt.
c
c Output: arpa, arpb, arpc, arpd, ntype, qu, qv, nerr.
c
c Calls: aptvdic, aptvaxc
c
c
c Glossary:
c
c arpa Output Parallelogram area at vertex "a" (nopt = 0 or 2).
c Size np.
c
c arpb Output Parallelogram area at vertex "b" (nopt = 0 or 2).
c Size np.
c
c arpc Output Parallelogram area at vertex "c" (nopt = 0 or 2).
c Size np.
c
c arpd Output Parallelogram area at vertex "d" (nopt = 0 or 2).
c Size np.
c
c au, av Input The u and v coordinates of vertex "a" of quadrangle.
c Size np.
c
c bu, bv Input The u and v coordinates of vertex "b" of quadrangle.
c Size np.
c
c cu, cv Input The u and v coordinates of vertex "c" of quadrangle.
c Size np.
c
c du, dv Input The u and v coordinates of vertex "d" of quadrangle.
c Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if nopt is not 0, 1 or 2.
c
c nopt Input Output option:
c 0 or 2 to return vertex parallelogram areas.
c 1 or 2 to return shape type code ntype.
c 2 to return areas, ntype, and bowtie intersections.
c
c np Input Size of arrays au, av, bu, bv, cu, cv, du, dv,
c and if nopt = 0 or 2, arpa, arpb, arpc and arpd,
c and if nopt = 1 or 2, ntype.
c
c ntype Output Shape type of quadrangle (nopt = 1 or 2). Size np.
c 0 if all vertices are convex. Regular.
c 11 if only vertex "a" is concave. Boomerang.
c 12 if only vertex "b" is concave. Boomerang.
c 13 if only vertex "c" is concave. Boomerang.
c 14 if only vertex "d" is concave. Boomerang.
c 21 if only vertices "c" and "d" are concave. Bowtie.
c 22 if only vertices "d" and "a" are concave. Bowtie.
c 23 if only vertices "a" and "b" are concave. Bowtie.
c 24 if only vertices "b" and "c" are concave. Bowtie.
c 31 if only vertex "a" is convex. Inverted boomerang.
c 32 if only vertex "b" is convex. Inverted boomerang.
c 33 if only vertex "c" is convex. Inverted boomerang.
c 34 if only vertex "d" is convex. Inverted boomerang.
c 40 if all vertices are concave. Inverted regular.
c
c qu, qv Output The u and v coordinates of the intersection "q" of
c two opposite sides of a bowtied quadrangle.
c Returned only if nopt = 1 or 2, and ntype = 21
c to 24. Size np, only if nopt = 1 or 2.
c
c tol Input Numerical tolerance limit. With 64-bit floating
c point numbers, 1.e-5 to 1.e-11 is recommended.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Parallelogram area at vertex "a".
dimension arpa (1)
c---- Parallelogram area at vertex "b".
dimension arpb (1)
c---- Parallelogram area at vertex "c".
dimension arpc (1)
c---- Parallelogram area at vertex "d".
dimension arpd (1)
c---- Coordinate u of vertex "a".
dimension au (1)
c---- Coordinate v of vertex "a".
dimension av (1)
c---- Coordinate u of vertex "b".
dimension bu (1)
c---- Coordinate v of vertex "b".
dimension bv (1)
c---- Coordinate u of vertex "c".
dimension cu (1)
c---- Coordinate v of vertex "c".
dimension cv (1)
c---- Coordinate u of vertex "d".
dimension du (1)
c---- Coordinate v of vertex "d".
dimension dv (1)
c---- Shape type number of quadrangle.
dimension ntype (1)
c---- Coordinate u of bowtie point "q".
dimension qu (1)
c---- Coordinate v of bowtie point "q".
dimension qv (1)
c.... Local variables.
c---- Magnitude of edge vector "ab".
common /laptqvac/ ablen (64)
c---- Component u of edge vector "ab".
common /laptqvac/ abu (64)
c---- Component v of edge vector "ab".
common /laptqvac/ abv (64)
c---- Parallelogram area at vertex "a".
common /laptqvac/ arla (64)
c---- Parallelogram area at vertex "b".
common /laptqvac/ arlb (64)
c---- Parallelogram area at vertex "c".
common /laptqvac/ arlc (64)
c---- Parallelogram area at vertex "d".
common /laptqvac/ arld (64)
c---- Magnitude of edge vector "bc".
common /laptqvac/ bclen (64)
c---- Component u of edge vector "bc".
common /laptqvac/ bcu (64)
c---- Component v of edge vector "bc".
common /laptqvac/ bcv (64)
c---- A very big number.
common /laptqvac/ big
c---- Magnitude of edge vector "cd".
common /laptqvac/ cdlen (64)
c---- Component u of edge vector "cd".
common /laptqvac/ cdu (64)
c---- Component v of edge vector "cd".
common /laptqvac/ cdv (64)
c---- Magnitude of edge vector "da".
common /laptqvac/ dalen (64)
c---- Component u of edge vector "da".
common /laptqvac/ dau (64)
c---- Component v of edge vector "da".
common /laptqvac/ dav (64)
c---- Fractional distance along "ab" of "q".
common /laptqvac/ fab
c---- Fractional distance along "bc" of "q".
common /laptqvac/ fbc
c---- Fractional distance along "cd" of "q".
common /laptqvac/ fcd
c---- Fractional distance along "da" of "q".
common /laptqvac/ fda
c---- A very small number.
common /laptqvac/ fuz
c---- Index in shpal array.
common /laptqvac/ n
c---- First index of subset of data.
common /laptqvac/ n1
c---- Last index of subset of data.
common /laptqvac/ n2
c---- Index in external array.
common /laptqvac/ nn
c---- 1 if ntype values not all zero.
common /laptqvac/ nodd
c---- Size of current subset of data.
common /laptqvac/ ns
c---- Component u of temporary bowtie point.
common /laptqvac/ qabu
c---- Component v of temporary bowtie point.
common /laptqvac/ qabv
c---- Component u of temporary bowtie point.
common /laptqvac/ qbcu
c---- Component v of temporary bowtie point.
common /laptqvac/ qbcv
c---- Component u of temporary bowtie point.
common /laptqvac/ qcdu
c---- Component v of temporary bowtie point.
common /laptqvac/ qcdv
c---- Component u of temporary bowtie point.
common /laptqvac/ qdau
c---- Component v of temporary bowtie point.
common /laptqvac/ qdav
cbugc***DEBUG begins.
cbugc---- 0 if all quadrangles are convex.
cbug common /laptqvac/ noddt
cbug 9901 format (/ 'aptqvac finding quadrangle shapes. nopt=',i3 /
cbug & (i3,' au=',1pe22.14,' av=',1pe22.14 /
cbug & ' bu=',1pe22.14,' bv=',1pe22.14 /
cbug & ' cu=',1pe22.14,' cv=',1pe22.14 /
cbug & ' du=',1pe22.14,' dv=',1pe22.14))
cbug write (3, 9901) nopt, (n, au(n), av(n), bu(n), bv(n),
cbug & cu(n), cv(n), du(n), dv(n), n = 1, np)
cbug noddt = 0
cbugc***DEBUG ends.
c.... initialize.
c---- A very big number.
big = 1.e+99
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 ((nopt .lt. 0) .and. (nopt .gt. 2)) then
nerr = 2
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the edge vectors.
call aptvdic (au(n1), av(n1), bu(n1), bv(n1), ns, tol,
& abu, abv, ablen, nerr)
call aptvdic (bu(n1), bv(n1), cu(n1), cv(n1), ns, tol,
& bcu, bcv, bclen, nerr)
call aptvdic (cu(n1), cv(n1), du(n1), dv(n1), ns, tol,
& cdu, cdv, cdlen, nerr)
call aptvdic (du(n1), dv(n1), au(n1), av(n1), ns, tol,
& dau, dav, dalen, nerr)
c.... Find the parallelogram areas at the vertices. Positive if convex.
call aptvaxc (dau, dav, abu, abv, ns, tol, arla, nerr)
call aptvaxc (abu, abv, bcu, bcv, ns, tol, arlb, nerr)
call aptvaxc (bcu, bcv, cdu, cdv, ns, tol, arlc, nerr)
call aptvaxc (cdu, cdv, dau, dav, ns, tol, arld, nerr)
c.... See if the shape type code is to be returned.
c++++ Find ntype.
if ((nopt .eq. 1) .or. (nopt .eq. 2)) then
c.... Find and identify all positive convex quadrangles.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 0
else
ntype(nn) = 1
endif
c---- End of loop over subset of data.
120 continue
nodd = 0
c---- Loop over subset of data.
do 130 nn = n1, n2
nodd = nodd + ntype(nn)
c---- End of loop over subset of data.
130 continue
cbugc***DEBUG begins.
cbug noddt = noddt + nodd
cbugc***DEBUG ends.
c.... Find and identify all other quadrangles.
c---- Odd shapes exist.
if (nodd .gt. 0) then
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
c.... Test for a convex quadrangle.
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 0
else
ntype(nn) = -1
endif
c.... Test for a boomerang (one concave vertex).
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 11
endif
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 12
endif
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 13
endif
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 14
endif
c.... Test for a bowtie (two concave vertices).
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 21
endif
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 22
endif
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 23
endif
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 24
endif
c.... Test for a negative boomerang (three concave vertices).
if ((arla(n) .ge. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 31
endif
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .ge. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 32
endif
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .ge. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 33
endif
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .ge. 0.0) ) then
ntype(nn) = 34
endif
c.... Test for an inside-out convex quadrangle (four concave vertices).
if ((arla(n) .lt. 0.0) .and.
& (arlb(n) .lt. 0.0) .and.
& (arlc(n) .lt. 0.0) .and.
& (arld(n) .lt. 0.0) ) then
ntype(nn) = 40
endif
c---- End of loop over subset of data.
140 continue
c---- Tested nodd.
endif
c---- Tested nopt.
endif
c.... See if the vertex parallelogram areas are to be returned.
c++++ Find parallelogram areas.
if ((nopt .eq. 0) .or. (nopt .eq. 2)) then
c---- Loop over subset of data.
do 150 n = 1, ns
nn = n + n1 - 1
arpa(nn) = arla(n)
arpb(nn) = arlb(n)
arpc(nn) = arlc(n)
arpd(nn) = arld(n)
c---- End of loop over subset of data.
150 continue
c---- Tested nopt.
endif
c.... See if the edge intersections of bowtied quadrangles are to be
c.... returned.
c++++ Find bowtie point.
if ((nopt .eq. 2) .and. (nodd .gt. 0)) then
c---- Loop over subset of data.
do 160 n = 1, ns
nn = n + n1 - 1
fab = arld(n) / (arla(n) - arlb(n) + fuz)
fbc = arla(n) / (arlb(n) - arlc(n) + fuz)
fcd = arlb(n) / (arlc(n) - arld(n) + fuz)
fda = arlc(n) / (arld(n) - arla(n) + fuz)
qabu = au(nn) + fab * abu(n)
qabv = av(nn) + fab * abv(n)
qbcu = bu(nn) + fbc * bcu(n)
qbcv = bv(nn) + fbc * bcv(n)
qcdu = cu(nn) + fcd * cdu(n)
qcdv = cv(nn) + fcd * cdv(n)
qdau = du(nn) + fda * dau(n)
qdav = dv(nn) + fda * dav(n)
if (ntype(nn) .eq. 21) then
qu(nn) = qbcu
qv(nn) = qbcv
else
qu(nn) = -big
qv(nn) = -big
endif
if (ntype(nn) .eq. 22) then
qu(nn) = qcdu
qv(nn) = qcdv
endif
if (ntype(nn) .eq. 23) then
qu(nn) = qdau
qv(nn) = qdav
endif
if (ntype(nn) .eq. 24) then
qu(nn) = qabu
qv(nn) = qabv
endif
c---- End of loop over subset of data.
160 continue
c---- Tested nopt and nodd.
endif
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptqvac results. noddt=',i3)
cbug 9903 format (i3,' arpa=',1pe22.14,' arpb=',1pe22.14 /
cbug & ' arpc=',1pe22.14,' arpd=',1pe22.14)
cbug 9904 format (i3,' ntype=',i3)
cbug 9905 format (i3,' arpa=',1pe22.14,' arpb=',1pe22.14 /
cbug & ' arpc=',1pe22.14,' arpd=',1pe22.14 /
cbug & ' qu= ',1pe22.14,' qv= ',1pe22.14,' ntype=',i3)
cbug write ( 3, 9902) noddt
cbug if (nopt .eq. 0) then
cbug write ( 3, 9903) (n, arpa(n), arpb(n),
cbug & arpc(n), arpd(n), n = 1, np)
cbug endif
cbug if (nopt .eq. 1) then
cbug write ( 3, 9904) (n, ntype(n), n = 1, np)
cbug endif
cbug if (nopt .eq. 2) then
cbug write ( 3, 9905) (n, arpa(n), arpb(n),
cbug & arpc(n), arpd(n), qu(n), qv(n), ntype(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptqvac. (+1 line.)
end
UCRL-WEB-209832