!
! 
! List of subroutines contained in Orblib.f90 in alphabetical order:
! Most of them are rewritten from PB's ORB2.bas or ORB3.bas
! 
!     ABG2xy
!     ABG2LonLat
!     Atan2f
!     Checkflipped
!     Cross
!     Drawnode
!     DrawElement
!     DrawFault
!     Exists
!     Flipped
!     GetElement
!     Modeparms
!     OnAnyE
!     Nearest_
!     NEXTto
!     OnAnyF
!     Pixels
!     RaiseOne
!     Scaler
!     SetBins
!     Shade
!     Sterradians
!     XandY
!     xy2abg
!     xy2lonlat
!     unitize
!     UnitVec
!     winframe
!
! Following subroutines modified from orbweave
!
! Convert any unit vector (Alpha, Beta, Gamma) to the (x,y) system of current window
! Alpha: points to lat=0, lon=0
! Beta:  points to lat=0, lon=90E
! Gamma: points to north pole
! (x,y) system of current window is independent of window rotation
subroutine ABG2xy(tempvec, visible, x, y)
! input: tempvec(3)
! output: visible, x,y
! require: global winout(3)
use global
real tempvec(3)
real x,y
real dot
integer visible
dot = tempvec(1)*winout(1) + tempvec(2)*winout(2) + tempvec(3)*winout(3)
if (dot < 0) then
   visible = 0
   x = 0
   y = 0
else
   visible = 1
   x = tempvec(1)*winright(1) + tempvec(2)*winright(2) + tempvec(3)*winright(3)
   y = tempvec(1)*winup(1) + tempvec(2)*winup(2) + tempvec(3)*winup(3)
endif
end subroutine ABG2xy

! Convert any vector in (Alpha, Beta, Gamma) to N. Longitude and E. Latitude, in Radians
! Alpha: points to lat=0, lon=0
! Beta:  points to lat=0, lon=90E
! Gamma: points to north pole
! where (alpha^2 + beta^2 +gamma^2 = 1)
subroutine ABG2lonlat(tempvec, lon, lat)
! input: tempvec(3)
! output: lon, lat
real tempvec(3)
real lon, lat
real equapart
lon = atan2f(tempvec(2), tempvec(1))
equapart = sqrt(tempvec(1)*tempvec(1) + tempvec(2)*tempvec(2))
lat = atan2f(tempvec(3), equapart)
end subroutine ABG2lonlat

! arc-tangent function with two arguments, returning result from -pi to pi
real function atan2f(y,x)
! input: x,y
real x,y
if((y==0.0).and.(x==0.0)) then
    atan2f = 0.0
elseif (x==0.0) then
    if (y > 0) then
	   atan2f = 1.570795
	else
	   atan2f = -1.570795
	end if
else
    atan2f = atan2(y,x)
endif
return
end function atan2f

! Check flipped of all elements, if element is inverted, reverse
! nodes(2, n), nodes(3, n), where n is element#
! this routine is called following Make_Global_Grid, which seems to give
! elements not in consistent counterclockwise order
! 
subroutine CheckFlipped
! Adjust: EMemo, nodes
use global
integer(4) i, tmp
do i = 1, NUMEL
   call flipped(i)
   if (EMemo(i)) then
       tmp = nodes(2, i)
	   nodes(2,i) = nodes(3, i)
	   nodes(3,i) = tmp
	   EMemo(i) = 0
   end if
end do
end subroutine

! cross product of two vectors
subroutine cross(v1, v2, vcross)
real v1(3), v2(3), vcross(3)
vcross(1) = v1(2)*v2(3) - v1(3)*v2(2)
vcross(2) = v1(3)*v2(1) - v1(1)*v2(3)
vcross(3) = v1(1)*v2(2) - v1(2)*v2(1)
end subroutine cross

! Draw single node
subroutine drawnode(icolor,n)
! input: n -- node #
! require: global ifrontcolor,R2C, nodeABG, subroutine (ABG2xy, pixels)
use global
use dflib
real tempvec(3)
real x, y
integer visible
integer col, row
integer(4) n
integer dummy, col0, row0, col1, row1
integer(4) iret,icolor
!! important !!
iret = SETCOLORRGB(icolor)
!!
tempvec(1) = nodeABG(1,n)
tempvec(2) = nodeABG(2,n)
tempvec(3) = nodeABG(3,n)
call ABG2xy (tempvec, visible, x, y)
if (visible.AND.ShowNodes) then
   if ((x*x + y*y) < R2C) then
      call pixels(x,y, col, row, visible)
      if (visible > 0) then
	     ! circle radius: 3 pixels
         col0 = col - 3
	     row0 = row - 3
	     col1 = col + 3
	     row1 = row + 3
	     dummy = ELLIPSE($GBORDER, col0, row0, col1, row1)
      end if
   endif
endif
end subroutine drawnode

