subroutine polysizing (coefexp, np, list, x, y, z, nx, ny, nz, ier)
real coefexp
integer np, ier
integer :: list(1:np)
real :: x(1:np), y(1:np), z(1:np), nx(1:np), ny(1:np), nz(1:np)

!***********************************************************
!   This subroutine take and input polygon and expands it or shrinkes it,
! returning a new polygon with the same number of points. Wher ethe polygon
! collapses on itself, multiple points are mapped to the same coordinate.
! the point order does not change.

! On input:
!       COEFEXP = distance the polygon is to shrink in fraction of unit sphere radius
!       NP = number of points in polygon array
!       LIST = array of index values for order of polygon
!       X,Y,Z = Cartesian coordinates of points of input polygon

! Input parameter COEFEXP is altered on return
! All other input parameters are not altered by this function.

! On output:
!       NX,NY,NZ = Cartesian coordinates of points of resulting polygon
!       COEFEXP = distance the polygon has yet to shrink in fraction of unit sphere radius

!       IER = Error indicator:
!             IER = -1 returned before full distance reached. Locally, full distance was reached
!             IER = 0 if no errors were encountered.
!             IER = 1 polygon with only one valid segment
!             IER = 2 polygon with no valid segments
!***********************************************************

logical :: crossing
integer :: idxb, idx, idxa, jdx, dir_sign, idx1, idx2, jdx1, kdx
integer :: n_cross, cross_list(1:4*np), n_pass
real :: ax(1:np), ay(1:np), az(1:np), bx(1:np), by(1:np), bz(1:np)
real :: cx(5), cy(5), cz(5), dx(2), dy(2), dz(2)
real :: dist1, dist2, dist, min_dist, rem_dist
real :: ix(1:2*np), iy(1:2*np), iz(1:2*np)

real :: max_dist
! idxb        = index of point before first segment
! idx         = index of point to be moved (in the middle) (first of multiple points)
! idxa        = index of point after second segment
! jdx         = index counter
! dir_sign    = sign of coefexp
! idx1        = idx + 1 or 1 if crossing polygon end/beginning
! idx2        = idx1 + 1 or 1 if crossing polygon end/beginning
! jdx1        = jdx + 1 or 1 if crossing polygon end/beginning
! n_first_dup = number of duplicate points when idx = 1
! n_dup       = number of duplicate points between two valid segments
! pts_to_move = count of points left to move
! ptb         = point before the point to be moved
! pt          = point to be moved
! pta         = point after point to be moved
! vgb         = great circle vector of segment before point
! vg          = direction vector for moving point (great circle vector of perpendicular segment)
! vga         = great circle vector of segment after point
! vgn         = vector with magnitude to move point
! ptn         = point after move
! ax, ay, az  = polygon array placed in CCW polygon order (temporary use)
! bx, by, bz  = polygon array placed in CCW polygon order (temporary use)
! cx, cy, cz  = five point arrays in x, y, z for checking (and finding) line crossing.
! dx, dy, dz  = two point arrays in x, y, z for distance checking
! ix, iy, iz  = arrays in x, y, z for intersection points
! dist1, dist2 = distance between first (and second) point and intersection point in radians
! dist        = distance between two points in radians
! min_dist    = minimum distance to point intersection

logical find_cross
real arc_len

! find_cross = function to check (and find) great circle crossing
! arc_len    = function to find great circle distance between two points in radians

! debugging
real clat(5), clon(5), sc
sc = atan(1.)/45.

! polygon direction is that from before to after with the inside to the left.

! place polygon into working array in native order
do idx = 1,np
    ax(idx) = x(list(idx))
    ay(idx) = y(list(idx))
    az(idx) = z(list(idx))
end do

! move points toward interior or exterior (depending on sign of COEFEXP)
! movement is along great circle from b(efore) to a(fter)
rem_dist = coefexp
dir_sign = sign(1.0,coefexp)

! debug check
!write(*,*) "#coefexp, dist1"
!do idx = 1, np
!    dx(1) = ax(idx)
!    dy(1) = ay(idx)
!    dz(1) = az(idx)
!    dx(2) = bx(idx)
!    dy(2) = by(idx)
!    dz(2) = bz(idx)
!    dist1 = arc_len(dx, dy, dz)
!   write(*,*) coefexp, dist1
!end do
!stop

