C$Procedure MKMASK ( Set plate mutual visibility flag ) SUBROUTINE MKMASK ( NVOX, . VOXSIZ, VOXORI, . NVXTOT, VXPTR, . NVXLST, VXLIST, . EXTENT, . LSIDE, . NP, VERT, NORM, CENT, . V1, V2 , V3, . PMASK, MAXELE ) C$ Abstract C C For each plate in the model, this routine determines if other C plates can be seen from it. C C$ Copyright C C Copyright (1999), California Institute of Technology. C U.S. Government sponsorship acknowledged. C C$ Required_Reading C C None. C C$ Keywords C C plate mask C C$ Declarations IMPLICIT NONE INCLUDE 'pltmax.inc' INTEGER NVOX ( 3 ) DOUBLE PRECISION VOXSIZ DOUBLE PRECISION VOXORI ( 3 ) INTEGER NVXTOT INTEGER VXPTR ( * ) INTEGER NVXLST INTEGER VXLIST ( * ) DOUBLE PRECISION EXTENT ( 6 ) DOUBLE PRECISION LSIDE ( * ) INTEGER NP DOUBLE PRECISION VERT ( 3, 3, * ) DOUBLE PRECISION NORM ( 3, * ) DOUBLE PRECISION CENT ( 3, * ) INTEGER V1 ( * ) INTEGER V2 ( * ) INTEGER V3 ( * ) INTEGER PMASK ( * ) DOUBLE PRECISION MAXELE ( * ) C$ Brief_I/O C C VARIABLE I/O DESCRIPTION C -------- --- -------------------------------------------------- C NP I Number of plates. C VERT I Array containing plate vertex coordinates. C NORM I Array containing plate outward normal vectors. C CENT I Array containing plate center coordinates. C V1 I Array of plate vertex IDs (first vertex). C V2 I Array of plate vertex IDs (second vertex). C V3 I Array of plate vertex IDs (third vertex). C MAXPLT I Declared size of array/matrix variables. C PMASK O Array of mask flag values. C MXELEV O Array of maximum occluded elevations. C C$ Detailed_Input C C NP Number of plates to be checked for mutual visibility. C C VERT Array of vertex coordinates of each plate's vertices. C C NORM Array of double precision 3-vectors containing each C plate's outward-pointing normal vector. C C CENT Array of double precision 3-vectors containing the C location of each plate's center. C C V1 C C V2 C C V3 C C$ Detailed_Output C C PMASK Mask flag values: 0 if no plate is visible and 1 C if it is. C C MAXELE Maximum elevation of any obstruction visible from C plate-A. C C$ Parameters C C None. C C$ Files C C None. C C$ Exceptions C C Error free. C C$ Particulars C C Later. C C$ Examples C C Later. C C$ Restrictions C C None. C C$ Author_and_Institution C C J.A. Bytof (JPL) C C$ Literature_References C C None. C C$ Version C C- SPICELIB Version 2.1.0, XX-XXX-2004 (EDW) C C Under repair: C C Added error check on array sizes. C C Removed TRIPUP call, ANORM, ACENT, OVERID, C VISLST variables. C Restuctured logic to remove uneeded variables and C loops. C C- SPICELIB Version 2.0.2, 26-AUG-2002 (BVS) C C Removed WRITE. C C- SPICELIB Version 2.0.1, 27-JUN-2002 (JAB) C C Added a write statement, so that the user can C following the progress of horizon mask processing. C C- SPICELIB Version 2.0.0, 03-JAN-1999 (JAB) C C Made this routine independent of the plate manager, C by passing in vertex ID arrays instead of calling C a plate manager routine. C C- SPICELIB Version 1.0.0, 28-JAN-1999 (JAB) C C C-& C$ Index_Entries C C plate mask C C-& DOUBLE PRECISION HALFPI DOUBLE PRECISION VDOT DOUBLE PRECISION VNORM C C Local Variables. C LOGICAL RETURN LOGICAL INPRZM LOGICAL ABOVE ( 3 ) INTEGER I INTEGER IA INTEGER IB INTEGER NVIS INTEGER AVID ( 3 ) INTEGER BVID ( 3 ) DOUBLE PRECISION OVERMX DOUBLE PRECISION TMPMAX DOUBLE PRECISION AVERT ( 3, 3 ) DOUBLE PRECISION BVERT ( 3, 3 ) LOGICAL INLIST INTEGER PINHFS INTEGER MAXVOX PARAMETER ( MAXVOX = 600000 ) DOUBLE PRECISION ANGRAD DOUBLE PRECISION APLSUM DOUBLE PRECISION BASE DOUBLE PRECISION CTRCON DOUBLE PRECISION CTRELV DOUBLE PRECISION HYP DOUBLE PRECISION KCELV DOUBLE PRECISION KCTOFF ( 3 ) DOUBLE PRECISION KHGHT DOUBLE PRECISION KPERP ( 3 ) DOUBLE PRECISION KPRPLN DOUBLE PRECISION KVERT ( 3 ) DOUBLE PRECISION MAXCTR DOUBLE PRECISION MAXELV ( MAXPLT ) DOUBLE PRECISION PLCON ( MAXPLT ) DOUBLE PRECISION SPHRAD DOUBLE PRECISION VOXCTR ( 3, MAXVOX ) DOUBLE PRECISION VOXRAD INTEGER APPSUM INTEGER EXTPYR ( MAXPLT ) INTEGER J INTEGER K INTEGER NNE INTEGER PLTSET ( MAXPLT ) INTEGER VOXLOC ( 3 ) INTEGER VOXNPL ( MAXVOX ) INTEGER NABOVE INTEGER NBELOW LOGICAL CHKTRI C C Standard SPICE error handling. C IF ( RETURN() ) THEN RETURN END IF CALL CHKIN ( 'MKMASK' ) IF ( NP .GT. MAXPLT ) THEN CALL SETMSG ( 'Index larger than array. '// . 'Index = #1. Array size = #2.') CALL ERRINT ( '#1', NP ) CALL ERRINT ( '#2', MAXPLT ) CALL SIGERR ( 'SPICE(ARRAYTOOSMALL)' ) CALL CHKOUT ( 'MKMASK' ) RETURN END IF C C Check the number of voxels in the spatial index. C IF ( NVXTOT .GT. MAXVOX ) THEN CALL SETMSG ( 'Spatial index too large. ' // . 'Size = #1. Max allowed array size = #2.' ) CALL ERRINT ( '#1', NVXTOT ) CALL ERRINT ( '#2', MAXVOX ) CALL SIGERR ( 'SPICE(ARRAYTOOSMALL)' ) CALL CHKOUT ( 'MKMASK' ) RETURN END IF C C Compute the the center and radius of each voxel. The radius C is the distance from the center to a corner of the voxel. C NNE = 0 DO I = 1, NVXTOT C C Get the Ith voxel's grid coordinates. C CALL ID2VOX ( I, NVOX, VOXLOC ) DO J = 1, 3 VOXCTR(J,I) = VOXORI(J) + (VOXLOC(J) - 0.5)*VOXSIZ END DO C C Save the voxel's plate count. C VOXNPL(I) = VXLIST ( VXPTR(I) ) IF ( VOXNPL(I) .GT. 0 ) THEN NNE = NNE + 1 END IF END DO C C Initialize PLTSET. C CALL CLEARI ( NP, PLTSET ) C C Initialize EXTPYR. C CALL CLEARI ( NP, EXTPYR ) VOXRAD = (SQRT(3.D0)/2.D0) * VOXSIZ C C Compute the plane constant for each plate. C DO IA = 1, NP PLCON(IA) = VDOT ( NORM(1,IA), VERT(1,1,IA) ) END DO C C Outer loop over all plates. C APPSUM = 0 APLSUM = 0 DO IA = 1, NP C C Default mask flag value. C PMASK(IA) = 0 C C Get plate-A's vertex IDs. C AVID(1) = V1( IA ) AVID(2) = V2( IA ) AVID(3) = V3( IA ) C C Determine which voxels contain applicable plates. C For each voxel, test whether the voxel is entirely C in the half space below the current plate. If so, C the voxel is not applicable. C C Get the unit normal and plane constant for this plate. C C C Initialize the max center elevation for this plate. C MAXCTR = 0.D0 NABOVE = 0 NBELOW = 0 DO I = 1, NVXTOT IF ( VOXNPL(I) .GT. 0 ) THEN C C The voxel is non-empty. C C Compute the plane constant for the voxel's center. C If the constant is less than PLCON-VOXRAD, the entire C voxel is below the plate. Otherwise, the voxel C is "applicable." C CTRCON = VDOT( NORM(1,IA), VOXCTR(1,I) ) IF ( CTRCON .GE. ( PLCON(IA) - VOXRAD ) ) THEN C C The voxel is applicable. C C Traverse the plate list for the current voxel. C Flip the applicability switch for each plate. C If we can prove that a plate is in the exterior C of plate IA's pyramid, flip the pyramid exterior C switches for that plate. Also find an upper C elevation bound for each plate. C DO J = 1, VOXNPL(I) K = VXLIST( VXPTR(I) + J ) C C It's possible that this plate is on the plate list C for a voxel that has already been processed. If C so, PLTSET(K) will already be set to +/-IA. If C PLTSET(K) is -IA, that means we've already found C that plate K cannot contribute to the elevation C mask. In this case, we skip over plate K. C IF ( PLTSET(K) .NE. -IA ) THEN C C If we already know that plate K is outside of C plate IA's pyramid, we don't need to do this C computation. However, if plate K belongs to a C voxel we've already processed, but we haven't C found that plate K is outside the pyramid, it C may pay off to re-do the computation using the C current voxel. C IF ( PLTSET(K) .NE. IA ) THEN C C We haven't processed plate K before. C C Determine cheaply computed bounds on the C elevation of plate K as seen from plate IA. C For these computations, we'll refer to the C orthogonal projection onto the plane C containing plate IA of the center of plate C K. Call this vector KCPROJ. We'll also the C need the the offset of the center of plate K C from the center of plate IA; call this C KCTOFF. C CALL VSUB ( CENT(1,K), CENT(1,IA), KCTOFF ) C C Find the component of KCTOFF perpendicular C to the normal to plate IA. Let KPRPLN be C the length of this component. C CALL VPERP ( KCTOFF, NORM(1,IA), KPERP ) KPRPLN = VNORM ( KPERP ) C C Find the height of the center of plate K C above the plane. C CALL VSUB ( KCTOFF, KPERP, KVERT ) IF ( VDOT( NORM(1,IA), KVERT ) . .GT. 0.D0 ) THEN KHGHT = VNORM(KVERT) ELSE KHGHT = - VNORM(KVERT) END IF C C Find the elevation of the center of plate K C relative to the center of plate IA. This C elevation is a lower bound (not necessarily C a good one) for the elevation mask value for C plate IA. C CTRELV = ATAN2 ( KHGHT, KPRPLN ) C C Update the maximum plate center elevation. C The elevation mask for plate IA is at least C as great as this value. C MAXCTR = MAX ( CTRELV, MAXCTR ) C C Now find a reasonably sharp upper bound on C the elevation of plate K. C C Plate K is enclosed in a sphere centered at C plate K's center and having radius LSIDE(K). C All points of plate IA are contained in a C circle centered at the "center" of plate IA, C having radius LSIDE(IA), and lying in the C plane containing plate IA. So the maximum C elevation of any point in the sphere C enclosing plate K, as seen from the closest C point on plate IA's bounding circle to C KCPROJ, is an upper bound for the elevation C of plate K with respect to plate IA. C C Set the radius of the bounding sphere. C SPHRAD = LSIDE(K) C C Let VIEWPT be a point located on KPERP at a C distance of LSIDE(IA) away from plate IA's C center. Let BASE be the distance from C VIEWPT to KCPROJ. C BASE = VNORM(KPERP) - LSIDE(IA) C C If the sphere bounding plate K is too large, C that is, if its projection onto the plane C containing plate IA overlaps the bounding C circle around plate IA, we have a degenerate C geometric case. Set the upper elevation C limit to pi/2. C IF ( BASE .LE. SPHRAD ) THEN MAXELV( K ) = HALFPI() ELSE C C The bounding sphere around plate K is not C "over" the bounding circle around plate C IA. We can conclude that plate K is in C the exterior of the pyramid above plate C IA. Set EXTPYR(K) to indicate this. C EXTPYR(K) = IA C C Find the angular radius of plate K's C bounding sphere as seen VIEWPT. The C distance from VIEWPT to the center of the C sphere is the hypoteneuse of a right C triangle whose base is BASE and whose C other leg has length |KHGHT|. C HYP = SQRT ( BASE*BASE + KHGHT*KHGHT ) ANGRAD = ASIN ( SPHRAD / HYP ) C WRITE (*,*) 'ANGRAD = ', DPR()*ANGRAD C C The maximum elevation of the bounding C sphere is the sum of the elevation of the C center relative to VIEWPT and of the C angular radius. Our construction should C prevent this sum from exceeding 2*pi, but C we'll take no chances. C KCELV = ATAN2 ( KHGHT, BASE ) C WRITE (*,*) 'KCELV = ', DPR()*KCELV MAXELV(K) = MIN ( KCELV+ANGRAD, HALFPI() ) C C MAXELV(K) is set for the non-degenerate C case. C END IF END IF C C Now it's time to decide what to do about plate C K. If the upper elevation bound for plate K C is below MAXCTR, then plate K plays no role C in determining the horizon mask for plate IA. C We need not consider plate K any longer in this C iteration of the outer loop. Signify this C situation by setting PLTSET to -IA. C C Note that while this step cuts down on the C number of plates that can be re-processed during C this loop, MAXCTR tends to grow as the loop C progresses, so it's possible that some plates C that can be excluded are not excluded during C this step. We should again check MAXELV(K) C against MAXCTR after this loop completes. C IF ( MAXELV(K) .LT. MAXCTR ) THEN PLTSET(K) = -IA ELSE C C Indicate that plate K is applicable to plate C IA. C C WRITE (*,*) 'turning on ', K PLTSET(K) = IA END IF END IF C C We've processed plate K. C END DO C C We've processed the Jth plate in the list for voxel I. C END IF C C We've processed voxel I if this voxel is applicable. C END IF C C We've processed voxel I if this voxel is non-empty. C END DO C C We've processed all voxels in the spatial index. C C At this point, C C PLTSET(K) is set to +/- IA for every plate K in the set of C applicable plates. C C EXTPYR(K) is set to IA for every plate known to be outside C of the triangular pyramid above plate IA, for every K in C the set of applicable plates. C C MAXELV(K) is set to an upper elevation bound for every C plate C K in the set of applicable plates. C C Always turn off plate IA in the applicability list. C PLTSET(IA) = 0 TMPMAX = -99999.D0 OVERMX = -99999.D0 CALL CHKINI () DO IB = 1, NP C C Consider only plates that are on the applicability C list for plate IA. C IF ( ( PLTSET(IB) .EQ. IA ) . .AND. ( MAXELV(IB) .GT. MAXCTR ) ) THEN C C Get plate-B's vertex IDs. C BVID(1) = V1( IB ) BVID(2) = V2( IB ) BVID(3) = V3( IB ) C C Give some vectors shorter names. C CALL VEQU ( VERT(1,1,IA), AVERT(1,1) ) CALL VEQU ( VERT(1,2,IA), AVERT(1,2) ) CALL VEQU ( VERT(1,3,IA), AVERT(1,3) ) CALL VEQU ( VERT(1,1,IB), BVERT(1,1) ) CALL VEQU ( VERT(1,2,IB), BVERT(1,2) ) CALL VEQU ( VERT(1,3,IB), BVERT(1,3) ) C C Are any vertices of plate-B above plate-A's C "local horizon"? If yes, return a list of them. C If plate-A is connected to plate-B, common vertices C are not listed. C NVIS = 0 ABOVE(1) = .FALSE. ABOVE(2) = .FALSE. ABOVE(3) = .FALSE. DO I = 1, 3 IF ( .NOT. INLIST( 3, AVID, BVID(I)) ) THEN IF ( PINHFS( BVERT(1,I), . CENT(1,IA), . NORM(1,IA) ) .EQ. 1 ) THEN NVIS = NVIS + 1 ABOVE(I) = .TRUE. END IF END IF END DO IF ( NVIS .GT. 0 ) THEN PMASK(IA) = 1 C C Does any part of plate-B lie within the right C triangular prism having plate-A at one end? C If the answer is yes, set the maximum elevation C for plate-A at 90 degrees. If no, determine the C maximum elevation by checking every case. C CHKTRI = .TRUE. IF ( EXTPYR(IB) .NE. IA ) THEN C C Plate IB might intersect the pyramid. Check. C IF ( INPRZM ( AVERT , BVERT, . NORM(1,IA), CENT(1,IA) ) ) THEN TMPMAX = 90.D0 CHKTRI = .FALSE. END IF ELSE END IF IF ( CHKTRI ) THEN C C We may need to find the elevation of plate IB C relative to plate IA. However, if our elevation C bound for plate IB is less than the current C elevation mask value for plate IA, plate IB can't C influence the mask value and thus needn't be C considered. C C Note that OVERMX is stored in units of degrees. C convert MAXELV to degrees for comparison. C IF ( MAXELV(IB) .GE. MAXCTR ) THEN C C We do need to find the elevation of plate IB C relative to plate IA. C CALL TRAZEL ( AVERT, NORM(1,IA), BVERT, . BVID, ABOVE, TMPMAX ) C C Store the three processed vertices. C CALL CHKREP( V1( IB ), V2( IB ), V3( IB ) ) ELSE C C We should never get here. C CALL SIGERR ( 'BUG' ) END IF END IF C C We've checked the elevation of plate IB relative to C plate IA. C ELSE TMPMAX = 0.D0 END IF OVERMX = MAX ( OVERMX, TMPMAX ) END IF C C End of B-plate loop C C IB = IB + 1 END DO C C Assign the max elevation element. This value C is greater-than C MAXELE(IA) = OVERMX IF ( (IA-(IA/1000)*1000) .EQ. 0 ) THEN WRITE (*,*) ' processed ', IA, ' of ', NP,' plates' END IF END DO CALL CHKOUT ( 'MKMASK' ) RETURN END