! Draw one element (n)
! 
subroutine drawelement(icolor,n)
! input: icolor, n(element#)
! require: global switches (ColorIn, contour);
!          ifrontcolor,nodeABG, nodes, eqcm, subroutine (ABG2xy,setcolorRGB, moveto,lineto)
!          shade!!!
use global
use dflib
integer(4) n,k1,k2,k3
integer visible
integer ecol(3), erow(3)
integer ColC, RowC
integer(4) j, jtmp
real tempvec(3), R21, R22, R23
real xt(3),yt(3)
real f1, f2, f3
!
real flow, fhigh
integer ilow, ihigh
integer Clist(10), Rlist(10)  ! temporary array, be cautious that size may be exceeded!!!
integer npnts, numcon, index
integer iin, iout
real xlow, ylow, xmid, ymid, xhigh, yhigh
real fcon, fmid, sbase, sfar
real xbase, ybase, xfar, yfar
real xobase, yobase, xofar, yofar
integer inner, inold
!
integer(2) status
integer(4) iret,icolor
type (xycoord) xy
!****************************
!INTEGER(2) style
! style = #FFFF ! solid line
! style = #EEEE ! dashed line; 11/16 bits drawn
! style = #AAAA ! dotted line
! CALL SETLINESTYLE(style)
!N.B. Doesn't work well, because same element side is common to two
!     elements and is drawn twice, with dashes not lining up.
!     A better way to keep the basemap visible is to draw it
!     last (on top), not before.
!***************************
!
iret = setcolorRGB(icolor)
!
k1 = nodes(1,n)
tempvec(1) = nodeABG(1,k1)
tempvec(2) = nodeABG(2,k1)
tempvec(3) = nodeABG(3,k1)
call ABG2xy(tempvec, visible, xt(1),yt(1))
if (visible > 0) then
else
  return
end if
R21 = xt(1)*xt(1) + yt(1)*yt(1)
k2 = nodes(2,n)
tempvec(1) = nodeABG(1,k2)
tempvec(2) = nodeABG(2,k2)
tempvec(3) = nodeABG(3,k2)
call ABG2xy(tempvec, visible, xt(2), yt(2))
if (.not.visible) return
R22 = xt(2)*xt(2) + yt(2)*yt(2)
k3 = nodes(3,n)
tempvec(1) = nodeABG(1,k3)
tempvec(2) = nodeABG(2,k3)
tempvec(3) = nodeABG(3,k3)
call ABG2xy(tempvec, visible, xt(3),yt(3))
if (visible > 0) then
else
   return
end if
R23 = xt(3)*xt(3) + yt(3)*yt(3)
if ((R21 < R2C).or.(R22 < R2C).or.(R23 < R2C)) then
   do j = 1, 3
        call pixels(xt(j),yt(j), ecol(j),erow(j), visible)
        if (.not.visible) return
   end do
   ! for gap viewing purpose
   if(ColorIn) then
        call flipped(n)
        if(EMemo(n) == 1) then
		   jtmp = ecol(2)
		   ecol(2) = ecol(3)
		   ecol(3) = jtmp
		   jtmp = erow(2)
		   erow(2) = erow(3)
		   erow(3) = jtmp
		end if
        call RaiseOne(ecol, erow)
   end if
   ! draw contour given option
   if(contour) then
      if(idata > 6) then
	      ! add mus here
	  else
	      f1 = eqcm(idata, k1)
		  f2 = eqcm(idata, k2)
		  f3 = eqcm(idata, k3)
	  end if
	  if(logs) then
	      if(f1>0) then
		     f1 = log10(f1)
		  else
		     f1 = botF
		  end if
		  if(f2>0) then
		     f2 = log10(f2)
		  else
		     f2 = botf
		  end if
		  if(f3>0) then
		     f3 = log10(f3)
		  else
		     f3 = botf
		  end if
	  end if
	  flow = f1
	  if(f2 < flow) flow = f2
	  if(f3 < flow) flow = f3
	  fhigh = f1
	  if(f2 > fhigh) fhigh = f2
	  if(f3 > fhigh) fhigh = f3
      IF (dfc > 0.0) THEN
	      ilow = int(2.01 + (flow - botf)/dfc)
	      ihigh = int(1.99 + (fhigh - botf)/dfc)
      ELSE
          ilow  = 1
          ihigh = 1
      END IF
	  if((ilow >= ihigh).or.(ilow > hicolor).or.(ihigh < 2)) then
	      ! uniform color, no contours
		  Clist(1:3) = ecol(1:3)
		  Rlist(1:3) = erow(1:3)
		  npnts = 3
		  index = ihigh
		  call shade(npnts, ecol, erow, Clist, Rlist, index)
	  else
	      ! contours must be located
		  numcon = ihigh - ilow
		  ! rotate to standard view ??
		  if(f1==flow) then
		       xlow = ecol(1)
			   ylow = erow(1)
			   if(f2<=f3) then
			        fmid = f2
					xmid = ecol(2)
					ymid = erow(2)
					xhigh = ecol(3)
					yhigh = erow(3)
			   else
			        fmid = f3
					xmid = ecol(3)
					ymid = erow(3)
					xhigh = ecol(2)
					yhigh = erow(2)
			   end if
		  elseif(f2==flow) then
		       xlow = ecol(2)
			   ylow = erow(2)
			   if(f1<=f3) then
			        fmid = f1
					xmid = ecol(1)
					ymid = erow(1)
					xhigh = ecol(3)
					yhigh = erow(3)
			   else
			        fmid = f3
					xmid = ecol(3)
					ymid = erow(3)
					xhigh = ecol(1)
					yhigh = erow(1)
			   end if
		  else
		       ! f3 is the low value
			   xlow = ecol(3)
			   ylow = erow(3)
			   if(f1<=f2) then
			        fmid = f1
					xmid = ecol(1)
					ymid = erow(1)
					xhigh = ecol(2)
					yhigh = erow(2)
			   else
			        fmid = f2
					xmid = ecol(2)
					ymid = erow(2)
					xhigh = ecol(1)
					yhigh = erow(1)
			   end if
	      end if
		  do j = 1, numcon
		       ! iin > = 2 !!
		       iin = ilow + j - 1
			   iout = iin + 1
			   fcon = botf + (iin-1)*dfc
			   sbase = (fcon - flow)/(fhigh - flow)
			   xbase = int(xlow + sbase * (xhigh - xlow))
			   ybase = int(ylow + sbase * (yhigh - ylow))
			   if(fcon < fmid) then
			       inner = 1
               else
			       inner = 0
			   end if
			   if(inner) then
			       sfar = (fcon - flow)/(fmid - flow)
				   xfar = int(xlow + sfar*(xmid-xlow))
				   yfar = int(ylow + sfar*(ymid-ylow))
			   else
				   sfar = (fcon - fmid)/(fhigh - fmid)
				   xfar = int(xmid + sfar*(xhigh - xmid))
				   yfar = int(ymid + sfar*(yhigh - ymid))
			   end if
			   if(j==1) then  ! 1st polygon
			       Clist(1) = xlow
				   Rlist(1) = ylow
				   Clist(2) = xbase
				   Rlist(2) = ybase
				   Clist(3) = xfar
				   Rlist(3) = yfar
				   if (inner) then
				        npnts = 3
				   else
				        npnts = 4
						Clist(4) = xmid
						Rlist(4) = ymid
				   end if
				   index = iin
				   call shade(npnts, ecol, erow, Clist, Rlist, index)
			   else
			       ! one more polygon in a series
				   Clist(1) = xobase
				   Rlist(1) = yobase
				   Clist(2) = xbase
				   Rlist(2) = ybase
				   Clist(3) = xfar
				   Rlist(3) = yfar
				   if(inner == inold) then
				        npnts = 4
						Clist(4) = xofar
						Rlist(4) = yofar
				   else
				        npnts = 5
						Clist(4) = xmid
						Rlist(4) = ymid
						Clist(5) = xofar
						Rlist(5) = yofar
				   end if
				   index = iin
				   call shade(npnts, ecol, erow, Clist, Rlist, index)
			   end if
			   if(j==numcon) then
			       ! last polygon to fill in
				   Clist(1) = xbase
				   Rlist(1) = ybase
				   Clist(2) = xhigh
				   Rlist(2) = yhigh
				   if(inner) then
				        npnts = 4
						Clist(3) = xmid
						Rlist(3) = ymid
						Clist(4) = xfar
						Rlist(4) = yfar
				   else
				        npnts = 3
						Clist(3) = xfar
						Rlist(3) = yfar
				   end if
				   index = iout
				   call shade(npnts, ecol, erow, Clist, Rlist, index)
			   end if
			   xobase = xbase
			   yobase = ybase
			   xofar = xfar
			   yofar = yfar
			   inold = inner
		  end do ! loop of numcon
      end if ! determine contour
   end if ! contour
   ! draw element side
   iret = SETCOLORRGB(icolor)
   call moveto(ecol(1), erow(1), xy)
   status = lineto(ecol(2),erow(2))
   call moveto(ecol(2), erow(2), xy)
   status = lineto(ecol(3),erow(3))
   call moveto(ecol(3), erow(3), xy)
   status = lineto(ecol(1), erow(1))
   ! check flipped status of element, update EMemo
   call flipped(n)
   ! add a schematic triangle at the center point
   if (doicon) then
	  ColC = sum(ecol) / 3.0
      RowC = sum(erow) / 3.0
      if (.not.EMemo(n)) then
	     ecol(1) = ColC - 2*sqrt(3.0)
         erow(1) = RowC + 2
         ecol(2) = ColC + 2*sqrt(3.0)
         erow(2) = RowC + 2
         ecol(3) = ColC
         erow(3) = RowC - 4
         call moveto(ecol(1), erow(1), xy)
         status = lineto(ecol(2),erow(2))
         call moveto(ecol(2), erow(2), xy)
         status = lineto(ecol(3),erow(3))
         call moveto(ecol(3), erow(3), xy)
         status = lineto(ecol(1), erow(1))
      else
	     ecol(1) = ColC - 2*sqrt(3.0)
		 erow(1) = RowC - 2
		 ecol(2) = ColC + 2*sqrt(3.0)
		 erow(2) = RowC + 2
		 call moveto(ecol(1), erow(1), xy)
		 status = lineto(ecol(2),erow(2))
		 ecol(1) = ColC + 2*sqrt(3.0)
		 erow(1) = RowC - 2
		 ecol(2) = ColC - 2*sqrt(3.0)
		 erow(2) = RowC + 2
		 call moveto(ecol(1), erow(1), xy)
		 status = lineto(ecol(2), erow(2))
      end if
   end if ! icon at the center of element
end if
! restore back to default color
!iret = SETCOLORRGB(ifrontcolor)
!****************************
!style = #FFFF ! solid line
! CALL SETLINESTYLE(style)
!N.B. See comments above
!***************************
end subroutine drawelement

!
! draw single fault
! dipping symbols varies with dip angle
!
! Modified on March 10, 2005, to show correct dip symbols following
! Bird and Kagan [2004] table 5:
! normal_dip_degree = 55.0; thrust_dip_degree = 25.0; subduction_dip_degree = 14.
! 
! Steep fault using tick: changed from > 55.0 to >= 55.0
! shallow dip: filled if <= 18.0, changed from if <=27.5
!
subroutine drawfault(icolor,n)
! input: n(fault#)
! require: global ifrontcolor, NODEF, nodeABG, FDIP
use dflib
use dfwin
use global
implicit none
integer(4) n
integer(4) k1,k2
integer i
real tempvec(3)
real xt(2),yt(2)
real R21,R22
real dx, dy, size
real s, f1, f2, dip, adip, df1ds, df2ds, dxds, dyds, arg, normal, xp, yp
real x,y
real atan2f
integer visible
integer ecol(2), erow(2)
integer colA, rowA, colB, rowB, colC, rowC, colfill, rowfill
type (xycoord) xy
!!!!!!!!!!!!!!!!!!!!!!!!!
type (xycoord) poly(4)
!!!!!!!!!!!!!!!!!!!!!!!!!
integer(2) status
integer(4) icolor,iret
!
iret = SETCOLORRGB(icolor) 
!
k1 = nodef(1,n)
tempvec(1) = nodeABG(1,k1)
tempvec(2) = nodeABG(2,k1)
tempvec(3) = nodeABG(3,k1)
call ABG2xy(tempvec, visible, xt(1),yt(1))
if (visible > 0) then
else
   return
end if
R21 = xt(1)*xt(1) + yt(1)*yt(1)
k2 = nodef(2,n)
tempvec(1) = nodeABG(1,k2)
tempvec(2) = nodeABG(2,k2)
tempvec(3) = nodeABG(3,k2)
call ABG2xy(tempvec, visible, xt(2),yt(2))
if (.not.visible) return
R22 = xt(2)*xt(2) + yt(2)*yt(2)
if ((R21 > R2C).and.(R22 > R2C)) return
call pixels(xt(1),yt(1),ecol(1), erow(1), visible)
if (.not.visible) return
call pixels(xt(2),yt(2),ecol(2), erow(2), visible)
if (.not.visible) return
call moveto(ecol(1),erow(1),xy)
status = lineto(ecol(2),erow(2))
if ((abs(ecol(1)-ecol(2)) > abs(erow(1)-erow(2)))) then
    call moveto(ecol(1), erow(1) + 1, xy)
	status = lineto(ecol(2),erow(2) + 1)
	call moveto(ecol(1), erow(1) - 1, xy)
	status = lineto(ecol(2),erow(2) - 1)
else
    call moveto(ecol(1) + 1, erow(1), xy)
	status = lineto(ecol(2) + 1,erow(2))
	call moveto((ecol(1) - 1), erow(1), xy)
	status = lineto(ecol(2) - 1,erow(2))
end if
! plot dip symbols
dx = xt(2) - xt(1)
dy = yt(2) - yt(1)
size = 0.1 * sqrt(dx*dx + dy*dy)
do i = 1, 4
  s = FIPoint(i)
  f1 = 1 - s
  f2 = s
  dip = fdip(1,n) * f1 + fdip(2,n) * f2
  adip = 90.0 - abs(dip - 90)
  if (adip < 76.0) then
      x = f1*xt(1) + f2*xt(2)
	  y = f1*yt(1) + f2*yt(2)
	  df1ds = -1.0
	  df2ds = 1.0
	  dxds = df1ds * xt(1) + df2ds * xt(2)
	  dyds = df1ds * yt(1) + df2ds * yt(2)
	  arg = atan2f(dyds, dxds)
	  if (dip < 90.0) then
	     normal = arg - 1.5714
	  else
	     normal = arg + 1.5714
	  end if
	  ! if (adip > 55.0) then ! steep fault use tick
	  ! changed to >= 55.0
	  if(adip >= 55.0) then
	     call pixels(x, y, colA, rowA, visible)
         if (.not.visible) return
		 xp = x + size * cos(normal)
		 yp = y + size * sin(normal)
		 call pixels(xp,yp, colB, rowB, visible)
         if (.not.visible) return
		 call moveto(colA,rowA,xy)
		 status = lineto(colB, rowB)
	  elseif (adip > 35.0) then ! intermediate dip; use open box
	     xp = x + 0.5 * size * cos(arg)
		 yp = y + 0.5 * size * sin(arg)
		 call pixels(xp,yp, colA, rowA, visible)
         if (.not.visible) return
		 xp = xp + size * cos(normal)
		 yp = yp + size * sin(normal)
		 call pixels(xp, yp, colB, rowB, visible)
         if (.not.visible) return
		 call moveto(colA, rowA, xy)
		 status = lineto(colB, rowB)
		 xp = xp - size * cos(arg)
		 yp = yp - size * sin(arg)
		 call pixels(xp, yp, colA, rowA, visible)
         if (.not.visible) return
		 status = lineto(colA, rowA)
		 xp = xp - size * cos(normal)
		 yp = yp - size * sin(normal)
		 call pixels(xp, yp, colA, rowA, visible)
         if (.not.visible) return
		 status = lineto(colA, rowA)
	  else ! shallow dip; triangular (WAS: "filled if <= 27.5 degrees" in OrbWeaver and Shells, up to 2004.12.)
	       ! (Now, changed to "filled if <= 19 degrees" for OrbWin and Shells 2005+, per Bird & Kagan [2004, BSSA].)
	     xp = x + 0.5 * size * cos(arg)
		 yp = y + 0.5 * size * sin(arg)
		 call pixels(xp, yp, colA, rowA, visible)
         if (.not.visible) return
		 call moveto(colA, rowA, xy)
		 xp = x + size * cos(normal)
		 yp = y + size * sin(normal)
		 call pixels(xp, yp, colB, rowB, visible)
         if (.not.visible) return
		 status = lineto(colB, rowB)
		 xp = x - 0.5 * size * cos(arg)
		 yp = y - 0.5 * size * sin(arg)
		 call pixels(xp, yp, colC, rowC, visible)
         if (.not.visible) return
		 status = lineto(colC, rowC)
		 status = lineto(colA, rowA)
		 ! if (adip <= 27.5) then ! subduction zone, fill interior of tick
		   if (adip <= 19.0) then ! 
		    poly(1)%xcoord = colA
			poly(1)%ycoord = rowA
		    poly(2)%xcoord = colB
			poly(2)%ycoord = rowB
		    poly(3)%xcoord = colC
			poly(3)%ycoord = rowC
		    poly(4)%xcoord = colA
			poly(4)%ycoord = rowA
!			colfill = (colA + colB + colC) / 3
!			rowfill = (rowA + rowB + rowC) / 3
!           iret = floodfillRGB(colfill, rowfill, icolor)
            status = POLYGON($GFILLINTERIOR, poly, INT2(3))
         end if
      end if
  end if 
end do ! next i
! restore to default graphic color
! iret = SETCOLORRGB(ifrontcolor)
end subroutine drawfault

! Check whether a node already exists at location, if so, return node number, otherwise
! return 0. Test begin with n1 to allow repeat calls. n1 = 1, search nodes 1 to numnod
! otherwise, search n1 to numnod. This may speed up the search. 
subroutine exists(n1, x, y, n)
! input: n1, x, y (location x,y)
! output: node# n
! require: nodeABG, NUMNOD, Tolerance
use global
real x,y
real R2min, Rmin, dA, dB, dG, R2
real tempvec(3)
integer outside
integer(4) n1,n, i,iclose
call xy2ABG(x,y, outside, tempvec)
if (outside == 1) then
    n = 0
else
    R2min = 3.E+10
	do i = n1, numnod
	   dA = tempvec(1) - nodeABG(1,i)
	   dB = tempvec(2) - nodeABG(2,i)
	   dG = tempvec(3) - nodeABG(3,i)
	   R2 = dA*dA + dB*dB + dG*dG
	   if (R2 < R2min) then
	       R2min = R2
		   iclose = i
       else
	   end if
    end do
    Rmin = sqrt(R2min)
	if (Rmin <= tolerance) then
	   n = iclose
    else
	   n = 0
	end if
end if
end subroutine exists


! Check flip or not of an element n
subroutine flipped(n)
! input: nodes, nodeABG
! output: EMemo
use global
real v12(3), v23(3), v31(3), vcross(3)
integer(4) n, n1, n2, n3
real areasign
n1 = nodes(1, n)
n2 = nodes(2, n)
n3 = nodes(3, n)
v12(1) = nodeABG(1, n2) - nodeABG(1, n1)
v12(2) = nodeABG(2, n2) - nodeABG(2, n1)
v12(3) = nodeABG(3, n2) - nodeABG(3, n1)
v23(1) = nodeABG(1, n3) - nodeABG(1, n2)
v23(2) = nodeABG(2, n3) - nodeABG(2, n2)
v23(3) = nodeABG(3, n3) - nodeABG(3, n2)
v31(1) = nodeABG(1, n1) - nodeABG(1, n3)
v31(2) = nodeABG(2, n1) - nodeABG(2, n3)
v31(3) = nodeABG(3, n1) - nodeABG(3, n3)
call cross(v12, v23, vcross)
areasign = vcross(1) * nodeABG(1,n1) + vcross(2) * nodeABG(2, n1) + vcross(3) *nodeABG(3, n1)
if (areasign > 0) then
   EMemo(n) = 0  ! no problem with element
else
   EMemo(n) = 1  ! there is problem with element
end if       
end subroutine flipped



! Find triangular element closest to (x, y), ported from basic script
subroutine GetElement(x, y, jj, R2min)
! input: x, y
! output: jj(element#), R2min
use Global
integer(4) jj, i
integer(4) n1, n2, n3
integer visible
real x, y, x1, y1, x2, y2, x3, y3
real xc, yc, dx, dy, R2
real R2min
real tempvec(3), tempv2(3), tempv3(3)
R2min = 3.3E+10
do i = 1, NUMEL
   n1 = Nodes(1,i)
   n2 = Nodes(2,i)
   n3 = Nodes(3,i)
   tempvec(1) = nodeABG(1, n1)
   tempvec(2) = nodeABG(2, n1)
   tempvec(3) = nodeABG(3, n1)
   call ABG2xy(tempvec, visible, x1, y1)
   if (visible > 0) then
      tempv2(1) = nodeABG(1, n2)
	  tempv2(2) = nodeABG(2, n2)
	  tempv2(3) = nodeABG(3, n2)
	  call ABG2xy(tempv2, visible, x2, y2)
	  if (visible > 0) then
	      tempv3(1) = nodeABG(1, n3)
		  tempv3(2) = nodeABG(2, n3)
		  tempv3(3) = nodeABG(3, n3)
		  call ABG2xy(tempv3, visible, x3, y3)
		  if (visible > 0) then
		     xc = (x1 + x2 + x3)/3.0
			 yc = (y1 + y2 + y3)/3.0
			 dx = x - xc
			 dy = y - yc
			 R2 = dx*dx + dy*dy
			 if (R2 < R2min) then
			    R2min = R2
				jj = i
			 end if
		  end if
	  end if
   end if
end do
end subroutine GetElement


! Table parameters describing current screen mode
! Work on current window
! output: global HiRow, HiCol, Lines
subroutine modeparms
use dflib
use global
type (windowconfig) wc
logical status
status = getwindowconfig(wc)
HiRow = wc%numypixels
HiCol = wc%numxpixels
Lines = wc%numtextrows
end subroutine modeparms

! Given node# n, search range (n1 to nlast), find corresponding element and vertice#
subroutine OnAnyE(n, n1, nlast, ie, je)
!input: n, n1, nlast
!output: ie (element), je(vertice)
use global
integer(4) n, n1, nlast
integer(4) ie, je
integer(4) i, j
ie = 0
je = 0
do i = n1, nlast
   do j = 1, 3
      if (nodes(j, i) == n) then
	      ie = i
		  je = j
		  return
      end if
   end do
end do
return
end subroutine OnAnyE 

!
! Find nearest nodes close to mouse location (up to 4 nodes)
!
subroutine Nearest_(x, y, numnear, nearones)
! input: x, y
! output: numnear, nearones(4)
! require: numnod, nodeABG(numnod)
use global
real x, y
real tempvec(3)
real R2toN, R2
real AN, BN, GN
integer numnear
integer outside
integer(4) i
integer(4) nearones(4)
call xy2ABG(x, y, outside, tempvec)
if(outside) then
   numnear = 0
   nearones(1:4) = 0 
else
   R2toN = 9.99E+10
   do i = 1, NUMNOD
      R2 = (tempvec(1) - nodeABG(1, i))**2 + (tempvec(2) - nodeABG(2, i))**2 + &
	       (tempvec(3) - nodeABG(3, i))**2 
      if(R2 < R2toN) then
           AN = nodeABG(1, i)
		   BN = nodeABG(2, i)
		   GN = nodeABG(3, i)
		   R2toN = R2   
	  end if
   end do
   numnear = 0
   Nearones(1:4) = 0
   do i = 1, NUMNOD
       if(nodeABG(1, i) == AN) then
             if(nodeABG(2, i) == BN) then
			    if(nodeABG(3, i) == GN) then
				     numnear = numnear + 1
					 nearones(numnear) = i
				end if
 			 end if
		end if
		if(numnear == 4) goto 1000
   end do
end if
1000 return
end subroutine Nearest_

!
! determine whether there are more element adjacent to side j of triangular continuum element i
!     j = 1 --- the side opposite node # nodes(1, i)
!     j = 2 --- the side opposite node # nodes(2, i)
!     j = 3 --- the side opposite node # nodes(3, i)
! if a fault element is adjacent, its number is kfault, otherwise, kfault = 0
! if another triangular continuum element is adjacent (even across fault element kfault), then
! its number is kele, otherwise, kele = 0
!
subroutine NEXTto(i, j, Kfault, jf, Kele, je)
! input: i(element#), j(side#)
! output: Kfault, jf, Kele, je
use global
integer(4) n1, n2
integer(4) k
integer(4) m1, m2, m3, m4
integer(4) i, Kele, Kfault, Klow, Khigh
integer j, jf, je
logical foundf
! two node numbers along the side of interest, counterclockwise
n1 = nodes(mod(j,3)+1,i)
n2 = nodes(mod(j+1,3)+1,i)
! check for adjacent fault element first
foundf = .false.
do k = 1, NFL
   m1 = nodef(1, k)
   m2 = nodef(2, k)
   m3 = nodef(3, k)
   m4 = nodef(4, k)
   if ((m1 == n2).and.(m2 == n1)) then
       foundf = .true.
	   Kfault = k
	   jf = 1
	   goto 1
   end if
   if ((m3 == n2).and.(m4 == n1)) then
       foundf = .true.
	   Kfault = k
	   jf = 2
	   goto 1
   end if
end do
1   if(.not.foundf) Kfault = 0
    ! if there is fault, replace 2 node numbers that we search for
	if(foundf) then
	    if (jf == 1) then
		    n1 = m3
			n2 = m4
		else ! jf = 2 
		    n1 = m1
			n2 = m2
		end if
	end if
	! search for adjacent triangular continuum element
	Kele = 0
	Klow = i
	Khigh = i
2	Klow = Klow - 1
	if (Klow > 0) then
	    do k = 1, 3
		   m1 = nodes(mod(k,3)+1, Klow)
		   m2 = nodes(mod(k+1,3)+1,Klow)
		   if ((m1 == n2).and.(m2==n1)) then
		       Kele = Klow
			   je = k
			   return
		   end if
		end do
	end if
	Khigh = Khigh + 1
	if (Khigh <= NUMEL) then
	    do k = 1, 3
		     m1 = nodes(mod(k,3)+1, Khigh)
			 m2 = nodes(mod(k+1,3)+1,Khigh)
			 if((m1==n2).and.(m2==n1)) then
			    Kele = Khigh
				je = k
				return
			 end if
		end do
	end if
if((Klow > 1).or.(Khigh < NUMEL)) goto 2
end subroutine NEXTto

!
! Given node# n, search range (n1, to nlast), locate any fault including the specificed node n
!
subroutine OnAnyF(n, n1, nlast, ie, je)
! input: n, n1, nlast
! output: ie(fault element), je(vertice)
use Global
integer(4) n, n1, nlast, ie, je
integer(4) i, j
ie = 0
je = 0
do i = n1, nlast
   do j = 1, 4
      if (nodef(j, i) == n) then
	      ie = i
		  je = j
		  return
	  end if
   end do
end do 
end subroutine OnAnyF

! convert (x,y) coordinates to pixels; prevent overflows
subroutine pixels(x,y, col, row, visible)
! input: x, y
! output: col, row
! modify: visible (usually 1 = T on entry, but may switch to 0 = F)
! require: global scales (2,2)
use global
real x, y, xcol, yrow  ! note: real(8) require 1.D0 format
integer col, row, visible
xcol = scales(1,1)*x + scales(1,2)
yrow = scales(2,1)*y + scales(2,2)
if ((abs(xcol) < 30000).and.(abs(yrow) < 30000)) then
   !safe to convert real to integer, and call C-routines (for integers up to ~32000):
   visible = 1
   col = xcol
   row = yrow
else
   visible = 0
   col = 0
   row = 0
end if
end subroutine pixels

!
! subroutine for raising color-code of each pixel
! Note: General logic follows PB basic script. But 
!       SETPIXELRGB & GETPIXELRGB are used!!
! 
subroutine RaiseOne(vcol, vrow)
! Require: Hirow, Hicol
use dfwin
use dflib
use global
real x(4), y(4), r
integer vcol(3), vrow(3)
integer row1, row2, col1, col2
integer row, col, edge
integer i,j, index
integer(4) hue, istatus
logical usewhich(3,2)
row1 = 3.2e+6
row2 = -3.2e+6
do j = 1, 3
   if(vrow(j) > row2) row2 = vrow(j)
   if(vrow(j) < row1) row1 = vrow(j)
end do
if(row1 < 0) row1 = 0
if(row1 > hirow) row1 = hirow
if(row2 < 0) row2 = 0
if(row2 > hirow) row2 = hirow
if(row2 <= row1) return
if(vrow(1) /= vrow(2)) then
   if(vrow(2) > vrow(1)) then
      usewhich(1,1) = .true.
	  usewhich(1,2) = .false.
   else
      usewhich(1,1) = .false.
	  usewhich(1,2) = .true.
   end if
else
   usewhich(1,1) = .false.
   usewhich(1,2) = .false.
end if
if(vrow(2) /= vrow(3)) then
   if(vrow(3) > vrow(2)) then
      usewhich(2,1) = .true.
	  usewhich(2,2) = .false.
   else
      usewhich(2,1) = .false.
	  usewhich(2,2) = .true.
   end if
else
   usewhich(2,1) = .false.
   usewhich(2,2) = .false.
end if
if(vrow(3) /= vrow(1)) then
   if(vrow(1) > vrow(3)) then
      usewhich(3,1) = .true.
	  usewhich(3,2) = .false.
   else
      usewhich(3,1) = .false.
	  usewhich(3,2) = .true.
   end if
else
   usewhich(3,1) = .false.
   usewhich(3,2) = .false.
end if
x(1:3) = real(vcol(1:3))
x(4) = x(1)
y(1:3) = real(vrow(1:3))
y(4) = y(1)
do row = row1, row2 -1
   r = real(row)
   col1 = 0
   col2 = hicol
   do i = 1,3
      j = i + 1
	  if(usewhich(i,1)) then
	      edge = 1 + int(x(i) + (x(j)-x(i))*(r-y(i))/(y(j)-y(i)))
		  if(edge > col1) col1 = edge
	  end if
	  if(usewhich(i,2)) then
	      edge = int(x(i) + (x(j)-x(i))*(r-y(i))/(y(j)-y(i)))
		  if(edge < col2) col2 = edge
	  end if
   end do
   do col = col1, col2
      hue = GETPIXELRGB(col, row)
	  do k = 1, hicolor
	     if(hue == colorarray(k)) then
		    index = k
			goto 1
         else
		    index = 0
		 end if
      end do
1     if(index == Hicolor) then
         index = 1
      else
         index = index + 1
      end if
      istatus = SETPIXELRGB(col, row, colorarray(index))
   end do
end do
end subroutine Raiseone

!
! Uses WindowHeight, HiCol, and HiRow (from MODULE Global) to compute
! R2C, Scales(3), and UnScale(3) (in global) which
! convert back and forth from (x,y) to (col,row). 
!
subroutine scaler
use global
R2C = (0.5 * WindowHeight) ** 2 + (0.666667 * WindowHeight) **2
scales(1, 1) = 0.75 * real(HiCol) / WindowHeight
scales(1, 2) = 0.5 * real(HiCol)
scales(2, 1) = -real(HiRow) / WindowHeight
scales(2, 2) = 0.5 * real(HiRow)
unscale(1, 1) = 1.33333 * WindowHeight / real(HiCol)
unscale(1, 2) = -0.66667 * WindowHeight
unscale(2, 1) = -WindowHeight / real(HiRow)
unscale(2, 2) = 0.5 * WindowHeight
End subroutine scaler

!
! Figure out best assignment of (Hicolor -1) colors to nodal data values
! logs indicates that base-10 logarithmic scale should be used
! if idata = 1, 2, 3, 4, 5, 6 then data is taken from eqcm(idata, nodes)
! if idata = 7, 8, 9, then data is taken from mus(idata-6, elements)
!
!
subroutine setbins
! input: contour, idata
! modify: oldbotf, oldidata, oldtopf
! output: botf, dfc, logs, topf
! require: hicolor, eqcm, numnod
use global
logical got_one, visible
real maxf, minf, x, y
REAL, DIMENSION(3) :: tempvec
integer(4) col, i, n, row
if(idata > 6) then
   ! put mus here
else ! Determine range of data VISIBLE IN THE CURRENT WINDOW,
     ! and use that range as the basis for dividing the color bar.
     !(Actually, what is currently programmed is to get the 
     ! range from the nodal values visible in the window,
     ! neglecting possible linear gradients to more extreme
     ! nodal values just slightly outside the window.)
!GPBhere
    got_one = .FALSE.
    do n = 1, numnod
        tempvec(1) = nodeABG(1,n)
        tempvec(2) = nodeABG(2,n)
        tempvec(3) = nodeABG(3,n)
        call ABG2xy (tempvec, visible, x, y)
        if (visible > 0) then ! this JUST means it is in the near hemisphere!
           !still need to check coordinates (x, y) against window frame.
            CALL pixels(x, y, col, row, visible)
            IF ((col >= 0).AND.(col <= HiCol)) THEN
                IF ((row >= 0).AND.(row <= HiRow)) THEN
                    IF (got_one) THEN ! this is NOT the 1st node in the window
                        maxf = MAX(maxf, eqcm(idata, n))
                        minf = MIN(minf, eqcm(idata, n))
                    ELSE ! this IS the first node found in the window
                        maxf = eqcm(idata, n)
                        minf = maxf
                    END IF
                    got_one = .TRUE.
                END IF ! row is within window
            END IF ! column is within window
        end if ! visible
    end do ! loop on all nodes
    IF (.NOT.got_one) THEN ! empty window (no nodes yet, or extreme zoom)
        maxf = 0.0
        minf = 0.0
    END IF
end if ! idata > 6, or NOT

! automatically let colors span the data range
if((minf>0).and.(maxf<.001)) then
    !use log scale to bring out range of mu_ (Restore v.2):
     logs = .true.
	 botf = log10(minf)
	 topf = log10(maxf)
	 dfc = (topf-botf)/real(hicolor - 1)
else
   !use linear scale, allowing negative values:
     logs = .false.
	 botf = minf
	 topf = maxf
end if
!Note: It is not a problem if topf = botf = 0.0
dfc = (topf-botf)/real(hicolor - 1)
! NeedToDraw omitted here
oldidata = idata
oldbotf = botf
oldtopf = topf
end subroutine setbins

!
! shading subroutine used in DrawElement
! Note: using Polygon fill to speed things up!
!       Different way!
!
subroutine shade(npnts, ecol, erow, Clist, Rlist, index)
! input: npnts, Clist(10), Rlist(10), ecol(3), erow(3), index
use dflib
use global
integer(4) index, temp_index
integer(4) npnts
integer(4) Clist(10), Rlist(10)
integer(4) ecol(3), erow(3)
real seedC, seedR, R2, area,t
real longest, wide
integer(4) iret
integer(4) jt, jp
integer(4) c1, c2, rw1, rw2
integer(4) Hue, icolor
integer(2) status
type(xycoord) xy
type(xycoord), allocatable:: poly(:)
! for out of range, set to lowest color
! contour use color (2-Hicolor)
temp_index = MAX(1, MIN(HiColor, index))
!!
colorpick = colorarray(temp_index)
iret = setcolorRGB(colorpick)
!
seedC = 0.
seedR = 0.
allocate(poly(npnts))
do jt = 1, npnts
   seedC = seedC + real(Clist(jt))
   seedR = seedR + real(Rlist(jt))
   if(jt < npnts) then
       jp = jt + 1
   else
       jp = 1
   end if
   c1 = Clist(jt)
   c2 = Clist(jp)
   rw1 = Rlist(jt)
   rw2 = Rlist(jp)
   call moveto(c1, rw1, xy)
   status = lineto(c2,rw2)
   poly(jt)%xcoord = c1
   poly(jt)%ycoord = rw1
end do
seedC = seedC/real(npnts)
seedR = seedR/real(npnts)
Hue = GETPIXELRGB(int(seedC), int(seedR))
!! for efficiency purpose, use following
status = POLYGON($GFILLINTERIOR, poly, int2(npnts))
deallocate(poly)
end subroutine shade

!
! Calculate solid angle of spherical triangle in sterradians
! v1, v2, v3: unit vectors from center of unit sphere to three vertices
! if v1, v2, v3: counterclockwise, area is positive
!                clockwise,                negative
! ported from PB basic script
!
subroutine sterradians(v1, v2, v3, SolidAngle)
real v1(3), v2(3), v3(3)
real(8) vcross(3), ra(3), rb(3), rc(3), cb(3), ba(3), center(3)
real(8) la, lb, lc, adot, across, aside, bdot, bcross, bside, cdot, ccross, cside
real(8) s, denominator, sinhalf, coshalf, half, a, b, c
real(8) dot
real SolidAngle
ra(1:3) = DBLE(v1(1:3))
rb(1:3) = DBLE(v2(1:3))
rc(1:3) = DBLE(v3(1:3))
! ensure they are truly unit vectors after double precision
la= DSQRT(ra(1)**2 + ra(2)**2 + ra(3)**2)
ra(1:3) = ra(1:3)/la
lb = DSQRT(rb(1)**2 + rb(2)**2 + rb(3)**2)
rb(1:3) = rb(1:3)/lb
lc = DSQRT(rc(1)**2 + rc(2)**2 + rc(3)**2)
rc(1:3) = rc(1:3)/lc
! compute length of triangle side
! a side
adot = rb(1)*rc(1) + rb(2)*rc(2) + rb(3)*rc(3)
vcross(1) = rb(2)*rc(3) - rb(3)*rc(2)
vcross(2) = rb(3)*rc(1) - rb(1)*rc(3)
vcross(3) = rb(1)*rc(2) - rb(2)*rc(1)
across = DSQRT(vcross(1)**2 + vcross(2)**2 + vcross(3)**2)
if(across == 0.D0) then
   SolidAngle = 0.
   return
end if
! aside = atan2d(across, adot)
aside = DBLE(atan2f(real(across), real(adot)))
! bside
bdot = rc(1)*ra(1) + rc(2)*ra(2) + rc(3)*ra(3)
vcross(1) = rc(2)*ra(3) - rc(3)*ra(2)
vcross(2) = rc(3)*ra(1) - rc(1)*ra(3)
vcross(3) = rc(1)*ra(2) - rc(2)*ra(1)
bcross = DSQRT(vcross(1)**2 + vcross(2)**2 + vcross(3)**2)
if(bcross == 0.D0) then
   SolidAngle = 0.
   return
end if
! bside = atan2d(bcross, bdot)
bside = DBLE(atan2f(real(bcross), real(bdot)))
! cside
cdot = ra(1)*rb(1) + ra(2)*rb(2) + ra(3)*rb(3)
vcross(1) = ra(2)*rb(3) - ra(3)*rb(2)
vcross(2) = ra(3)*rb(1) - ra(1)*rb(3)
vcross(3) = ra(1)*rb(2) - ra(2)*rb(1)
ccross = DSQRT(vcross(1)**2 + vcross(2)**2 + vcross(3)**2)
if(ccross == 0.D0) then
   SolidAngle = 0.
   return
end if
! cside = atan2d(ccross, cdot)
cside = DBLE(atan2f(real(ccross), real(cdot)))
! ---------------------------------------------------------------
! determine half-angle, then double them
 s = (aside + bside + cside)/2.D0
 denominator = sin(bside)*sin(cside)
 sinhalf = dsqrt(dsin(s-bside)*dsin(s-cside)/denominator)
 coshalf = dsqrt(dsin(s)*dsin(s-aside)/denominator)
 half = DBLE(atan2f(real(sinhalf), real(coshalf)))
 a = 2.0D0*half
 !
 denominator = sin(cside)*sin(aside)
 sinhalf = dsqrt(dsin(s-cside)*dsin(s-aside)/denominator)
 coshalf = dsqrt(dsin(s)*dsin(s-bside)/denominator)
 half = DBLE(atan2f(real(sinhalf), real(coshalf)))
 b = 2.0D0*half
 !
 denominator = sin(aside)*sin(bside)
 sinhalf = dsqrt(dsin(s-aside)*dsin(s-bside)/denominator)
 coshalf = dsqrt(dsin(s)*dsin(s-cside)/denominator)
 half = DBLE(atan2f(real(sinhalf), real(coshalf)))
 c = 2.0D0*half
! ---------------------------------------------------------------
SolidAngle = real((a + b + c) - 3.1415926535D0)
! check sign change
! first, get outward normal vector usign cross product
ba(1:3) = rb(1:3)- ra(1:3)
cb(1:3) = rc(1:3) - rb(1:3)
vcross(1) = ba(2)*cb(3) - ba(3)*cb(2)
vcross(2) = ba(3)*cb(1) - ba(1)*cb(3)
vcross(3) = ba(1)*cb(2) - ba(2)*cb(1)
! second, get rough center positon by averaging nodes
center(1:3) = ra(1:3) + rb(1:3) + rc(1:3)
! last, use product to check their general directions
dot = center(1)*vcross(1) + center(2)*vcross(2) + center(3)*vcross(3)
if(dot < 0.D0) SolidAngle = - SolidAngle
end subroutine sterradians


! convert pixel coordinates to (x,y) using predefined transformation
subroutine XandY (col, row, x,y)
! input: col, row
! output: x,y
! require:  global array unscale(2,2)
use global
integer col, row
real x,y
x = unscale(1,1) * col + unscale(1,2)
y = unscale(2,1) * row + unscale(2,2)
end subroutine XandY


! 
! Draws a line
!

! convert (x,y) of orthographic projection plane to (alpha, beta, gamma)
! orthogonal cartesian coordinates
subroutine xy2abg(x,y, outside, tempv)
! input: x, y
! output: outside, tempv
! Require: global array winright(3), winup(3), winout(3)
use global
real x, y, r2, outofplane
real tempv(3)
integer outside
r2 = x**2 + y**2
if (r2 >= 1.0) then
  outside = 1  
  tempv(1) = 0.
  tempv(2) = 0.
  tempv(3) = 0.
else
  outside = 0
  outofplane = sqrt(1 - r2)
  tempv(1) = x*winright(1) + y*winup(1) + outofplane*winout(1)
  tempv(2) = x*winright(2) + y*winup(2) + outofplane*winout(2)
  tempv(3) = x*winright(3) + y*winup(3) + outofplane*winout(3)
endif
end subroutine xy2abg


! Converts (x,y) in projection plane to (lat, lon) if x^2 + y^2 < 1
! else declare outside
subroutine xy2lonlat(x,y,outside, lon, lat)
! input: x,y
! output: outside, lon,lat
! require: ABG2lonlat, global
use global
real x,y
real lon, lat
integer outside
real r2, outofplane
real tempvec(3), tempv2(3)
r2= x*x + y*y
if (r2>=1) then
   outside = 1
   lon = 0.
   lat = 0.
else
   outside = 0
   outofplane = sqrt(1-r2)
   tempvec(1) = x*winright(1) + y*winup(1) + outofplane*winout(1)
   tempvec(2) = x*winright(2) + y*winup(2) + outofplane*winout(2)
   tempvec(3) = x*winright(3) + y*winup(3) + outofplane*winout(3)
   call ABG2lonlat(tempvec,lon,lat)
endif
end subroutine xy2lonlat

! Convert any vector in 3-component double-precision form to a unit vector
subroutine unitize(x,y,z)
! modify: x, y, z
real(8) x, y, z
real(8) r
r = dsqrt(x*x + y*y + z*z)
if (r > 0.D0) then
   x = x/r
   y = y/r
   z = z/r
else
   x = 1.D0
   y = 0.D0
   z = 0.D0
endif
end subroutine unitize

!
!!Converts any vector to a unit vector.
!
subroutine UnitVec(tempvec)
real r2, r
real tempvec(3)
R2 = tempvec(1) * tempvec(1) + tempvec(2) * tempvec(2) + tempvec(3) * tempvec(3)
if (r2 > 0) Then
        r = sqrt(r2)
		tempvec(1:3) = tempvec(1:3)/r
else
        tempvec(1) = 1.
		tempvec(2) = 0.
		tempvec(3) = 0.
end If
end subroutine Unitvec


! determine unit vectors of window frame 
!
subroutine winframe(winlat, winlon, winright, winout, winup)
! input: winlat, winlon
! output: winright, winout, winup
! require: subroutine cross
real winlat, winlon
real winright(3), winout(3), winup(3)
winout(1) = cos(winlat) * cos(winlon)
winout(2) = cos(winlat) * sin(winlon)
winout(3) = sin(winlat)
winright(1) = cos(winlon + 1.570796)
winright(2) = sin(winlon + 1.570796)
winright(3) = 0.0
call cross(winout, winright, winup)
end subroutine winframe