! ! 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 implicit none real, intent(in) :: tempvec(3) logical, intent(out) :: visible real, intent(out) :: x, y real dot dot = tempvec(1)*winout(1) + tempvec(2)*winout(2) + tempvec(3)*winout(3) if (dot < 0) then visible = .FALSE. x = 0 y = 0 else visible = .TRUE. 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 implicit none real, intent(in) :: tempvec(3) real, intent(out) :: lon, lat real equapart real ATan2F 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 implicit none real, intent(in) :: y, x ! The input "tangent" for this "arctangent" function would be (y / x) if it were a single REAL. if ((y == 0.0).and.(x == 0.0)) then ATan2F = 0.0 else if (x == 0.0) then if (y > 0) then ATan2F = 1.57079632 else ATan2F = -1.57079632 end if else ATan2F = ATan2(y,x) endif return end function ATan2F ! Check for flipped (reversed) condition 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 implicit none 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) implicit none real, intent(in) :: v1(3), v2(3) real, intent(out) :: 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 implicit none integer(4), intent(in) :: iColor, n real tempvec(3) real x, y logical visible integer col, row integer dummy, col0, row0, col1, row1 integer(2) dummy_2, col0_2, row0_2, col1_2, row1_2 integer(4) iRet !! 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) then ! circle radius: 3 pixels col0 = col - 3 row0 = row - 3 col1 = col + 3 row1 = row + 3 col0_2 = col0; row0_2 = row0; col1_2 = col1; row1_2 = row1 ! required INT(2) arguments dummy_2 = ELLIPSE($GBORDER, col0_2, row0_2, col1_2, row1_2) end if endif endif end subroutine DrawNode ! Draw one element (n) ! subroutine DrawElement(iColor, n) ! Input: iColor {for outline and icon; not for any "contour" filling}, n (element#) ! Requires: Global switches (colorIn, contour); ! ifrontcolor,nodeABG, nodes, eqcm, subroutine (ABG2xy,setcolorRGB, moveto,lineto) ! Shade!!! use global use dflib implicit none integer(4), intent(in) :: iColor, n integer(4) k1, k2, k3 logical visible integer ecol(3), erow(3), eleCol(5), eleRow(5) integer ColC, RowC integer(4) j, jtmp integer(2) status_2, arg1_2, arg2_2 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, which probably never has more than 5 entries 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 LOGICAL inner, inOld, pentagon_somewhere integer(2) status integer(4) iRet 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 (.not.visible) return 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 (.not.visible) return 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/overlap-viewing purposes: 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(s) if this logical switch is on: IF (contour) THEN IF (iData > 6) THEN ! add mus here ELSE IF (iData >= 1) THEN ! use eqcm array f1 = eqcm(iData, k1) f2 = eqcm(iData, k2) f3 = eqcm(iData, k3) ELSE IF (iData == 0) THEN ! color depends on continuum_LRi(n) f1 = continuum_LRi(n) f2 = f1 f3 = f1 END IF if (logs) then if (f1 > 0.0) then f1 = log10(f1) else f1 = botF end if if (f2 > 0.0) then f2 = log10(f2) else f2 = botf end if if (f3 > 0.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.0 + (fLow - botF) / dFC) iHigh = INT(2.0 + (fHigh - botF) / dFC) ELSE iLow = 2 ! because color index #1 = gray is not part of the color spectrum iHigh = 2 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 else if (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 pentagon_somewhere = (fMid > fLow).AND.(fHigh > fMid) inOld = .TRUE. 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)) inner = (fCon < fMid) if (inner) then IF (fMid > fLow) THEN sFar = (fCon - fLow) / (fMid - fLow) ELSE sFar = 0.0 END IF xFar = INT(xLow + sFar*(xMid-xLow)) yFar = INT(yLow + sFar*(yMid-yLow)) else IF (fHigh > fMid) THEN sFar = (fCon - fMid) / (fHigh - fMid) ELSE sFar = 0.0 END IF 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 (pentagon_somewhere.AND.(inner.NEQV.inOld)) THEN npnts = 5 cList(4) = xMid rList(4) = yMid cList(5) = xofar rList(5) = yofar ELSE ! normal case npnts = 4 cList(4) = xofar rList(4) = 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! j = 1, numCon END IF ! this element has internal contours END IF ! contour (color-in elements?) ! draw element sides: iRet = SETCOLORRGB(iColor) arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(3); arg2_2 = erow(3) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(3); arg2_2 = erow(3) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! check flipped status of element, update EMemo (.TRUE. if flipped) call flipped(n) ! add a schematic triangle (or X) at the center point if (doicon) then ColC = sum(ecol) / 3.0 RowC = sum(erow) / 3.0 if (.not.EMemo(n)) then ! normal icon(s) = (older) small unfilled triangle, or (newer) color-filled square IF (plot_eleCenter_icons) THEN ! in Global; updated by Redraw IF (appType == 3) THEN ! special element-center icon for Restore3+ mode: displays mu(s) as color(s)! ! define limits of square (in terms of pixel rows/columns): eleCol(1) = ColC - 6 eleRow(1) = RowC + 6 eleCol(2) = ColC + 6 eleRow(2) = eleRow(1) eleCol(3) = eleCol(2) eleRow(3) = RowC - 6 eleCol(4) = ColC - 6 eleRow(4) = eleRow(3) eleCol(5) = eleCol(1) eleRow(5) = eleRow(1) IF ((element_data(2, n) > 9999.0D0).OR.(element_data(1, n) == element_data(3, n))) THEN ! !produce ONE square of uniform color: cList(1:5) = eleCol(1:5) rList(1:5) = eleRow(1:5) nPnts = 5 index = NINT(2.0D0 + 13.0D0 * (element_data(1, n) - lowest_mu) / (highest_mu - lowest_mu)) IF (iColor == iBackColor) index = 1 ! gray; used for deleting an existing element CALL Shade (nPnts, eCol, eRow, cList, rList, index) ELSE ! the box must be divided into two colors: younger mu on left, older mu on right: !produce LEFT rectangle of color corresponding to younger mu: cList(1:5) = eleCol(1:5) rList(1:5) = eleRow(1:5) cList(2:3) = ColC nPnts = 5 index = NINT(2.0D0 + 13.0D0 * (element_data(1, n) - lowest_mu) / (highest_mu - lowest_mu)) IF (iColor == iBackColor) index = 1 ! gray; used for deleting an existing element CALL Shade (nPnts, eCol, eRow, cList, rList, index) !GPBhere !produce RIGHT rectangle of color corresponding to older mu: cList(1:5) = eleCol(1:5) rList(1:5) = eleRow(1:5) cList(1) = ColC cList(4:5) = ColC nPnts = 5 index = NINT(2.0D0 + 13.0D0 * (element_data(3, n) - lowest_mu) / (highest_mu - lowest_mu)) IF (iColor == iBackColor) index = 1 ! gray; used for deleting an existing element CALL Shade (nPnts, eCol, eRow, cList, rList, index) END IF ! draw the box around the color(s) iRet = SETCOLORRGB(iColor) arg1_2 = eleCol(1); arg2_2 = eleRow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) ! do not draw while going to initial position arg1_2 = eleCol(2); arg2_2 = eleRow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! draw line to point2 arg1_2 = eleCol(3); arg2_2 = eleRow(3) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! draw line to point3 arg1_2 = eleCol(4); arg2_2 = eleRow(4) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! draw line to point4 arg1_2 = eleCol(5); arg2_2 = eleRow(5) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! draw line to point5 = point1 ELSE ! appType /= 3 (probably 1 = thin-shell, or 2 = Restore2 mode) :: draw the standard (old-school) triangle ecol(1) = ColC - 3.464 ! -2*sqrt(3.0), precomputed for speed erow(1) = RowC + 2 ecol(2) = ColC + 3.464 ! +2*sqrt(3.0), precomputed for speed erow(2) = RowC + 2 ecol(3) = ColC erow(3) = RowC - 4 iRet = SETCOLORRGB(iColor) arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(3); arg2_2 = erow(3) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(3); arg2_2 = erow(3) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) END IF ! special mu-colored icon? Or, ordinary green triangle? END IF ! plot_eleCenter_icons? ELSE ! BAD element icon = X ecol(1) = ColC - 5.196 ! -3*sqrt(3.0), precomputed for speed erow(1) = RowC - 3 ecol(2) = ColC + 5.196 ! +3*sqrt(3.0), precomputed for speed erow(2) = RowC + 3 arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ecol(1) = ColC + 5.196 ! +3*sqrt(3.0), precomputed for speed erow(1) = RowC - 3 ecol(2) = ColC - 5.196 ! -3*sqrt(3.0), precomputed for speed erow(2) = RowC + 3 arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) end if end if ! icon at the center of element end if end subroutine DrawElement ! ! Color-selecting routine used for the cores of fault traces ! when LRDraw = .TRUE. ! INTEGER*4 FUNCTION iColorFromLRi (LRi) USE Global ! for botF, dFC, hiColor IMPLICIT NONE INTEGER*4, INTENT(IN) :: LRi REAL :: f INTEGER*4 :: index f = 1.0 * LRi ! convert to floating-point IF (dFC /= 0.0) THEN index = INT(2.0 + (f - botF) / dFC) ! NOTE that index #1 denotes dark-gray, index = MIN(hiColor, MAX(2, index)) ! which is not considered part of the color spectrum. ELSE index = 2 END IF iColorFromLRi = colorArray(index) END FUNCTION iColorFromLRi ! ! Draw single fault element; ! dip symbols vary 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#) ! Requires: Global's nodeF, nodeABG, fDip use dflib use dfwin use global implicit none integer(4), intent(in) :: iColor, n integer(4) i, k1, k2, red_I4, yellow_I4 integer(2) status_2, arg1_2, arg2_2 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 logical visible integer ecol(2), erow(2) integer colA, rowA, colB, rowB, colC, rowC, colfill, rowfill type (xycoord) xy !!!!!!!!!!!!!!!!!!!!!!!!! type (xycoord) poly(4) !!!!!!!!!!!!!!!!!!!!!!!!! integer(4) iRet !-------------------------- iRet = SETCOLORRGB(iColor) !This single function-call is usually all we need, and produces a fault trace (triple-width) ! with dip symbols that are all the same color. !However, note that special plotting convention for fault elements in the ! Lithospheric Rheology (LR#) map modes (either single-element, or block-set) ! require another SETCOLORRGB call to produce a variable-color trace line bordered by red, ! so as to set the trace apart from any neighboring continuum elements of the same LR#. !-------------------------- 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 (.not.visible) return 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 arg1_2 = ecol(1); arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) ! <============================ position the pen arg1_2 = ecol(2); arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! <============================ draw the trace line !Now, thicken the trace line by drawing 2 parallel lines, 1 pixel away if ((abs(ecol(1)-ecol(2)) > abs(erow(1)-erow(2)))) then ! trace line is roughly horizontal (W <--> E): arg1_2 = ecol(1); arg2_2 = erow(1) + 1 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) + 1 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1); arg2_2 = erow(1) - 1 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) - 1 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) else ! trace line is roughly vertical (N <--> S): arg1_2 = ecol(1) + 1; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) + 1; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1) - 1; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) - 1; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) end if !Consider whether it is necessary to outline fault trace with additional parallel traces (2 & 3 pixels away, on each side)? IF (LRDraw_checked) THEN !red_I4 = 255 ! red iRet = SETCOLORRGB(iFaultColor) ! from Global; usually red, unless user changed it. !---------------------------------------------------------------- if ((abs(ecol(1)-ecol(2)) > abs(erow(1)-erow(2)))) then ! trace line is roughly horizontal (W <--> E): arg1_2 = ecol(1); arg2_2 = erow(1) + 2 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) + 2 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1); arg2_2 = erow(1) - 2 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) - 2 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1); arg2_2 = erow(1) + 3 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) + 3 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1); arg2_2 = erow(1) - 3 ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2); arg2_2 = erow(2) - 3 ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) else ! trace line is roughly vertical (N <--> S): arg1_2 = ecol(1) + 2; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) + 2; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1) - 2; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) - 2; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1) + 3; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) + 3; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = ecol(1) - 3; arg2_2 = erow(1) ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = ecol(2) - 3; arg2_2 = erow(2) ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) end if !---------------------------------------------------------------- END IF ! LRDraw_checked; special two-color fault traces ! 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 arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = colB; arg2_2 = rowB ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) 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 arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) arg1_2 = colB; arg2_2 = rowB ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) xp = xp - size * cos(arg) yp = yp - size * sin(arg) call pixels(xp, yp, colA, rowA, visible) if (.not.visible) return arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) xp = xp - size * cos(normal) yp = yp - size * sin(normal) call pixels(xp, yp, colA, rowA, visible) if (.not.visible) return arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) 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 arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments call moveto(arg1_2, arg2_2, xy) xp = x + size * cos(normal) yp = y + size * sin(normal) call pixels(xp, yp, colB, rowB, visible) if (.not.visible) return arg1_2 = colB; arg2_2 = rowB ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) 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 arg1_2 = colC; arg2_2 = rowC ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) arg1_2 = colA; arg2_2 = rowA ! required INT(2) arguments status_2 = lineto(arg1_2, arg2_2) ! 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_2 = POLYGON($GFILLINTERIOR, poly, INT2(3)) end if end if end if end do ! next i 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 implicit none integer(4), intent(in) :: n1 real, intent(in) :: x, y integer(4), intent(out) :: n real R2min, Rmin, dA, dB, dG, R2 real tempvec(3) logical outside integer(4) i, iclose call xy2ABG(x,y, outside, tempvec) if (outside) 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 whether any continuum-element centroid exists at location? ! If so, return element number nE; otherwise return 0. ! SUBROUTINE ExistsCentroid(x, y, nE) ! Input: x, y ! Output: element# nE ! Requires: nodeABG, nodes, numEl, nomNod, tolerance USE Global IMPLICIT NONE REAL, INTENT(IN) :: x, y INTEGER(4), INTENT(OUT) :: nE REAL*4 R2min, Rmin, dA, dB, dG, R2 REAL*4 centroid(3), tempvec(3) LOGICAL outside INTEGER(4) i, iClose, n1, n2, n3 !- - - - - - - - - - - - - - - - - CALL xy2ABG(x,y, outside, tempvec) IF (outside) THEN nE = 0 ELSE R2min = 3.E+10 ! hoping that this will be replaced by a smaller number... iClose = 0 DO i = 1, numEl n1 = nodes(1, i) n2 = nodes(2, i) n3 = nodes(3, i) centroid(1) = (nodeABG(1, n1) + nodeABG(1, n2) + nodeABG(1, n3)) / 3.0 centroid(2) = (nodeABG(2, n1) + nodeABG(2, n2) + nodeABG(2, n3)) / 3.0 centroid(3) = (nodeABG(3, n1) + nodeABG(3, n2) + nodeABG(3, n3)) / 3.0 dA = tempvec(1) - centroid(1) dB = tempvec(2) - centroid(2) dG = tempvec(3) - centroid(3) R2 = dA*dA + dB*dB + dG*dG IF (R2 < R2min) THEN R2min = R2 iClose = i END IF END DO Rmin = sqrt(R2min) IF (Rmin <= tolerance) THEN nE = iClose ELSE nE = 0 END IF END IF ! outside, or not END SUBROUTINE ExistsCentroid ! ! Check whether any fault-element midpoint exists at location? ! If so, return fault number nF; otherwise return 0. ! SUBROUTINE ExistsMidpoint(x, y, nF) ! Input: x, y ! Output: element# nE ! Requires: nodeABG, nodes, numEl, nomNod, tolerance USE Global IMPLICIT NONE REAL, INTENT(IN) :: x, y INTEGER(4), INTENT(OUT) :: nF REAL*4 R2min, Rmin, dA, dB, dG, R2 REAL*4 midpoint(3), tempvec(3) LOGICAL outside INTEGER(4) i, iClose, n1, n2 !- - - - - - - - - - - - - - - - - CALL xy2ABG(x,y, outside, tempvec) IF (outside) THEN nF = 0 ELSE R2min = 3.E+10 ! hoping that this will be replaced by a smaller number... iClose = 0 DO i = 1, nFl n1 = nodeF(1, i) n2 = nodeF(2, i) midpoint(1) = (nodeABG(1, n1) + nodeABG(1, n2)) / 2.0 midpoint(2) = (nodeABG(2, n1) + nodeABG(2, n2)) / 2.0 midpoint(3) = (nodeABG(3, n1) + nodeABG(3, n2)) / 2.0 dA = tempvec(1) - midpoint(1) dB = tempvec(2) - midpoint(2) dG = tempvec(3) - midpoint(3) R2 = dA*dA + dB*dB + dG*dG IF (R2 < R2min) THEN R2min = R2 iClose = i END IF END DO Rmin = sqrt(R2min) IF (Rmin <= tolerance) THEN nF = iClose ELSE nF = 0 END IF END IF ! outside, or not END SUBROUTINE ExistsMidpoint ! Check whether element #n is upside-down (flipped)? subroutine Flipped(n) ! input: n (which element?) ! input: nodes, nodeABG (<--arrays in Global) ! output: EMemo (<--array in Global) use Global implicit none integer(4), intent(in) :: n real v12(3), v23(3), v31(3), vcross(3) integer(4) 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 implicit none real(4), intent(in) :: x, y integer(4), intent(out) :: jj real(4), intent(out) :: R2min integer(4) i integer(4) n1, n2, n3 logical visible real x1, y1, x2, y2, x3, y3 real xc, yc, dx, dy, R2 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) then tempv2(1) = nodeABG(1, n2) tempv2(2) = nodeABG(2, n2) tempv2(3) = nodeABG(3, n2) call ABG2xy(tempv2, visible, x2, y2) if (visible) then tempv3(1) = nodeABG(1, n3) tempv3(2) = nodeABG(2, n3) tempv3(3) = nodeABG(3, n3) call ABG2xy(tempv3, visible, x3, y3) if (visible) 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 implicit none 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 vertex# subroutine OnAnyE(n, n1, nlast, ie, je) !input: n, n1, nlast !output: ie (element), je(vertice) use global implicit none integer(4), intent(in) :: n, n1, nlast integer(4), intent(out) :: 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 5 nodes) ! subroutine Nearest_(x, y, numnear, nearones) ! input: x, y ! output: numNear, nearOnes(5) ! requires: numNod, nodeABG(numnod) use global implicit none real(4), intent(in) :: x, y integer(4), intent(out) :: numNear integer(4), intent(out) :: nearOnes(5) real tempvec(3) real R2toN, R2 real AN, BN, GN logical outside integer(4) i call xy2ABG(x, y, outside, tempvec) if (outside) then numnear = 0 nearones(1:5) = 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:5) = 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 == 5) 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 implicit none integer(4), intent(in) :: i, j integer(4), intent(out) :: Kfault, jf, Kele, je integer(4) n1, n2 integer(4) k integer(4) m1, m2, m3, m4 integer(4) Klow, Khigh 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 (vertex) use Global implicit none integer(4), intent(in) :: n, n1, nlast integer(4), intent(out) :: 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, visible(T/F) ! require: global scales(2,2) use global implicit none real(4), intent(in) :: x, y integer(4), intent(out) :: col, row logical, intent(out) :: visible real xcol, yrow ! note: real(8) require 1.D0 format 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 = .TRUE. col = xcol row = yrow else visible = .FALSE. 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) ! Changes color of one pixel, incrementally (according to predefined discrete color scale) ! Requires: Hirow, Hicol (<--in Global) use dfwin use dflib use global implicit none integer(4), intent(in) :: vcol(3), vrow(3) real x(4), y(4), r integer row1, row2, col1, col2 integer row, col, edge integer i, index, j, k integer(4) hue, istatus integer(2) col_2, row_2 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 col_2 = col; row_2 = row ! required INT(2) arguments hue = GETPIXELRGB(col_2, row_2) 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_2, row_2, 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; in radii) to (col,row; in pixels). ! ! Notes: The original version of this routine (in OrbWin 1.0, 1.1) ! assumed that monitor aspect ratio was 4:3 (e.g., 1024:768). ! This failed when wider monitors became common (e.g., 16:9, for HD movies). ! This revised version (for OrbWin 2.0) assumes that pixels are square. ! SUBROUTINE Scaler USE Global IMPLICIT NONE !(Most variables defined, and reside, in my MODULE Global.) REAL :: aspect_ratio ! (e.g., 4:3 = 1.33333, or 1920:1200 = 1.6, or 16:9 = 1.7778, .... aspect_ratio = REAL(HiCol) / REAL(HiRow) ! referring to the full-screen monitor size, which is the whole (virtual) size of the graphics window. !Note that R2C is the square of the distance (in pixels) from center of screen to corner of screen; ! thus, it is a test-value to see whether a point might *POSSIBLY* be visible? R2C = (0.5 * WindowHeight) ** 2 + (0.5 * WindowHeight * aspect_ratio) **2 scales(1, 1) = REAL(HiCol) / (WindowHeight * aspect_ratio) ! col = scales(1, 1) * x + scales(1, 2) scales(1, 2) = 0.5 * REAL(HiCol) scales(2, 1) = -REAL(HiRow) / WindowHeight ! row = scales(2, 1) * y + scales(2, 2) scales(2, 2) = 0.5 * REAL(HiRow) unscale(1, 1) = aspect_ratio * WindowHeight / REAL(HiCol) ! x = unscale(1, 1) * col + unscale(1, 2) unscale(1, 2) = -0.5 * aspect_ratio * WindowHeight unscale(2, 1) = -WindowHeight / REAL(HiRow) ! y = unscale(2, 1) * row + unscale(2, 2) unscale(2, 2) = 0.5 * WindowHeight !Old code follows (for historical reference): !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: Global's contour, iData ! Modify: oldBotF, oldIData, oldTopF ! Output: botF, dFC, logs, topF ! Requires: hiColor, eqcm, numNod use global implicit none 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.) 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) 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) ! because color index #1 = gray is not part of the spectrum 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) ! because color index #1 = gray is not part of the spectrum ! NeedToDraw omitted here oldidata = idata oldbotf = botf oldtopf = topf end subroutine SetBins ! ! Figure out best assignment of (hiColor -1) colors to element LR# values ! when iData = 0, indicating that data are taken from continuum_LRi(1:mxel); ! HOWEVER, also consider value range in fault_LRi, as affecting the total ! range of assigned colors. ! subroutine SetBinsLR ! input: contour, iData (from Global) ! modify: oldBotF, oldIData, oldTopF ! output: botF, dFC, logs, topF ! requires: hiColor, eqcm, numEl use global implicit none logical got_one, visible real maxF, minF, x, y integer(4) col, i, LRi, n, n1, n2, n3, row REAL*4 tempVec(3) ! Determine range of data VISIBLE IN THE CURRENT WINDOW, ! and use that range as the basis for dividing the color bar. got_one = .FALSE. DO i = 1, numEl LRi = continuum_LRi(i) n1 = nodes(1, i) n2 = nodes(2, i) n3 = nodes(3, i) tempVec(1) = (nodeABG(1, n1) + nodeABG(1, n2) + nodeABG(1, n3)) / 3.0 tempVec(2) = (nodeABG(2, n1) + nodeABG(2, n2) + nodeABG(2, n3)) / 3.0 tempVec(3) = (nodeABG(3, n1) + nodeABG(3, n2) + nodeABG(3, n3)) / 3.0 call ABG2xy (tempVec, visible, x, y) if (visible) 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, 1.0 * continuum_LRi(i)) minf = MIN(minf, 1.0 * continuum_LRi(i)) ELSE ! this IS the first node found in the window maxf = 1.0 * continuum_LRi(i) 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 elements DO i = 1, nFl LRi = fault_LRi(i) n1 = nodeF(1, i) n2 = nodeF(2, i) tempVec(1) = (nodeABG(1, n1) + nodeABG(1, n2)) / 2.0 tempVec(2) = (nodeABG(2, n1) + nodeABG(2, n2)) / 2.0 tempVec(3) = (nodeABG(3, n1) + nodeABG(3, n2)) / 2.0 call ABG2xy (tempVec, visible, x, y) if (visible) 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, 1.0 * fault_LRi(i)) minf = MIN(minf, 1.0 * fault_LRi(i)) ELSE ! this IS the first node found in the window maxf = 1.0 * fault_LRi(i) 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 faults IF (.NOT.got_one) THEN ! empty window (no elements visible in window?) maxF = 0.0 minF = 0.0 END IF ! Automatically let colors span the data range, ! using linear scale, and thus allowing values of zero: logs = .FALSE. botF = minF topF = maxF !Note: It is not a problem if topf = botf = 0.0 dFC = (topF - botF) / REAL(hiColor - 1) ! because color index #1 = gray is not part of the spectrum oldIData = iData oldBotF = botF oldTopF = topF end subroutine SetBinsLR ! ! shading subroutine used in DrawElement ! Note: using Polygon fill to speed things up! ! SUBROUTINE Shade(nPnts, eCol, eRow, cList, rList, index) ! input: nPnts, eCol(3), eRow(3), cList(10), rList(10), index USE DFLib USE Global IMPLICIT NONE INTEGER(4), INTENT(IN) :: nPnts INTEGER(4), INTENT(IN) :: eCol(3), eRow(3) INTEGER(4), INTENT(IN) :: cList(10), rList(10) INTEGER(4), INTENT(IN) :: index INTEGER(4) temp_index INTEGER(4) iRet INTEGER(4) jt, jp INTEGER(4) c1, c2, rw1, rw2 INTEGER(4) hue, iColor INTEGER(2) arg1_2, arg2_2, status_2 TYPE(xyCoord) xy TYPE(xyCoord), ALLOCATABLE:: poly(:) !- - - - - - - - - - - - - - - - - - - - - - - ! If color out-of-range, set to lowest or highest color: temp_index = MAX(1, MIN(hiColor, index)) !N.B. Color 1 (gray) is used for element-center icons of deleted elements, even though not part of normal spectrum. colorPick = colorArray(temp_index) iRet = SetColorRGB(colorPick) ALLOCATE ( poly(nPnts) ) DO jt = 1, nPnts IF (jt < nPnts) THEN jp = jt + 1 ELSE jp = 1 END IF c1 = cList(jt) c2 = cList(jp) rw1 = rList(jt) rw2 = rList(jp) arg1_2 = c1; arg2_2 = rw1 ! required INT(2) arguments CALL MoveTo (arg1_2, arg2_2, xy) arg1_2 = c2; arg2_2 = rw2 ! required INT(2) arguments status_2 = LineTo(arg1_2, arg2_2) poly(jt)%xCoord = c1 poly(jt)%yCoord = rw1 END DO status_2 = POLYGON($GFILLINTERIOR, poly, int2(nPnts)) DEALLOCATE (poly) END SUBROUTINE Shade ! ! Calculate solid angle of spherical triangle in steradians ! 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 Steradians(v1, v2, v3, solidAngle) USE DSphere ! provided by Peter Bird as file DSphere.f90 ! Specifically, CALL DSpherical_Area() replaces old code here ! which unfortunately produced "NaN" results in some cases. IMPLICIT NONE REAL(4), INTENT(IN) :: v1(3), v2(3), v3(3) REAL(4), INTENT(OUT) :: solidAngle REAL(8) :: a, across, adot, answer, aside, b, ba(3), bcross, bdot, bside, & & c, cb(3), ccross, cdot, center(3), coshalf, cside, & & denominator, dot, half, la, lb, lc, & & ra(3), rb(3), rc(3), s, sinhalf, vcross(3) REAL(4) :: ATan2F ra(1:3) = DBLE(v1(1:3)) rb(1:3) = DBLE(v2(1:3)) rc(1:3) = DBLE(v3(1:3)) ! Ensure they are still truly unit vectors after conversion to REAL(8): 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 !------------------------------------------------------- CALL DSpherical_Area(ra, rb, rc, answer) ! in MODULE DSphere !N.B. answer will be negative if ra-->rb-->rc-->ra is a CLOCKWISE circuit. solidAngle = answer ! REAL(4) = REAL(8), for export !------------------------------------------------------- ! N.B. The reason for the existence of this routine is just to provide a ! REAL(4) interface for DSpherical_Area which is native REAL(8). END SUBROUTINE Steradians ! convert pixel coordinates to (x,y) using predefined transformation subroutine XandY (col, row, x,y) ! input: col, row ! output: x,y ! requires: Global array unscale(2,2) use global implicit none integer(4), intent(in) :: col, row real(4), intent(out) :: 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 implicit none real(4), intent(in) :: x, y logical, intent(out) :: outside real(4), intent(out) :: tempv(3) real r2, outofplane r2 = x**2 + y**2 if (r2 >= 1.0) then outside = .TRUE. tempv(1) = 0. tempv(2) = 0. tempv(3) = 0. else outside = .FALSE. 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 implicit none real(4), intent(in) :: x, y logical, intent(out) :: outside real(4), intent(out) :: lon, lat real r2, outofplane real tempvec(3), tempv2(3) r2= x*x + y*y if (r2>=1) then outside = .TRUE. lon = 0. lat = 0. else outside = .FALSE. 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 implicit none real(8), intent(inout) :: 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) implicit none real(4), intent(inout) :: tempvec(3) real r2, r 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 ! requires: subroutine Cross implicit none real(4), intent(in) :: winlat, winlon real(4), intent(out) :: 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