! check for locations where the arcs of point movement cross.
!crossing = .false.
crossing = .true.
n_pass = 400
do while( crossing )
    ! move points
    call move_points(rem_dist, np, ax, ay, az, bx, by, bz, ier)
    n_pass = n_pass - 1
    if(n_pass .eq. 0 ) then
        ! place new polygon into old polygon's place
        do idx = 1,np
            ax(idx) = bx(idx)
            ay(idx) = by(idx)
            az(idx) = bz(idx)
        end do
        exit
    end if
    ! prove there is a crossing to continue
    crossing = .false.
    ! a very large number
    min_dist = 1.0e10
    do idx = 1,np
        ! check each line from old point to new point against every other created line
        cx(1) = ax(idx)
        cy(1) = ay(idx)
        cz(1) = az(idx)
        cx(2) = bx(idx)
        cy(2) = by(idx)
        cz(2) = bz(idx)
        
        ! check against all other points
        do jdx = idx+1,np
            cx(3) = ax(jdx)
            cy(3) = ay(jdx)
            cz(3) = az(jdx)
            cx(4) = bx(jdx)
            cy(4) = by(jdx)
            cz(4) = bz(jdx)
            if( find_cross (cx,cy,cz) ) then
                ! set flag that crossing occured
                crossing = .true.
                ! it crosses, find distance to crossing
                ! first point to intersection
                dx(1) = cx(1)
                dy(1) = cy(1)
                dz(1) = cz(1)
                dx(2) = cx(5)
                dy(2) = cy(5)
                dz(2) = cz(5)
                dist1 = arc_len(dx, dy, dz)
                ! second point to intersection
                dx(1) = cx(3)
                dy(1) = cy(3)
                dz(1) = cz(3)
                dist2 = arc_len(dx, dy, dz)
                dist = min(dist1, dist2)
                if( dist .lt. min_dist ) then
                    ! this is the minimum crosser
                    min_dist = dist
                    n_cross = 1
                    cross_list(1) = idx
                    cross_list(2) = jdx
                    ix(1) = cx(5)
                    iy(1) = cy(5)
                    iz(1) = cz(5)
                else if( dist .eq. min_dist ) then
                    ! tied with the minimum crosser
                    n_cross = n_cross + 1
                    !write(*,*) "# idx, jdx, n_cross, min_dist:", idx, jdx, n_cross, min_dist
                    cross_list(n_cross*2-1) = idx
                    cross_list(n_cross*2) = jdx
                    ix(n_cross) = cx(5)
                    iy(n_cross) = cy(5)
                    iz(n_cross) = cz(5)
                end if
            end if
        end do
    end do

    !write(*,*) "# Ready to move Minimum distance"
    ! check for crosser and distance to move points
    if( crossing ) then
        ! two movement arcs cross, move all points the minimum distance
        call move_points(dir_sign*min_dist, np, ax, ay, az, bx, by, bz, ier)
        !write(*,*) "# points moved"
        ! then move the crossing points to common point
        do idx = 1,n_cross
            ! check for crossing from np to 1
            if( (cross_list(idx*2)-cross_list(idx*2-1)) .gt. np/2 ) then
                ! indexes go too far, must be wrapped, run across break (end to beginning
                ! from beginning to first point
                do jdx = 1, cross_list(idx*2-1)
                    bx(jdx) = ix(idx)
                    by(jdx) = iy(idx)
                    bz(jdx) = iz(idx)
                end do
                ! from second point to end
                do jdx = cross_list(idx*2), np
                    bx(jdx) = ix(idx)
                    by(jdx) = iy(idx)
                    bz(jdx) = iz(idx)
                end do
            else
                ! loop from first point to second point getting all points in between
                do jdx = cross_list(idx*2-1), cross_list(idx*2)
                    bx(jdx) = ix(idx)
                    by(jdx) = iy(idx)
                    bz(jdx) = iz(idx)
                end do
            end if
        end do
        rem_dist = rem_dist - dir_sign*min_dist
        if( sign(1.0,rem_dist) .ne. dir_sign ) then
            ! moved points past input distance. This is an error
            !write(*,*) "# n_pass, coefexp, min_dist, rem_dist:", n_pass, coefexp, min_dist, rem_dist
            write(*,*) "# ERROR: interpolate:polysizing: past input distance"
            stop
        end if
    else
        ! no crossing, so move went the entire distance
        rem_dist = 0.0
    end if

    !write(*,*) "# Checking for self intersection"
    ! check for self intersection of new polygon and remove smaller loops
    do idx = 1,np - 4
        ! check each line from old point to new point against every other created line
        if( idx .eq. np ) then
            idx1 = 1
        else
            idx1 = idx + 1
            if( idx1 .eq. np ) then
                idx2 = 1
            end if
            idx2 = idx1 + 1
        end if
        cx(1) = bx(idx)
        cy(1) = by(idx)
        cz(1) = bz(idx)
        cx(2) = bx(idx1)
        cy(2) = by(idx1)
        cz(2) = bz(idx1)
        
        ! check against all other segments that can intersect
        do jdx = idx2,np-2
            if( jdx .eq. np ) then
                jdx1 = 1
            else
                jdx1 = jdx + 1
            end if
            cx(3) = bx(jdx)
            cy(3) = by(jdx)
            cz(3) = bz(jdx)
            cx(4) = bx(jdx1)
            cy(4) = by(jdx1)
            cz(4) = bz(jdx1)
            if( find_cross (cx,cy,cz) ) then
                ! found self intersection, count around and find shortest loop
                ! since idx starts at 1, jdx1 will always be greater than idx
                ! move all points in short loop to intersection point
                if( (jdx1 - idx) .le. (idx + np - jdx1) ) then
                    ! short loop in positive index direction
                    do kdx = idx1, jdx
                        ! move all points in short loop to intersection point
                        bx(kdx) = cx(5)
                        by(kdx) = cy(5)
                        bz(kdx) = cz(5)
                    end do
                else
                    ! short loop crosses break in polygon
                    do kdx = idx, 1, -1
                        ! move all points in short loop to intersection point
                        bx(kdx) = cx(5)
                        by(kdx) = cy(5)
                        bz(kdx) = cz(5)
                    end do
                    do kdx = np, jdx1, -1
                        ! move all points in short loop to intersection point
                        bx(kdx) = cx(5)
                        by(kdx) = cy(5)
                        bz(kdx) = cz(5)
                    end do
                end if
            end if
        end do
    end do    

    ! place new polygon into old polygon's place
    do idx = 1,np
        ax(idx) = bx(idx)
        ay(idx) = by(idx)
        az(idx) = bz(idx)
    end do

    ! check for maximum distance moved by any one point
    do idx = 1,np
        dx(1) = x(list(idx))
        dy(1) = y(list(idx))
        dz(1) = z(list(idx))
        dx(2) = bx(idx)
        dy(2) = by(idx)
        dz(2) = bz(idx)
        if( (coefexp .lt. 0.0) .and. (-arc_len(dx,dy,dz) .lt. coefexp) ) then
            ! one point moved greater than maximum distance, return for incremental mirroring
            ier = -1  ! indicates early termination
            crossing = .false.  ! used to exit outer loop
            exit  ! use to exit inner loop
        end if
    end do

    ! debug
    n_pass = n_pass - 1
    if(n_pass .eq. 0 ) then
        exit
    end if

end do

!write(*,*) "# Number of Passes left (n_pass):", n_pass

! copy new points back into new array
do idx = 1,np
    nx(list(idx)) = bx(idx)
    ny(list(idx)) = by(idx)
    nz(list(idx)) = bz(idx)
end do

! this set value left to move polygon for early return use
coefexp = rem_dist

! write(*,*) "# before return:", n_pass

return
end

subroutine move_points (coefexp, np, x, y, z, nx, ny, nz, ier)
real coefexp
integer np, ier
real :: x(1:np), y(1:np), z(1:np), nx(1:np), ny(1:np), nz(1:np)

!***********************************************************
!   This subroutine takes an counter clock wise ordered input polygon
! (inside to the left) and moves the points either in or out depending
! on the coefficient. No checks are done for points crossing ie. self
! intersection.

! On input:
!       COEFEXP = distance the polygon is to shrink in fraction of unit sphere radius
!       NP = number of points in polygon array
!       X,Y,Z = Cartesian coordinates of points of input polygon

! Input parameters are not altered by this function.

! On output:
!       NX,NY,NZ = Cartesian coordinates of points of resulting polygon

!       IER = Error indicator:
!             IER = 0 if no errors were encountered.
!             IER = 1 polygon with only one valid segment
!             IER = 2 polygon with no valid segments
!***********************************************************

integer :: idxb, idx, idxa, jdx, n_first_dup, n_dup, pts_to_move
integer :: dup_list(np)
real :: ptb(3), pt(3), pta(3), vgb(3), vg(3), vga(3), vgn(3), ptn(3)

! idxb        = index of point before first segment
! idx         = index of point to be moved (in the middle) (first of multiple points)
! idxa        = index of point after second segment
! n_first_dup = number of duplicate points when idx = 1
! n_dup       = number of duplicate points between two valid segments
! pts_to_move = count of points left to move
! ptb         = point before the point to be moved
! pt          = point to be moved
! pta         = point after point to be moved
! vgb         = great circle vector of segment before point
! vg          = direction vector for moving point (great circle vector of perpendicular segment)
! vga         = great circle vector of segment after point
! vgn         = vector with magnitude to move point
! ptn         = point after move

! polygon direction is that from before to after with the inside to the left.

! move points toward interior or exterior (depending on sign of COEFEXP
! loop through all bounding polygon line segments
! move pt to ptn by moving it along the line bisecting the angle between the two segments
! negative values for coefexp move to the left (shrinking the polygon)
! positive values for coefexp move to the right (expanding the polygon)

n_first_dup = 0
n_dup = 0
pts_to_move = np
idxb = np
idx = 1
idxa = 2

! find first valid segment
do while(.TRUE.)
    ! assign polygon point to vector
    ptb(1) = x(idxb)
    ptb(2) = y(idxb)
    ptb(3) = z(idxb)

    pt(1) = x(idx)
    pt(2) = y(idx)
    pt(3) = z(idx)

    ! previous side, in  counter clock wise direction is great circle bisecting
    ! angle between lower indexed polygon segment and this segment.
    ! great circle of previous segment
    call cross_product(ptb, pt, vgb)
    call unit_vector(vgb, vgb, ier)
    if( ier .gt. 0 ) then
        ! the two points were identical, point cannot be moved until a direction is found
        ! this is the first point of the polygon, so move backward to find valid segment
        n_first_dup = n_first_dup + 1
        idxb = idxb - 1
        if( idxb .lt. 1 ) then
            !failed to find a valid first segment, exit with error
            ier = 2
            return
        end if
    else
        exit
    end if
end do

! valid first segment found find second valid segment
do while (pts_to_move .gt. 0)
    ! set index for point after
    if( idx .eq. np ) then
        idxa = 1
    else
        idxa = idx + 1
    end if
    ! assign polygon point to vector
    pta(1) = x(idxa)
    pta(2) = y(idxa)
    pta(3) = z(idxa)

    ! great circle of following segment
    call cross_product(pt, pta, vga)
    call unit_vector(vga, vga, ier)
    if( ier .gt. 0 ) then
        ! the two points were identical, point cannot be moved until a direction is found
        ! move forward to find valid segment
        n_dup = n_dup + 1
        idx = idx + 1
        idxa = idxa + 1
        if( idxa .gt. np ) then
            ! no valid segment was found, degenerate polygon
            ier = 1 
            exit
        end if
        cycle
    end if

    ! direction to move point (sum two vectors and normalize)
    ! add two segment normal unit vectors to get direction
    call add_vector(vgb, vga, vg)
    call unit_vector(vg, vg, ier)
    if( ier .gt. 0 ) then
        ! vectors are in exactly the opposite direction, degenerate polygon.
        ! no valid segment was found, degenerate polygon
        ier = 1 
        exit
    end if

    ! find new point location (see sign convention above)
    call smult_vector(tan(-coefexp), vg, vgn)
    call add_vector(pt, vgn, ptn)
    call unit_vector(ptn, ptn, ier)

    ! transfer to new polygon arrays
    nx(idx) = ptn(1)
    ny(idx) = ptn(2)
    nz(idx) = ptn(3)
    ! point moved, decrement count
    pts_to_move = pts_to_move - 1

    ! check for duplicate points and move them also
    do while( n_first_dup .gt. 0 )
        ! transfer to new polygon arrays
        jdx = np - n_first_dup + 1
        nx(jdx) = ptn(1)
        ny(jdx) = ptn(2)
        nz(jdx) = ptn(3)
        ! point moved, decrement counts
        n_first_dup = n_first_dup - 1
        pts_to_move = pts_to_move - 1
    end do

    ! check for duplicate points and move them also
    do while( n_dup .gt. 0 )
        ! transfer to new polygon arrays
        jdx = idx - n_dup
        nx(jdx) = ptn(1)
        ny(jdx) = ptn(2)
        nz(jdx) = ptn(3)
        ! point moved, decrement counts
        n_dup = n_dup - 1
        pts_to_move = pts_to_move - 1
    end do

    ! move indexes and previously calculated values to next pair of segments
    ! present value is now before point
    idxb = idx
    ptb = pt
    ! advance value is now point to be moved
    idx = idxa
    pt = pta
    ! great circle vector for these two points has already been calculated
    vgb = vga

end do
return
end

logical function find_cross (px,py,pz)
real    px(5), py(5), pz(5)

! This subroutine takes four points, the first two of which describe a
! great circle arc A, and points 3 and 4 which describe a great circle
! arc B. The arc is the short distance, not the long distance. The
! function returns TRUE if the two arcs have an intersection point and 
! FALSE if they do not. The intersection point is returned.

! On input:
!       PX,PY,PZ = Arrays of length 5 containing the Cartesian
!           coordinates of the four points to be tested
!       These are not altered by the routine

! On return:

!       TRUE if the two arcs have an intersection point
!            and the intersection point is placed in P?(5)
!       FALSE if they do not
!***********************************************************

integer idx, ier
real    vp1(3), vp2(3), vp3(3), vp4(3), vf(3), vg(3), vi(3)
logical lft1, lft2

! Local parameters:

! put four intial points in to vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)
vp3(1) = px(3)
vp3(2) = py(3)
vp3(3) = pz(3)
vp4(1) = px(4)
vp4(2) = py(4)
vp4(3) = pz(4)

! find the great circle normal vector for the first two points
call cross_product(vp1, vp2, vg)
call unit_vector(vg, vg, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    find_cross = .FALSE.
    return
end if

! find the great circle normal vector for the last two points
call cross_product(vp3, vp4, vf)
call unit_vector(vf, vf, ier)
if( ier .gt. 0 ) then
    ! the two points were identical
    find_cross = .FALSE.
    return
end if

! check for position of first two points with relation to great circle formed by last two points
lft1 = (dot_product(vp1, vf) .gt. 0.0)
lft2 = (dot_product(vp2, vf) .gt. 0.0)
if( lft1 .eqv. lft2 ) then
    ! similar sign of the dot product means that both vectors are on the same side of the great circle.
    find_cross = .FALSE.
    return
end if

! check for position of last two points with relation to great circle formed by first two points
lft1 = (dot_product(vp3, vg) .gt. 0.0)
lft2 = (dot_product(vp4, vg) .gt. 0.0)
if( lft1 .eqv. lft2 ) then
    ! similar sign of the dot product means that both vectors are on the same side of the great circle.
    find_cross = .FALSE.
    return
end if

! if we reached here, then the two arcs must intersect
find_cross = .TRUE.
call cross_product(vf, vg, vi)
! check for right hemisphere
if( dot_product(vi, vp1) .lt. 0.0 ) then
    ! reverse hemisphere of intersection
    call cross_product(vg, vf, vi)
end if
call unit_vector(vi, vi, ier)
if( ier .gt. 0 ) then
    ! the points are really colinear or the same (degenerate case)
    find_cross = .FALSE.
end if

! copy intersection into return location
px(5) = vi(1)
py(5) = vi(2)
pz(5) = vi(3)

return
end

real function arc_len (px,py,pz)
real    px(2), py(2), pz(2)

! This subroutine takes two points on the surface of the unit sphere and
! finds the distance between them, expressed in radians. This can be
! converted into a great circle distance if the radius is known.

! On input:
!       PX,PY,PZ = Arrays of length 2 containing the Cartesian
!           coordinates of the two points to be tested
!       These are not altered by the routine

! On return:

!       the arc angle in radian from point 1 to point 2
!***********************************************************

real    vp1(3), vp2(3), result

! Local parameters:

! put two points into vector form
vp1(1) = px(1)
vp1(2) = py(1)
vp1(3) = pz(1)
vp2(1) = px(2)
vp2(2) = py(2)
vp2(3) = pz(2)

! find the arc angle for the first two points
result = dot_product(vp1,vp2)
if( (result .ge. 1.0) .or. (result .le. -1.0) ) then
    arc_len = 0.0
else
    arc_len = acos( dot_product(vp1,vp2))
end if

return
end