diff options
Diffstat (limited to 'sys/gio/ncarutil/isosrf.f')
| -rw-r--r-- | sys/gio/ncarutil/isosrf.f | 1696 |
1 files changed, 1696 insertions, 0 deletions
diff --git a/sys/gio/ncarutil/isosrf.f b/sys/gio/ncarutil/isosrf.f new file mode 100644 index 00000000..7be532ee --- /dev/null +++ b/sys/gio/ncarutil/isosrf.f @@ -0,0 +1,1696 @@ + SUBROUTINE ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,SLAB,TISO,IFLAG) +C +C +-----------------------------------------------------------------+ +C | | +C | Copyright (C) 1986 by UCAR | +C | University Corporation for Atmospheric Research | +C | All Rights Reserved | +C | | +C | NCARGRAPHICS Version 1.00 | +C | | +C +-----------------------------------------------------------------+ +C +C +C +C +C DIMENSION OF T(LU,LV,MW),EYE(3),SLAB(MUVWP2,MUVWP2) +C ARGUMENTS +C +C LATEST REVISION DECEMBER 1984 +C +C PURPOSE ISOSRF DRAWS AN APPROXIMATION OF AN ISO-VALUED +C SURFACE FROM A THREE-DIMENSIONAL ARRAY WITH +C HIDDEN LINES REMOVED. +C +C USAGE IF THE FOLLOWING ASSUMPTIONS ARE MET, USE +C +C CALL EZISOS (T,MU,MV,MW,EYE,SLAB,TISO) +C +C ASSUMPTIONS: +C -- ALL OF THE T ARRAY IS TO BE USED. +C -- IFLAG IS CHOSEN INTERNALLY. +C -- FRAME IS CALLED BY EZISOS. +C +C IF THE ASSUMPTIONS ARE NOT MET, USE +C +C CALL ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2, +C SLAB,TISO,IFLAG) +C +C ARGUMENTS +C +C ON INPUT T +C THREE DIMENSIONAL ARRAY OF DATA THAT DEFINES +C THE ISO-VALUED SURFACE. +C +C LU +C FIRST DIMENSION OF T IN THE CALLING PROGRAM. +C +C MU +C THE NUMBER OF DATA VALUES OF T TO BE +C PROCESSED IN THE U DIRECTION (THE FIRST +C SUBSCRIPT DIRECTION). WHEN PROCESSING THE +C ENTIRE ARRAY, LU = MU (AND LV = MV). +C +C LV +C SECOND DIMENSION OF T IN THE CALLING PROGRAM. +C +C MV +C THE NUMBER OF DATA VALUES OF T TO BE +C PROCESSED IN THE V DIRECTION (THE SECOND +C SUBSCRIPT DIRECTION). +C +C MV +C THE NUMBER OF DATA VALUES OF T TO BE +C PROCESSED IN THE W DIRECTION (THE THIRD +C SUBSCRIPT DIRECTION). +C +C EYE +C THE POSITION OF THE EYE IN THREE-SPACE. T IS +C CONSIDERED TO BE IN A BOX WITH OPPOSITE +C CORNERS (1,1,1) AND (MU,MV,MW). THE EYE IS +C AT (EYE(1),EYE(2),EYE(3)), WHICH MUST BE +C OUTSIDE THE BOX THAT CONTAINS T. WHILE GAINING +C EXPERIENCE WITH THE ROUTINE, A GOOD CHOICE +C FOR EYE MIGHT BE (5.0*MU,3.5*MV,2.0*MW). +C +C MUVWP2 +C THE MAXIMUM OF (MU,MV,MW)+2; THAT IS, +C MUVWP2 = MAX(MU,MV,MW)+2). +C +C SLAB +C A WORK SPACE USED FOR INTERNAL STORAGE. SLAB +C MUST BE AT LEAST MUVWP2*MUVWP2 WORDS LONG. +C +C TISO +C THE ISO-VALUE USED TO DEFINE THE SURFACE. THE +C SURFACE DRAWN WILL SEPARATE VOLUMES OF T THAT +C HAVE VALUES GREATER THAN OR EQUAL TO TISO FROM +C VOLUMES OF T THAT HAVE VALUES LESS THAN TISO. +C +C IFLAG +C THIS FLAG SERVES TWO PURPOSES. +C . FIRST, THE ABSOLUTE VALUE OF IFLAG +C DETERMINES WHICH TYPES OF LINES ARE DRAWN +C TO APPROXIMATE THE SURFACE. THREE TYPES +C OF LINES ARE CONSIDERED: LINES OF +C CONSTANT U, LINES OF CONSTANT V AND LINES +C OF CONSTANT W. THE FOLLOWING TABLE LISTS +C THE TYPES OF LINES DRAWN. +C +C LINES OF CONSTANT +C ----------------- +C IABS(IFLAG) U V W +C ----------- --- --- --- +C 1 NO NO YES +C 2 NO YES NO +C 3 NO YES YES +C 4 YES NO NO +C 5 YES NO YES +C 6 YES YES NO +C 0, 7 OR MORE YES YES YES +C +C . SECOND, THE SIGN OF IFLAG DETERMINES WHAT +C IS INSIDE AND WHAT IS OUTSIDE, HENCE, +C WHICH LINES ARE VISIBLE AND WHAT IS DONE +C AT THE BOUNDARY OF T. FOR IFLAG: +C +C POSITIVE T VALUES GREATER THAN TISO ARE +C ASSUMED TO BE INSIDE THE SOLID +C FORMED BY THE DRAWN SURFACE. +C NEGATIVE T VALUES LESS THAN TISO ARE +C ASSUMED TO BE INSIDE THE SOLID +C FORMED BY THE DRAWN SURFACE. +C IF THE ALGORITHM DRAWS A CUBE, REVERSE THE +C SIGN OF IFLAG. +C +C ON OUTPUT T,LU,MU,LV,MV,MW,EYE,MUVWP2,TISO AND IFLAG ARE +C UNCHANGED. SLAB HAS BEEN WRITTEN IN. +C +C NOTE . THIS ROUTINE IS FOR LOWER RESOLUTION ARRAYS +C THAN ISOSRFHR. 40 BY 40 BY 40 IS A +C PRACTICAL MAXIMUM. +C . TRANSFORMATIONS CAN BE ACHIEVED BY +C ADJUSTING SCALING STATEMENT FUNCTIONS IN +C ISOSRF, SET3D AND TR32. +C . THE HIDDEN-LINE ALGORITHM IS NOT EXACT, SO +C VISIBILITY ERRORS CAN OCCUR. +C . THREE-DIMENSIONAL PERSPECTIVE CHARACTER +C LABELING OF ISOSRF IS POSSIBLE BY USING +C THE UTILITY PWRZI. FOR A DESCRIPTION OF +C THE USAGE, SEE THE PWRZI DOCUMENTATION. +C +C ENTRY POINTS ISOSRF, EZISOS, SET3D, TRN32I, ZEROSC, +C STCNTR, DRCNTR, TR32, FRSTS, KURV1S, KURV2S, +C FRSTC, FILLIN, DRAWI, ISOSRB, MMASK +C +C COMMON BLOCKS ISOSR1, ISOSR2, ISOSR3, ISOSR4, ISOSR5, +C ISOSR6, ISOSR7, ISOSR8, ISOSR9, TEMPR, +C PWRZ1I +C +C REQUIRED LIBRARY THE ERPRT77 PACKAGE AND THE SPPS. +C ROUTINES +C +C I/O PLOTS SURFACE +C +C PRECISION SINGLE +C +C LANGUAGE FORTRAN 77 +C +C HISTORY DEVELOPED FOR USERS OF ISOSRFHR WITH SMALLER +C ARRAYS. +C +C ALGORITHM CUTS THROUGH THE THREE-DIMENSIONAL ARRAY ARE +C CONTOURED WITH A SMOOTHING CONTOURER WHICH ALSO +C MARKS A MODEL OF THE PLOTTING PLANE. INTERIORS +C OF BOUNDARIES ARE FILLED IN AND THE RESULT IS +C .OR.ED INTO ANOTHER MODEL OF THE PLOTTING PLANE +C WHICH IS USED TO TEST SUBSEQUENT CONTOUR LINES +C FOR VISIBILITY. +C +C TIMING VARIES WIDELY WITH SIZE OF T AND THE VOLUME OF +C THE SPACE ENCLOSED BY THE SURFACE DRAWN. +C +C **NOTE** SPACE REQUIREMENTS CAN BE REDUCED BY +C CHANGING THE SIZE OF THE ARRAYS ISCR, ISCA +C (FOUND IN COMMON ISOSR2), MASK(FOUND IN +C COMMON ISOSR5) AND THE VARIABLE NBPW +C (COMMON ISOSR5). +C ISCR AND ISCA NEED 128X128 BITS. SO ON A +C 64 BIT MACHINE ISCR, ISCA CAN BE +C DIMENSIONED TO (2,128). NBPW SET IN +C SUBROUTINE MMASK SHOULD CONTAIN THE +C NUMBER OF BITS PER WORD YOU WISH TO +C UTILIZE. +C THE DIMENSION OF MASK AND NMASK SHOULD +C EQUAL THE VALUE OF NBPW. +C LS SHOULD BE SET TO THE FIRST DIMENSION +C OF ISCA AND ISCR. +C +C EXAMPLES: +C ON A 60 BIT MACHINE: +C DIMENSION ISCA(4,128), ISCR(4,128) +C DIMENSION MASK(32) +C NBPW = 32 +C ON A 64 BIT MACHINE: +C DIMENSION ISCA(2,128), ISCR(2,128) +C DIMENSION MASK(64) +C NBPW = 64 +C +C INTERNAL PARAMETERS NAME DEFAULT FUNCTION +C ---- ------- -------- +C IREF 1 FLAG TO CONTROL DRAWING OF AXES. +C .IREF=NONZERO DRAW AXES. +C .IREF=ZERO DO NOT DRAW AXES. +C +C + SAVE + DIMENSION T(LU,LV,MW),EYE(3) ,SLAB(MUVWP2,MUVWP2) +C + COMMON /ISOSR1/ ISLBT ,U ,V ,W + COMMON /ISOSR2/ LX ,NX ,NY ,ISCR(8,128), + 1 ISCA(8,128) + COMMON /ISOSR3/ ISCALE ,XMIN ,XMAX ,YMIN , + 1 YMAX ,BIGD ,R0 + COMMON /ISOSR4/ RX ,RY + COMMON /ISOSR5/ NBPW ,MASK(16) ,GENDON + COMMON /ISOSR6/ IX ,IY ,IDX ,IDY , + 1 IS ,ISS ,NP ,CV , + 2 INX(8) ,INY(8) ,IR(500) ,NR + COMMON /ISOSR7/ IENTRY ,IONES + COMMON /ISOSR9/ BIG ,IXBIT +C + LOGICAL GENDON + DATA IREF/1/ +C + AVE(A,B) = (A+B)*.5 +C +C A.S.F. FOR SCALING +C + SU(UTEMP) = UTEMP + SV(VTEMP) = VTEMP + SW(WTEMP) = WTEMP +C +C +NOAO - Blockdata ISOSRB rewritten as run time initialization +C EXTERNAL ISOSRB + call isosrb +C -NOAO +C +C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR +C + CALL Q8QST4 ('NSSL','ISOSRF','ISOSRF','VERSION 12') + NERR = 0 +C +C 3-SPACE U,V,W,IU,IV,IW,ETC +C 2-SPACE X,Y,IX,IY,ETC +C +C INITIALIZE MASKS +C + IF (.NOT.GENDON) CALL MMASK +C +C SET SHIFT VALUE FOR X,Y PACKING +C +C IF YOUR MACHINE HAS MORE THAN 16 BITS PER WORD THIS CHECK MAY BE +C MODIFIED +C + IF (LU .LE. 256) GO TO 10 + NERR = NERR + 1 + CALL SETER('DIMENSION OF CUBE EXCEEDS 256',NERR,2) + RETURN + 10 DO 20 J=1,30 + IF (LU .LE. 2**(J-1)) GO TO 30 + 20 CONTINUE + 30 IXBIT = J + NU = MU + NUP2 = NU+2 + NV = MV + NVP2 = NV+2 + NW = MW + NWP2 = NW+2 + FNU = NU + FNV = NV + FNW = NW + SU1 = SU(1.) + SV1 = SV(1.) + SW1 = SW(1.) + SUNU = SU(FNU) + SVNV = SV(FNV) + SWNW = SW(FNW) + AVEU = AVE(SU1,SUNU) + AVEV = AVE(SV1,SVNV) + AVEW = AVE(SW1,SWNW) + EYEU = EYE(1) + EYEV = EYE(2) + EYEW = EYE(3) + NUVWP2 = MUVWP2 + TVAL = TISO + NFLAG = IABS(IFLAG) + IF (NFLAG.EQ.0 .OR. NFLAG.GE.8) NFLAG = 7 +C +C SET UP SCALING +C + FACT = -ISIGN(1,IFLAG) + CALL SET3D (EYE,1.,FNU,1.,FNV,1.,FNW) +C +C BOUND LOWER AND LEFT EDGE OF SLAB +C + EDGE = SIGN(BIG,FACT) + DO 40 IUVW=1,NUVWP2 + SLAB(IUVW,1) = EDGE + SLAB(1,IUVW) = EDGE + 40 CONTINUE +C +C SLICES PERPENDICULAR TO U. THAT IS, V W SLICES. T OF CONSTANT U. +C + IF (NFLAG .LT. 4) GO TO 100 + CALL ZEROSC + ISLBT = -1 +C +C BOUND UPPER AND RIGHT EDGE OF SLAB. +C + DO 50 IV=2,NVP2 + SLAB(IV,NWP2) = EDGE + 50 CONTINUE + DO 60 IW=2,NWP2 + SLAB(NVP2,IW) = EDGE + 60 CONTINUE +C +C GO THRU 3-D ARRAY IN U DIRECTION. IUEW=IU EITHER WAY. +C PICK IU BASED ON EYEU. +C + DO 90 IUEW=1,NU + IU = IUEW + IF (EYEU .GT. AVEU) IU = NU+1-IUEW + U = IU +C +C LOAD THIS SLICE OF T INTO SLAB. +C + DO 80 IV=1,NV + DO 70 IW=1,NW + SLAB(IV+1,IW+1) = T(IU,IV,IW) + 70 CONTINUE + 80 CONTINUE +C +C CONTOUR THIS SLAB. +C + CALL STCNTR (SLAB,NUVWP2,NVP2,NWP2,TVAL) +C +C CONSTRUCT VISIBILITY ARRAY. +C + CALL FILLIN + 90 CONTINUE +C +C SLICES PERPENDICULAR TO V. U W SLICES. T OF CONSTANT V. +C + 100 IF (MOD(NFLAG/2,2) .EQ. 0) GO TO 160 + CALL ZEROSC + ISLBT = 0 +C +C BOUND UPPER AND RIGHT EDGE OF SLAB. +C + DO 110 IU=2,NUP2 + SLAB(IU,NWP2) = EDGE + 110 CONTINUE + DO 120 IW=2,NWP2 + SLAB(NUP2,IW) = EDGE + 120 CONTINUE +C +C GO THRU T IN V DIRECTION. IVEW=IV EITHER WAY. +C + DO 150 IVEW=1,NV + IV = IVEW + IF (EYEV .GT. AVEV) IV = NV+1-IVEW + V = IV +C +C LOAD THIS SLICE OF T INTO SLAB. +C + DO 140 IU=1,NU + DO 130 IW=1,NW + SLAB(IU+1,IW+1) = T(IU,IV,IW) + 130 CONTINUE + 140 CONTINUE +C +C CONTOUR THIS SLAB. +C + CALL STCNTR (SLAB,NUVWP2,NUP2,NWP2,TVAL) +C +C CONSTRUCT VISIBILITY ARRAY. +C + CALL FILLIN + 150 CONTINUE +C +C SLICES PERPENDICULAR TO W. U V SLICES. T OF CONSTANT W. +C + 160 IF (MOD(NFLAG,2) .EQ. 0) GO TO 220 + CALL ZEROSC +C + ISLBT = 1 +C +C BOUND UPPER AND RIGHT EDGE OF SLAB. +C + DO 170 IU=2,NUP2 + SLAB(IU,NVP2) = EDGE + 170 CONTINUE + DO 180 IV=2,NVP2 + SLAB(NUP2,IV) = EDGE + 180 CONTINUE +C +C GO THRU T IN W DIRECTION. +C + DO 210 IWEW=1,NW + IW = IWEW + IF (EYEW .GT. AVEW) IW = NW+1-IWEW + W = IW +C +C LOAD THIS SLICE OF T INTO SLAB. +C + DO 200 IU=1,NU + DO 190 IV=1,NV + SLAB(IU+1,IV+1) = T(IU,IV,IW) + 190 CONTINUE + 200 CONTINUE +C +C CONTOUR THIS SLAB. +C + CALL STCNTR (SLAB,NUVWP2,NUP2,NVP2,TVAL) +C +C CONSTRUCT VISIBILITY ARRAY. +C + CALL FILLIN + 210 CONTINUE +C +C DRAW REFERENCE PLANE EDGES AND W AXIS. +C + 220 IF (IREF .EQ. 0) RETURN + CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2) + IF (EYEV .LT. SV1) GO TO 240 + CALL FRSTC (IFIX(XT),IFIX(YT),1) + DO 230 IU=2,NU + CALL TRN32I (SU(FLOAT(IU)),SV1,SW1,XT,YT,DUM,2) + CALL FRSTC (IFIX(XT),IFIX(YT),2) + 230 CONTINUE + GO TO 250 + 240 CALL PLOTIT (IFIX(XT),IFIX(YT),0) + CALL TRN32I (SUNU,SV1,SW1,XT,YT,DUM,2) + CALL PLOTIT (IFIX(XT),IFIX(YT),1) + 250 IF (EYEU .GT. SUNU) GO TO 270 + CALL FRSTC (IFIX(XT),IFIX(YT),1) + DO 260 IV=2,NV + CALL TRN32I (SUNU,SV(FLOAT(IV)),SW1,XT,YT,DUM,2) + CALL FRSTC (IFIX(XT),IFIX(YT),2) + 260 CONTINUE + GO TO 280 + 270 CALL PLOTIT (IFIX(XT),IFIX(YT),0) + CALL TRN32I (SUNU,SVNV,SW1,XT,YT,DUM,2) + CALL PLOTIT (IFIX(XT),IFIX(YT),1) + 280 IF (EYEV .GT. SVNV) GO TO 300 + CALL FRSTC (IFIX(XT),IFIX(YT),1) + DO 290 IUOW=2,NU + CALL TRN32I (SU(FLOAT(NU-IUOW+1)),SVNV,SW1,XT,YT,DUM,2) + CALL FRSTC (IFIX(XT),IFIX(YT),2) + 290 CONTINUE + GO TO 310 + 300 CALL PLOTIT (IFIX(XT),IFIX(YT),0) + CALL TRN32I (SU1,SVNV,SW1,XT,YT,DUM,2) + CALL PLOTIT (IFIX(XT),IFIX(YT),1) + 310 IF (EYEU .LT. SU1) GO TO 330 + CALL FRSTC (IFIX(XT),IFIX(YT),1) + DO 320 IVOW=2,NV + CALL TRN32I (SU1,SV(FLOAT(NV-IVOW+1)),SW1,XT,YT,DUM,2) + CALL FRSTC (IFIX(XT),IFIX(YT),2) + 320 CONTINUE + GO TO 340 + 330 CALL PLOTIT (IFIX(XT),IFIX(YT),0) + CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2) + CALL PLOTIT (IFIX(XT),IFIX(YT),1) + 340 IF (EYEU.LE.SU1 .OR. EYEV.LE.SV1) GO TO 360 + CALL FRSTC (IFIX(XT),IFIX(YT),1) + DO 350 IW=2,NW + CALL TRN32I (SU1,SV1,SW(FLOAT(IW)),XT,YT,DUM,2) + CALL FRSTC (IFIX(XT),IFIX(YT),2) + 350 CONTINUE +C +NOAO - Plotit buffer needs to be flushed before returning. + call plotit (0, 0, 2) +C -NOAO + RETURN + 360 CALL PLOTIT (IFIX(XT),IFIX(YT),0) + CALL TRN32I (SU1,SV1,SWNW,XT,YT,DUM,2) + CALL PLOTIT (IFIX(XT),IFIX(YT),1) +C +NOAO - Plotit buffer needs to be flushed before returning. + call plotit (0, 0, 2) +C -NOAO + RETURN + END + SUBROUTINE EZISOS (T,MU,MV,MW,EYE,SLAB,TISO) +C + SAVE + DIMENSION T(MU,MV,MW),EYE(3) +C + DATA ANG,PI/.35,3.141592/ +C +C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR +C + CALL Q8QST4 ('NSSL','ISOSRF','EZISOS','VERSION 12') +C +C ARGUMENTS DESCRIBED IN ISOSRF +C +C PICK TYPES OF LINES TO DRAW +C + NU = MU + NV = MV + NW = MW + TVAL = TISO + MAX = MAX0(NU,NV,NW)+2 + ATU = NU/2 + ATV = NV/2 + ATW = NW/2 + EYEU = EYE(1) + EYEV = EYE(2) + EYEW = EYE(3) + RU = EYEU-ATU + RV = EYEV-ATV + RW = EYEW-ATW + RU2 = RU*RU + RV2 = RV*RV + RW2 = RW*RW + DU = SQRT(RV2+RW2) + DV = SQRT(RU2+RW2) + DW = SQRT(RU2+RV2) + DR = 1./SQRT(RU2+RV2+RW2) +C +C COMPUTE THE ARCCOSINE +C + TU = DU*DR + ANGU = ATAN(ABS(SQRT(1.-TU*TU)/TU)) + IF (TU .LE. 0.) ANGU = PI-ANGU + TV = DV*DR + ANGV = ATAN(ABS(SQRT(1.-TV*TV)/TV)) + IF (TV .LE. 0.) ANGV = PI-ANGV + TW = DW*DR + ANGW = ATAN(ABS(SQRT(1.-TW*TW)/TW)) + IF (TW .LE. 0.) ANGW = PI-ANGW +C +C BREAK POINT IS ABOUT 20 DEGREES OR ABOUT .35 RADIANS +C + IFLAG = 0 + IF (ANGU .GT. ANG) IFLAG = IFLAG+4 + IF (ANGV .GT. ANG) IFLAG = IFLAG+2 + IF (ANGW .GT. ANG) IFLAG = IFLAG+1 +C +C FIND SIGN OF IFLAG +C + ICNT = 0 + IF (ABS(RU) .LE. ATU) GO TO 30 + IU = 1 + IF (EYEU .GT. ATU) IU = NU + DO 20 IW=1,NW + DO 10 IV=1,NV + IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2 + ICNT = ICNT+1 + 10 CONTINUE + 20 CONTINUE + 30 IF (ABS(RV) .LE. ATV) GO TO 60 + IV = 1 + IF (EYEV .GT. ATV) IV = NV + DO 50 IW=1,NW + DO 40 IU=1,NU + IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2 + ICNT = ICNT+1 + 40 CONTINUE + 50 CONTINUE + 60 IF (ABS(RW) .LE. ATW) GO TO 90 + IW = 1 + IF (EYEW .GT. ATW) IW = NW + DO 80 IV=1,NV + DO 70 IU=1,NU + IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2 + ICNT = ICNT+1 + 70 CONTINUE + 80 CONTINUE + 90 IFLAG = ISIGN(IFLAG,ICNT) + CALL ISOSRF (T,NU,NU,NV,NV,NW,EYE,MAX,SLAB,TVAL,IFLAG) +C +NOAO - Call to frame is suppressed. +C CALL FRAME +C -NOAO + RETURN + END + SUBROUTINE SET3D (EYE,ULO,UHI,VLO,VHI,WLO,WHI) + SAVE + COMMON /TEMPR/ RZERO +C + DIMENSION EYE(3) +C + COMMON /ISOSR3/ ISCALE ,XMIN ,XMAX ,YMIN , + 1 YMAX ,BIGD ,R0 + COMMON /PWRZ1I/ UUMIN ,UUMAX ,VVMIN ,VVMAX , + 1 WWMIN ,WWMAX ,DELCRT ,EYEU , + 2 EYEV ,EYEW +C +C + AVE(A,B) = (A+B)*.5 +C +C A.S.F. FOR SCALING +C + SU(UTEMP) = UTEMP + SV(VTEMP) = VTEMP + SW(WTEMP) = WTEMP +C +C CONSTANTS FOR PWRZ +C + UUMIN = ULO + UUMAX = UHI + VVMIN = VLO + VVMAX = VHI + WWMIN = WLO + WWMAX = WHI + EYEU = EYE(1) + EYEV = EYE(2) + EYEW = EYE(3) +C +C FIND CORNERS IN 2-SPACE FOR 3-SPACE BOX CONTAINING OBJECT +C + ISCALE = 0 + ATU = AVE(SU(UUMIN),SU(UUMAX)) + ATV = AVE(SV(VVMIN),SV(VVMAX)) + ATW = AVE(SW(WWMIN),SW(WWMAX)) + BIGD = 0. + IF (RZERO .LE. 0.) GO TO 10 +C +C RELETIVE SIZE FEATURE IN USE. +C GENERATE EYE POSITION THAT MAKES BOX HAVE MAXIMUM PROJECTED SIZE. +C + ALPHA = -(VVMIN-ATV)/(UUMIN-ATU) + VVEYE = -RZERO/SQRT(1.+ALPHA*ALPHA) + UUEYE = VVEYE*ALPHA + VVEYE = VVEYE+ATV + UUEYE = UUEYE+ATU + WWEYE = ATW + CALL TRN32I (ATU,ATV,ATW,UUEYE,VVEYE,WWEYE,1) + CALL TRN32I (UUMIN,VVMIN,ATW,XMIN,DUMM,DUMM,2) + CALL TRN32I (UUMAX,VVMIN,WWMIN,DUMM,YMIN,DUMM,2) + CALL TRN32I (UUMAX,VVMAX,ATW,XMAX,DUMM,DUMM,2) + CALL TRN32I (UUMAX,VVMIN,WWMAX,DUMM,YMAX,DUMM,2) + BIGD = SQRT((UUMAX-UUMIN)**2+(VVMAX-VVMIN)**2+(WWMAX-WWMIN)**2)*.5 + R0 = RZERO + GO TO 20 + 10 CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1) + CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMIN),X1,Y1,DUM,2) + CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMAX),X2,Y2,DUM,2) + CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMIN),X3,Y3,DUM,2) + CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMAX),X4,Y4,DUM,2) + CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMIN),X5,Y5,DUM,2) + CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMAX),X6,Y6,DUM,2) + CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMIN),X7,Y7,DUM,2) + CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMAX),X8,Y8,DUM,2) + XMIN = AMIN1(X1,X2,X3,X4,X5,X6,X7,X8) + XMAX = AMAX1(X1,X2,X3,X4,X5,X6,X7,X8) + YMIN = AMIN1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8) + YMAX = AMAX1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8) +C +C ADD RIGHT AMOUNT TO KEEP PICTURE SQUARE +C + 20 WIDTH = XMAX-XMIN + HIGHT = YMAX-YMIN + DIF = .5*(WIDTH-HIGHT) + IF (DIF) 30, 50, 40 + 30 XMIN = XMIN+DIF + XMAX = XMAX-DIF + GO TO 50 + 40 YMIN = YMIN-DIF + YMAX = YMAX+DIF + 50 ISCALE = 1 + CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1) + RETURN + END + SUBROUTINE TRN32I (U,V,W,XT,YT,ZT,IENT) +C +C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR- +C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE +C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15, +C 2, 193-204,1968. +C ARGUMENTS FOR SET +C U,V,W ARE THE 3-SPACE COORDINATES OF THE INTERSECTION +C OF THE LINE OF SIGHT AND THE IMAGE PLANE. THIS +C POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT. +C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION. +C +C TRN32 ARGUMENTS +C U,V,W ARE THE 3-SPACE COORDINATES OF A POINT TO BE +C TRANSFORMED. +C XT,YT THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR- +C MATION. WHEN ISCALE=0, XT AND YT ANR IN THE SAME +C UNITS AS U,V, AND W. WHEN ISCALE'0, XT AND YT +C ARE IN PLOTTER COORDINATES. +C ZT NOT USED. +C + SAVE + COMMON /PWRZ1I/ UUMIN ,UUMAX ,VVMIN ,VVMAX , + 1 WWMIN ,WWMAX ,DELCRT ,EYEU , + 2 EYEV ,EYEW + COMMON /ISOSR3/ ISCALE ,XMIN ,XMAX ,YMIN , + 1 YMAX ,BIGD ,R0 +C +C RANGE OF PLOTTER COORDINATES +C +C +C WARNING +C IF PLOTTER MAXIMUM VALUE RANGES (IN X OR Y DIRECTION) FALL BELOW +C 101, THEN CHANGES MUST BE MADE IN SUBROUTINE FRSTC. THE REQUIRED +C CHANGES ARE MARKED BY WARNING COMMENTS IN FRSTC. + DATA NLX,NBY,NRX,NTY/10,10,32760,32760/ + DATA PI/3.1411592/ +C +C STORE THE PARAMETERS OF THE SET CALL FOR USE +C WITH THE TRANSLATE CALL +C +C DECIDE IF SET OR TRANSLATE CALL +C + IF (IENT .NE. 1) GO TO 50 + AU = U + AV = V + AW = W + EU = XT + EV = YT + EW = ZT +C +C +C +C +C + DU = AU-EU + DV = AV-EV + DW = AW-EW + D = SQRT(DU*DU+DV*DV+DW*DW) + COSAL = DU/D + COSBE = DV/D + COSGA = DW/D +C +C COMPUTE THE ARCCOSINE +C + AL = ATAN(ABS(SQRT(1.-COSAL*COSAL)/COSAL)) + IF (COSAL .LE. 0.) AL = PI-AL + BE = ATAN(ABS(SQRT(1.-COSBE*COSBE)/COSBE)) + IF (COSBE .LE. 0.) BE = PI-BE + GA = ATAN(ABS(SQRT(1.-COSGA*COSGA)/COSGA)) + IF (COSGA .LE. 0.) GA = PI-GA + SINGA = SIN(GA) +C +C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF +C THE 2-SPACE. THE 3-SPACE W AXIS IS TRANSFORMED INTO THE +C 2-SPACE Y AXIS. IF THE LINE OF SIGHT IS CLOSE TO PARALLEL +C TO THE 3-SPACE W AXIS, THE 3-SPACE V AXIS IS CHOSEN (IN- +C STEAD OF THE 3-SPACE W AXIS) TO BE TRANSFORMED INTO THE +C 2-SPACE Y AXIS. +C + ASSIGN 90 TO JDONE + IF (ISCALE) 10, 30, 10 + 10 X0 = XMIN + Y0 = YMIN + X1 = NLX + Y1 = NBY + X2 = NRX-NLX + Y2 = NTY-NBY + X3 = X2/(XMAX-XMIN) + Y3 = Y2/(YMAX-YMIN) + X4 = NRX + Y4 = NTY + FACT = 1. + IF (BIGD .LE. 0.) GO TO 20 + X0 = -BIGD + Y0 = -BIGD + X3 = X2/(2.*BIGD) + Y3 = Y2/(2.*BIGD) + FACT = R0/D + 20 DELCRT = X2 + ASSIGN 80 TO JDONE + 30 IF (SINGA .LT. 0.0001) GO TO 40 + R = 1./SINGA + ASSIGN 70 TO JUMP + RETURN + 40 SINBE = SIN(BE) + R = 1./SINBE + ASSIGN 60 TO JUMP + RETURN +C +C******************** ENTRY TRN32 ************************ +C ENTRY TRN32 (U,V,W,XT,YT,ZT) +C + 50 UU = U + VV = V + WW = W + Q = D/((UU-EU)*COSAL+(VV-EV)*COSBE+(WW-EW)*COSGA) + GO TO JUMP,( 60, 70) + 60 UU = ((EW+Q*(WW-EW)-AW)*COSAL-(EU+Q*(UU-EU)-AU)*COSGA)*R + VV = (EV+Q*(VV-EV)-AV)*R + GO TO JDONE,( 80, 90) + 70 UU = ((EU+Q*(UU-EU)-AU)*COSBE-(EV+Q*(VV-EV)-AV)*COSAL)*R + VV = (EW+Q*(WW-EW)-AW)*R + GO TO JDONE,( 80, 90) + 80 XT = AMIN1(X4,AMAX1(X1,X1+X3*(FACT*UU-X0))) + YT = AMIN1(Y4,AMAX1(Y1,Y1+Y3*(FACT*VV-Y0))) + RETURN + 90 XT = UU + YT = VV + RETURN + END + SUBROUTINE ZEROSC + SAVE +C + COMMON /ISOSR2/ LX ,NX ,NY ,ISCR(8,128), + 1 ISCA(8,128) +C +C ZERO BOTH SCRENE MODELS. +C + DO 20 I=1,LX + DO 10 J=1,NY + ISCR(I,J) = 0 + ISCA(I,J) = 0 + 10 CONTINUE + 20 CONTINUE + RETURN + END + SUBROUTINE STCNTR (Z,L,M,N,CONV) +C + SAVE + DIMENSION Z(L,N) +C +C THIS ROUTINE FINDS THE BEGINNINGS OF ALL CONTOUR LINES AT LEVEL CONV. +C FIRST THE EDGES ARE SEARCHED FOR LINES INTERSECTING THE EDGE (OPEN +C LINES) THEN THE INTERIOR IS SEARCHED FOR LINES WHICH DO NOT INTERSECT +C THE EDGE (CLOSED LINES). BEGINNINGS ARE STORED IN IR TO PREVENT RE- +C TRACING OF LINES. IF IR IS FILLED, THE SEARCH IS STOPPED FOR THIS +C CONV. +C + COMMON /ISOSR6/ IX ,IY ,IDX ,IDY , + 1 IS ,ISS ,NP ,CV , + 2 INX(8) ,INY(8) ,IR(500) ,NR + COMMON /ISOSR7/ IENTRY ,IONES + COMMON /ISOSR9/ BIG ,IXBIT +C +C PACK X AND Y +C + IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1 +C + IENTRY = 0 + NP = 0 + CV = CONV +C +C THE FOLLOWING CODE SHOULD BE RE-ENABLED IF THIS ROUTINE IS USED FOR +C GENERAL CONTOURING +C +C ISS=0 +C DO 2 IP1=2,M +C I=IP1-1 +C IF(Z(I,1).GE.CV.OR.Z(IP1,1).LT.CV) GO TO 1 +C IX=IP1 +C IY=1 +C IDX=-1 +C IDY=0 +C IS=1 +C CALL DRLINE(Z,L,M,N) +C 1 IF(Z(IP1,N).GE.CV.OR.Z(I,N).LT.CV) GO TO 2 +C IX=I +C IY=N +C IDX=1 +C IDY=0 +C IS=5 +C CALL DRLINE(Z,L,M,N) +C 2 CONTINUE +C DO 4 JP1=2,N +C J=JP1-1 +C IF(Z(M,J).GE.CV.OR.Z(M,JP1).LT.CV) GO TO 3 +C IX=M +C IY=JP1 +C IDX=0 +C IDY=-1 +C IS=7 +C CALL DRLINE(Z,L,M,N) +C 3 IF(Z(1,JP1).GE.CV.OR.Z(1,J).LT.CV) GO TO 4 +C IX=1 +C IY=J +C IDX=0 +C IDY=1 +C IS=3 +C CALL DRLINE(Z,L,M,N) +C 4 CONTINUE +C + ISS = 1 + DO 40 JP1=3,N + J = JP1-1 + DO 30 IP1=2,M + I = IP1-1 + IF (Z(I,J).GE.CV .OR. Z(IP1,J).LT.CV) GO TO 30 + IXY = IPXY(IP1,J) + IF (NP .EQ. 0) GO TO 20 + DO 10 K=1,NP + IF (IR(K) .EQ. IXY) GO TO 30 + 10 CONTINUE + 20 NP = NP+1 + IF (NP .GT. NR) RETURN + IR(NP) = IXY + IX = IP1 + IY = J + IDX = -1 + IDY = 0 + IS = 1 + CALL DRCNTR (Z,L,M,N) + 30 CONTINUE + 40 CONTINUE + RETURN + END + SUBROUTINE DRCNTR (Z,L,MM,NN) + SAVE +C + DIMENSION Z(L,NN) +C +C THIS ROUTINE TRACES A CONTOUR LINE WHEN GIVEN THE BEGINNING BY STLINE. +C TRANSFORMATIONS CAN BE ADDED BY DELETING THE STATEMENT FUNCTIONS FOR +C FX AND FY IN DRLINE AND MINMAX AND ADDING EXTERNAL FUNCTIONS. +C X=1. AT Z(1,J), X=FLOAT(M) AT Z(M,J). X TAKES ON NON-INTEGER VALUES. +C Y=1. AT Z(I,1), Y=FLOAT(N) AT Z(I,N). Y TAKES ON NON-INTEGER VALUES. +C + COMMON /ISOSR6/ IX ,IY ,IDX ,IDY , + 1 IS ,ISS ,NP ,CV , + 2 INX(8) ,INY(8) ,IR(500) ,NR + COMMON /ISOSR9/ BIG ,IXBIT +C + LOGICAL IPEN ,IPENO +C + DATA IOFFP,SPVAL/0,0./ + DATA IPEN,IPENO/.TRUE.,.TRUE./ +C +C PACK X AND Y +C + IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1 + FX(X1,Y1) = X1 + FY(X1,Y1) = Y1 + C(P11,P21) = (P11-CV)/(P11-P21) +C + M = MM + N = NN + IF (IOFFP .EQ. 0) GO TO 10 + ASSIGN 100 TO JUMP1 + ASSIGN 150 TO JUMP2 + GO TO 20 + 10 ASSIGN 120 TO JUMP1 + ASSIGN 160 TO JUMP2 + 20 IX0 = IX + IY0 = IY + IS0 = IS + IF (IOFFP .EQ. 0) GO TO 30 + IX2 = IX+INX(IS) + IY2 = IY+INY(IS) + IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL + IPENO = IPEN + 30 IF (IDX .EQ. 0) GO TO 40 + Y = IY + ISUB = IX+IDX + X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX) + GO TO 50 + 40 X = IX + ISUB = IY+IDY + Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY) + 50 IF (IPEN) CALL FRSTS (FX(X,Y),FY(X,Y),1) + 60 IS = IS+1 + IF (IS .GT. 8) IS = IS-8 + IDX = INX(IS) + IDY = INY(IS) + IX2 = IX+IDX + IY2 = IY+IDY + IF (ISS .NE. 0) GO TO 70 + IF (IX2.GT.M .OR. IY2.GT.N .OR. IX2.LT.1 .OR. IY2.LT.1) GO TO 190 + 70 IF (CV-Z(IX2,IY2)) 80, 80, 90 + 80 IS = IS+4 + IX = IX2 + IY = IY2 + GO TO 60 + 90 IF (IS/2*2 .EQ. IS) GO TO 60 + GO TO JUMP1,(100,120) + 100 ISBIG = IS+(8-IS)/6*8 + IX3 = IX+INX(ISBIG-1) + IY3 = IY+INY(ISBIG-1) + IX4 = IX+INX(ISBIG-2) + IY4 = IY+INY(ISBIG-2) + IPENO = IPEN + IF (ISS .NE. 0) GO TO 110 + IF (IX3.GT.M .OR. IY3.GT.N .OR. IX3.LT.1 .OR. IY3.LT.1) GO TO 190 + IF (IX4.GT.M .OR. IY4.GT.N .OR. IX4.LT.1 .OR. IY4.LT.1) GO TO 190 + 110 IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL .AND. + 1 Z(IX3,IY3).NE.SPVAL .AND. Z(IX4,IY4).NE.SPVAL + 120 IF (IDX .EQ. 0) GO TO 130 + Y = IY + ISUB = IX+IDX + X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX) + GO TO 140 + 130 X = IX + ISUB = IY+IDY + Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY) + 140 GO TO JUMP2,(150,160) + 150 IF (.NOT.IPEN) GO TO 170 + IF (IPENO) GO TO 160 +C +C END OF LINE SEGMENT +C + CALL FRSTS (D1,D2,3) + CALL FRSTS (FX(XOLD,YOLD),FY(XOLD,YOLD),1) +C +C CONTINUE LINE SEGMENT +C + 160 CALL FRSTS (FX(X,Y),FY(X,Y),2) + 170 XOLD = X + YOLD = Y + IF (IS .NE. 1) GO TO 180 + NP = NP+1 + IF (NP .GT. NR) GO TO 190 + IR(NP) = IPXY(IX,IY) + 180 IF (ISS .EQ. 0) GO TO 60 + IF (IX.NE.IX0 .OR. IY.NE.IY0 .OR. IS.NE.IS0) GO TO 60 +C +C END OF LINE +C + 190 CALL FRSTS (D1,D2,3) + RETURN + END + SUBROUTINE TR32 (X,Y,MX,MY) + SAVE +C + COMMON /ISOSR1/ ISLBT ,U ,V ,W +C +C A.S.F. FOR SCALING +C + SU(UTEMP) = UTEMP + SV(VTEMP) = VTEMP + SW(WTEMP) = WTEMP +C + XX = X + YY = Y + IF (ISLBT) 10, 20, 30 + 10 CALL TRN32I (SU(U),SV(XX-1.),SW(YY-1.),XT,YT,DUM,2) + GO TO 40 + 20 CALL TRN32I (SU(XX-1.),SV(V),SW(YY-1.),XT,YT,DUM,2) + GO TO 40 + 30 CALL TRN32I (SU(XX-1.),SV(YY-1.),SW(W),XT,YT,DUM,2) + 40 MX = XT + MY = YT + RETURN + END + SUBROUTINE FRSTS (XX,YY,IENT) +C +C THIS IS A SPECIAL VERSION OF THE SMOOTHING DASHED LINE PACKAGE. LINES +C ARE SMOOTHED IN THE SAME WAY, BUT NO SOFTFARE DASHED LINES ARE USED. +C CONDITIONAL PLOTTING ROUTINES ARE CALL WHICH DETERMINE THE VISIBILITY +C OF A LINE SEGMENT BEFORE PLOTTING. +C + SAVE + DIMENSION XSAVE(70) ,YSAVE(70) ,XP(70) ,YP(70) , + 1 TEMP(70) +C + COMMON /ISOSR7/ IENTRY ,IONES +C + DATA NP/150/ + DATA L1/70/ + DATA TENSN/2.5/ + DATA PI/3.14159265358/ + DATA SMALL/128./ +C + AVE(A,B) = .5*(A+B) +C +C DECIDE IF FRSTS,VECTS,LASTS CALL +C + GO TO ( 10, 20, 40),IENT + 10 DEG = 180./PI + X = XX + Y = YY + LASTFL = 0 + SSLP1 = 0.0 + SSLPN = 0.0 + XSVN = 0.0 + YSVN = 0.0 +C +C INITIALIZE THE POINT AND SEGMENT COUNTER +C N COUNTS THE NUMBER OF POINTS/SEGMENT +C + N = 0 +C +C NSEG = 0 FIRST SEGMENT +C NSEG = 1 MORE THAN ONE SEGMENT +C + NSEG = 0 + CALL TR32 (X,Y,MX,MY) +C +C SAVE THE X,Y COORDINATES OF THE FIRST POINT +C XSV1 CONTAINS THE X COORDINATE OF THE FIRST POINT +C OF A LINE +C YSV1 CONTAINS THE Y COORDINATE OF THE FIRST POINT +C OF A LINE +C + XSV1 = MX + YSV1 = MY + GO TO 30 +C +C ************************* ENTRY VECTS ************************* +C ENTRY VECTS (XX,YY) +C + 20 X = XX + Y = YY +C +C VECTS SAVES THE X,Y COORDINATES OF THE ACCEPTED +C POINTS ON A LINE SEGMENT +C + CALL TR32 (X,Y,MX,MY) +C +CIF THE NEW POINT IS TOO CLOSE TO THE PREVIOUS POINT, IGNORE IT +C + IF (ABS(FLOAT(IFIX(XSVN)-MX))+ABS(FLOAT(IFIX(YSVN)-MY)) .LT. + 1 SMALL) RETURN + IFLAG = 0 + 30 N = N+1 +C +C SAVE THE X,Y COORDINATES OF EACH POINT OF THE SEGMENT +C XSAVE THE ARRAY OF X COORDINATES OF LINE SEGMENT +C YSAVE THE ARRAY OF Y COORDINATES OF LINE SEGMENT +C + XSAVE(N) = MX + YSAVE(N) = MY + XSVN = XSAVE(N) + YSVN = YSAVE(N) + IF (N .GE. L1-1) GO TO 50 + RETURN +C +C ************************* ENTRY LASTS ************************* +C ENTRY LASTS +C + 40 LASTFL = 1 +C +C LASTS CHECKS FOR PERIODIC LINES AND SETS UP +C THE CALLS TO KURV1S AND KURV2S +C +C IFLAG = 0 OK TO CALL LASTS DIRECTLY +C IFLAG = 1 LASTS WAS JUST CALLED FROM BY VECTS +C IGNORE CALL TO LASTS +C + IF (IFLAG .EQ. 1) RETURN +C +C COMPARE THE LAST POINT OF SEGMENT WITH FIRST POINT OF LINE +C + 50 IFLAG = 1 +C +C IPRD = 0 PERIODIC LINE +C IPRD = 1 NON-PERIODIC LINE +C + IPRD = 1 + IF (ABS(XSV1-XSVN)+ABS(YSV1-YSVN) .LT. SMALL) IPRD = 0 +C +C TAKE CARE OF THE CASE OF ONLY TWO DISTINCT P0INTS ON A LINE +C + IF (NSEG .GE. 1) GO TO 70 + IF (N-2) 160,150, 60 + 60 IF (N .GE. 4) GO TO 70 + DX = XSAVE(2)-XSAVE(1) + DY = YSAVE(2)-YSAVE(1) + SLOPE = ATAN2(DY,DX)*DEG+90. + IF (SLOPE .GE. 360.) SLOPE = SLOPE-360. + IF (SLOPE .LE. 0.) SLOPE = SLOPE+360. + SLP1 = SLOPE + SLPN = SLOPE + ISLPSW = 0 + SIGMA = TENSN + GO TO 110 + 70 SIGMA = TENSN + IF (IPRD .GE. 1) GO TO 90 + IF (NSEG .GE. 1) GO TO 80 +C +C SET UP FLAGS FOR A 1 SEGMENT, PERIODIC LINE +C + ISLPSW = 4 + XSAVE(N) = XSV1 + YSAVE(N) = YSV1 + GO TO 110 +C +C SET UP FLAGS FOR AN N-SEGMENT, PERIODIC LINE +C + 80 SLP1 = SSLPN + SLPN = SSLP1 + ISLPSW = 0 + GO TO 110 + 90 IF (NSEG .GE. 1) GO TO 100 +C +C SET UP FLAGS FOR THE 1ST SEGMENT OF A NON-PERIODIC LINE +C + ISLPSW = 3 + GO TO 110 +C +C SET UP FLAGS FOR THE NTH SEGMENT OF A NON-PERIODIC LINE +C + 100 SLP1 = SSLPN + ISLPSW = 1 +C +C CALL THE SMOOTHING ROUTINES +C + 110 CALL KURV1S (N,XSAVE,YSAVE,SLP1,SLPN,XP,YP,TEMP,S,SIGMA,ISLPSW) + IF (IPRD.EQ.0 .AND. NSEG.EQ.0 .AND. S.LT.70.) GO TO 170 + IENTRY = 1 +C +C DETERMINE THE NUMBER OF POINTS TO INTERPOLATE FOR EACH SEGMENT +C + IF (NSEG.GE.1 .AND. N.LT.L1-1) GO TO 120 + NPRIME = FLOAT(NP)-(S*FLOAT(NP))/(2.*32768.) + IF (S .GE. 32768.) NPRIME = .5*FLOAT(NP) + NPL = FLOAT(NPRIME)*S/32768. + IF (NPL .LT. 2) NPL = 2 + 120 DT = 1./FLOAT(NPL) + IF (NSEG .LE. 0) CALL FRSTC (IFIX(XSAVE(1)),IFIX(YSAVE(1)),1) + T = 0.0 + NSLPSW = 1 + IF (NSEG .GE. 1) NSLPSW = 0 + NSEG = 1 + CALL KURV2S (T,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP) +C +C SAVE SLOPE AT THE FIRST POINT OF THE LINE +C + IF (NSLPSW .GE. 1) SSLP1 = SLP + NSLPSW = 0 + XSOLD = XSAVE(1) + YSOLD = YSAVE(1) + DO 130 I=1,NPL + T = T+DT + TT = -T + IF (I .EQ. NPL) NSLPSW = 1 + CALL KURV2S (TT,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP) +C +C SAVE THE LAST SLOPE OF THIS LINE SEGMENT +C + IF (NSLPSW .GE. 1) SSLPN = SLP +C +C DRAW EACH PART OF THE LINE SEGMENT +C + CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2) + CALL FRSTC (IFIX(XS),IFIX(YS),2) + XSOLD = XS + YSOLD = YS + 130 CONTINUE + IF (IPRD .NE. 0) GO TO 140 +C +C CONNECT THE LAST POINT WITH THE FIRST POINT OF A PERIODIC LINE +C + CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2) + CALL FRSTC (IFIX(XSV1),IFIX(YSV1),2) +C +C BEGIN THE NEXT LINE SEGMENT WITH THE LAST POINT OF THIS SEGMENT +C + 140 XSAVE(1) = XS + YSAVE(1) = YS + N = 1 + 150 CONTINUE + 160 RETURN + 170 N = 0 + RETURN + END + SUBROUTINE FRSTC (MX,MY,IENT) + SAVE +C + COMMON /ISOSR2/ LX ,NX ,NY ,ISCR(8,128), + 1 ISCA(8,128) + COMMON /ISOSR4/ RX ,RY + COMMON /ISOSR5/ NBPW ,MASK(16) ,GENDON + LOGICAL GENDON + COMMON /ISOSR8/ NMASK(16) ,IXOLD ,IYOLD ,IBTOLD , + 1 HBFLAG ,IOSLSN ,LRLX ,IFSX , + 2 IFSY ,FIRST ,IYDIR ,IHX , + 3 IHB ,IHS ,IHV ,IVOLD , + 4 IVAL ,IHRX ,YCHANG ,ITPD , + 5 IHF + LOGICAL YCHANG ,HBFLAG ,FIRST ,IHF +C +C +C DRAW LINE TO THE POINT MX,MY +C +C ENTER THE POINT INTO THE CURRENT SCREEN, ISCR, IF THE POINT CONFORMS +C TO THE SHADING ALGORITHM. +C THE POINT IS NOT ENTERED WHEN; +C 1. IT IS THE SAME POINT USED IN THE LAST CALL, RESOLUTION PROBLEM +C 2. IT IS PART OF A HORIZONTAL LINE BUT NOT AN END POINT +C 3. THE ENTIRE CONTOUR RESTS ON A HORIZONTAL PLANE +C +C WHEN DRAWING A HORIZONTAL LINE THREE CONDITIONS EXIST; +C 1. WHEN THE LINE IS A HORIZONTAL STEP ENTER ONLY THE OUTSIDE POINT. +C A HORIZONTAL STEP IS DEFINED BY THE ENTERING AND EXITING Y +C DIRECTION THAT IS THE SAME. +C 2. ENTER BOTH END POINTS OF A HORIZONTAL TURNING POINT. A HORIZONTAL +C TURNING POINT IS A LINE WITH GREATER THAN 1 HORIZONTAL BITS +C AND THE ENTERING AND EXITING Y DIRECTION IS DIFFIRENT. +C 3. WHEN THE ENTIRE CONTOUR IS A HORIZONTAL LINE NO POINTS ARE +C ENTERED. THIS CONDITION IS DETECTED BY THE STATUS OF YCHANG. +C IF IT IS TRUE THEN THE CONTOUR IS NOT A SINGLE HORIZONTAL LINE. +C +C THE PREVIOUS POINT IS ERASED IF IT IS A VERTICAL TURNING POINT. +C A VERTICAL TURNING POINT IS A HORIZONTAL LINE WITH ONLY 1 POINT +C AND THE ENTERING AND EXITING Y DIRECTION DIFFERS.THIS DATA IS +C IN THE VARIABLES IOSLSN-OLD SLOPE AND ISLSGN-NEW SLOPE. +C THE CHANGE IN SLOPE MUST BE -1 TO 1 OR 1 TO -1. +C +C OTHERWISE THE POINT IS ENTERED INTO ISCR. +C +C THE TWO ENTRY POINTS ARE REQUIRED BY THE HARDWARE DRAWING ROUTINES. +C FIRSTC IS USED FOR THE FIRST POINT ON THE CONTOUR. THE REMAINING +C POINTS ON THE SAME CONTOUR ARE ENTERED VIA VECTC. +C + DATA IONE/1/ + AVE(A,B) = (A+B)*.5 +C +C COMPUTE VISIBILITY OF THIS POINT +C +C WARNING +C IF X OR Y PLOTTER MAXIMUM VALUE RANGES FALL BELOW 101 THEN THE +C FOLLOWING TWO STATEMENTS WHICH SET IX AND IY MUST BE CHANGED. +C REPLACE THE CONSTANT 1.0 BY 0.5 IN THE STATEMENTS WHERE THE +C MAXIMUM PLOTTER VALUE IS LESS THAN 101 FOR THAT DIRECTION. THE +C PLOTTER CORDINATE RANGES ARE SET IN SET32. +C + IX = FLOAT(MX-1)*RX+1.0 + NRLX = IX + IY = FLOAT(MY-1)*RY+1.0 + IBIT = NBPW-MOD(IX,NBPW) + IX = IX/NBPW+1 + IVNOW = IAND(ISHIFT(ISCA(IX,IY),1-IBIT),IONE) +C +C DECIDE IF FRSTC OR VECTC CALL +C + IF (IENT .NE. 1) GO TO 10 +C + XOLD = MX + YOLD = MY +C +C +C SET INITIAL VALUES +C + IHF = .FALSE. + IYDIR = 0 + ITPD = 0 + IVAL = 0 + IOSLSN = 0 + IFSX = NRLX + IFSY = IY + LASTV = IVNOW + HBFLAG = .FALSE. + YCHANG = .FALSE. + CALL PLOTIT (IFIX(XOLD),IFIX(YOLD),0) + GO TO 180 +C +C**************************** ENTRY VECTC **************************** +C ENTRY VECTC (MX,MY) +C + 10 XNOW = MX + YNOW = MY + JUMP = IVNOW*2+LASTV+1 + GO TO ( 20, 30, 40, 50),JUMP +C +C BOTH VISIBLE +C + 20 CALL PLOTIT (IFIX(XNOW),IFIX(YNOW),1) + GO TO 50 +C +C JUST TURNED VISIBLE +C + 30 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),0) + GO TO 50 +C +C JUST TURNED INVISIBLE +C + 40 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),1) +C +C BOTH INVISIBLE +C + 50 XOLD = XNOW + YOLD = YNOW + LASTV = IVNOW +C +C TEST FOR RESOLUTION PROBLEM +C + IF (NRLX.EQ.LRLX .AND. IY.EQ.IYOLD) RETURN +C +C TEST FOR HORIZONTAL BITS +C + IF (IYOLD .NE. IY) GO TO 70 +C +C HORIZONTAL BITS DETECTED. SET FLAG AND EXIT. +C THIS AND THE NEXT HORIZONTAL BIT TEST IS NECESSARY FOR ISCR TO +C CONFORM TO THE SHADING ALGORITHM IN SUBROUTINE FILLIN +C +C +C IF HORIZONTAL LINE PREVIOUSLY DETECTED EXIT +C + IF (.NOT.HBFLAG) GO TO 60 +C +C IF END OF CONTOUR ON A HORIZONTAL LINE BRANCH FOR SPECIAL PROCESSING. +C + IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210 + GO TO 200 +C +C SAVE SLOPE PRIOR TO HORIZONTAL LINE +C + 60 IHX = IXOLD + IHB = IBTOLD + IHS = IOSLSN + IOSLSN = 0 + HBFLAG = .TRUE. + IHRX = LRLX + IHV = IVOLD + IF (LRLX.EQ.IFSX .AND. IYOLD.EQ.IFSY) IHF = .TRUE. +C +C THIS IS THE SECOND TRAP FOR END OF CONTOUR ON A HORIZONTAL LINE. +C + IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210 + GO TO 200 +C +C COMPUTE THE SLOPE TO THIS POINT +C + 70 IF (IY-IYOLD) 80, 90,100 + 80 ISLSGN = 1 + GO TO 110 + 90 ISLSGN = 0 + GO TO 120 + 100 ISLSGN = -1 + 110 IF (IYDIR .EQ. 0) IYDIR = ISLSGN + 120 CONTINUE +C +C IF PROCESS REACHES THIS CODE THE CONTOUR IS NOT CONTAINED ON A SINGLE +C HORIZONTAL PLANE, SO RECORD THIS FACT BY SETTING Y CHANGE FLAG. +C + YCHANG = .TRUE. +C +C TEST FOR END OF HORIZONTAL LINE +C + IF (.NOT.HBFLAG) GO TO 160 + HBFLAG = .FALSE. +C +C HORIZONTAL LINE JUST ENDED +C +C TEST FOR REDRAW +C + ITEMP = IAND(ISCR(IXOLD,IYOLD),MASK(IBTOLD)) + IF ((IHV .EQ. 0) .AND. (ITEMP .EQ. 0)) GO TO 130 +C +C REDRAWING ERASE THIS POINT +C + ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD)) + ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB)) + GO TO 170 +C +C TEST FOR STEP PROBLEM +C + 130 IF (IHS .NE. ISLSGN) GO TO 140 +C +C STEP PROBLEM +C + GO TO 170 +C +C TURNING PROBLEM HORIZONTAL LINE IS A TURNING POINT +C + 140 CONTINUE +C +C ENTER THE TURNING POINT ONLY IF IT IS NOT THE SECOND SUCCEEDING +C EVENT IN A ROW +C + ICTPD = 1 + IF (IHRX .GT. NRLX) ICTPD = -1 + IF (ICTPD .NE. ITPD) GO TO 150 + ITPD = 0 +C +C ERASE THE FIRST POINT +C + ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB)) + GO TO 170 +C +C ENTER THE TURNING POINT +C + 150 CONTINUE + ITPD = ICTPD +C +C ENTER THE SECOND POINT +C + ISCR(IXOLD,IYOLD) = IOR(ISCR(IXOLD,IYOLD),MASK(IBTOLD)) + GO TO 170 +C +C CHECK IF PREVIOUS ENTRY WAS A VERTICAL TURNING POINT. +C IF SO ERASE IT. +C + 160 IF (ISLSGN.EQ.IOSLSN .OR. (IOSLSN.EQ.0 .OR. ISLSGN.EQ.0)) + 1 GO TO 170 + ITPD = 0 + ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD)) +C + 170 IOSLSN = ISLSGN +C +C CHECK IF THIS GRID POINT PREVIOUSLY ACTIVATED +C + IVAL = IAND(ISCR(IX,IY),MASK(IBIT)) +C +C IF GRID POINTS ACTIVATED BRANCH +C + IF (IVAL .NE. 0) GO TO 190 +C +C GRID POINT NOT ACTIVATED SET AND EXIT +C + 180 CONTINUE + ISCR(IX,IY) = IOR(ISCR(IX,IY),MASK(IBIT)) + GO TO 200 +C +C THIS POINT IS BEING REDRAWN SO ERASE IT. +C (THIS IS TO CONFORM WITH THE SHADING ALGORITHM, FILLIN. +C HOWEVER IF BACK TO STARTING POINT DO NOT ERASE +C + 190 IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) RETURN + ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT)) +C +C + 200 IXOLD = IX + LRLX = NRLX + IYOLD = IY + IBTOLD = IBIT + IVOLD = IVAL + RETURN +C +C PERFORM THIS OPERATION WHEN A CONTOUR STARTS OR ENDS ON A HORIZONTAL +C LINE. +C + 210 CONTINUE +C +C ERASE THE FIRST POINT OF A CONTOUR WHEN IT IS PART OF A HORIZONTAL +C LINE SEGMENT AND IS NOT THE ENDPOINT OF THE SEGMENT +C + IF (.NOT.IHF) GO TO 220 + ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT)) + 220 CONTINUE +C +C ERASE THE FIRST POINT OF A HORIZONTAL LINE SEGMENT WHEN IT ENDS +C THE CONTOUR AND IS NOT THE HIGHEST LINE SEG ON THS SIDE. +C + IF (.NOT.YCHANG) GO TO 230 + IF (IYDIR .NE. IHS) GO TO 200 + 230 ISCR(IHX,IY) = IAND(ISCR(IHX,IY),NMASK(IHB)) + GO TO 200 + END + SUBROUTINE FILLIN +C + SAVE + COMMON /ISOSR2/ LX ,NX ,NY ,ISCR(8,128), + 1 ISCA(8,128) + COMMON /ISOSR5/ NBPW ,MASK(16) ,GENDON + LOGICAL GENDON + COMMON /ISOSR7/ IENTRY ,IONES +C + IF (IENTRY .EQ. 0) RETURN +C +C THIS IS A SHADING ALGORITHM IT IS USED TO DETERMINE CONTOUR LINES +C THAT ARE HIDDEN BY THE PRESENT LINE. THE ALGORITHM PROCESSES +C HORIZONTAL ROWS. IT ASSUMES THAT THE BIT PATTERN PASSED TO IT +C HAS ONLY BITS SET TO MARK THE START AND END OF SHADING. THE +C ALGORITHM ALSO ASSUMES THAT WHEN AN ON BIT IS ENCOUNTERED THAT A +C CORRESPONDING OFF BIT IS INCLUDED IN THE SAME ROW. +C +C +C PULL OUT ROWS OF THE CONTOUR PATTERN +C + IBVAL = 0 + DO 80 IYNOW=1,NY + DO 40 IXNOW=1,LX +C +C IF NO ACTIVATED BITS BRANCH +C + ICRWD = ISCR(IXNOW,IYNOW) + IF (ICRWD .EQ. 0) GO TO 30 +C +C ACTIVATED BITS IN WORD SET SHADING FLAG +C +C CHECK BIT BY BIT FOR ON/OFF FLAGS +C + DO 20 IB=1,NBPW + IBIT = (NBPW+1)-IB +C +C +C PULL OUT THE CURRENT GRID POINT VALUE +C + IVAL = IAND(ICRWD,MASK(IBIT)) +C +C IF IVAL SET, THIS IS AN ON/OFF FLAG +C + IF (IVAL .EQ. 0) GO TO 10 +C +C FLAG BIT, ALWAYS SET +C + IBVAL = MOD(IBVAL+1,2) + GO TO 20 +C +C SHADE THE SCREEN ACCORDING TO THE STATUS OF IBVAL +C + 10 IF (IBVAL .NE. 0) ICRWD = IOR(ICRWD,MASK(IBIT)) +C + 20 CONTINUE +C +C ZERO OUT THE SCREEN +C + ISCR(IXNOW,IYNOW) = 0 + ISCA(IXNOW,IYNOW) = IOR(ICRWD,ISCA(IXNOW,IYNOW)) + GO TO 40 +C + 30 IF (IBVAL .NE. 0) ISCA(IXNOW,IYNOW) = IONES + 40 CONTINUE +C +C FIX FOR NONCORRECTABLE RUNAWAYS +C + IF (IBVAL .EQ. 0) GO TO 80 + IBVAL = 0 + DO 70 K=1,LX + ITEST = 0 + IF (IYNOW .EQ. 1) GO TO 50 + ITEST = ISCA(K,IYNOW-1) + IF (IYNOW .EQ. NY) GO TO 60 + 50 ITEST = IOR(ITEST,ISCA(K,IYNOW+1)) + 60 ISCA(K,IYNOW) = ITEST + 70 CONTINUE +C + 80 CONTINUE + RETURN + END + SUBROUTINE DRAWI (IXA,IYA,IXB,IYB) +C +C INCLUDED FOR USE BY PWRZ +C + SAVE + CALL FRSTC (IXA,IYA,1) + CALL FRSTC (IXB,IYB,2) + RETURN + END + SUBROUTINE MMASK +C +C MAKE THE MACHINE DEPENDENT MASKS USED IN THE CONTOUR DRAWING +C AND SHADING ALGORITHMS +C + SAVE + COMMON /ISOSR5/ NBPW ,MASK(16) ,GENDON + LOGICAL GENDON + COMMON /ISOSR7/ IENTRY ,IONES + COMMON /ISOSR8/ NMASK(16) ,IXOLD ,IYOLD ,IBTOLD , + 1 HBFLAG ,IOSLSN ,LRLX ,IFSX , + 2 IFSY ,FIRST ,IYDIR ,IHX , + 3 IHB ,IHS ,IHV ,IVOLD , + 4 IVAL ,IHRX ,YCHANG ,ITPD , + 5 IHF + COMMON /ISOSR9/ BIG ,IXBIT + LOGICAL YCHANG ,HBFLAG ,FIRST ,IHF + GENDON = .TRUE. + NBPW = 16 +C +C GET BIGGEST REAL NUMBER +C + BIG = R1MACH(2) +C +C MASKS TO SELECT A SPECIFIC BIT +C + DO 10 K=1,NBPW + MASK(K) = ISHIFT(1,K-1) + 10 CONTINUE +C +C GENERATE THE BIT PATTERN 177777 OCTAL +C + ITEMP1 = 0 + ITEMP = MASK(NBPW) + IST = NBPW-1 + DO 20 K=1,IST + ITEMP1 = IOR(ITEMP,ISHIFT(ITEMP1,-1)) + 20 CONTINUE + MFIX = IOR(ITEMP1,1) +C +C MASKS TO CLEAR A SPECIFIC BIT +C + DO 30 K=1,NBPW + NMASK(K) = IAND(ITEMP1,MFIX) + ITEMP1 = IOR(ISHIFT(ITEMP1,1),1) + 30 CONTINUE + IONES = MFIX + RETURN +C +C REVISION HISTORY--- +C +C JANUARY 1978 DELETED REFERENCES TO THE *COSY CARDS AND +C ADDED REVISION HISTORY +C JANUARY 1979 NEW SHADING ALGORITHM +C MARCH 1979 MADE CODE MACHINE INDEPENDENT AND CONFORM +C TO 66 FORTRAN STANDARD +C JUNE 1979 THIS VERSION PLACED ON ULIB. +C SEPTEMBER 1979 FIXED PROBLEM IN EZISOS DEALING WITH +C DETERMINATION OF VISIBILITY OF W PLANE. +C DECEMBER 1979 FIXED PROBLEM WITH PEN DOWN ON CONTOUR +C INITIALIZATION IN SUBROUTINE FRSTC +C MARCH CHANGED ROUTINE NAMES TRN32I AND DRAW TO +C TRN32I AND DRAWI TO BE CONSISTENT WITH THE +C USAGE OF THE NEW ROUTINE PWRZI. +C JUNE 1980 FIXED PROBLEM WITH ZERO INDEX COMPUTATION IN +C SUBROUTINE FRSTC. ADDED INPUT PARAMETER +C DIMENSION STATEMENT MISSING IN EZISOS. +C FIXED ERROR IN COMPUTATION OF ARCCOSINE +C IN EZISOS AND TRN32I. +C DECEMBER 1984 CONVERTED TO GKS LEVEL 0A AND STANDARD FORTRAN 77 +C----------------------------------------------------------------------- +C + END |