C graphics.f Old graphics programs that I wrote from 1970-72. C C ------------------------ C Wm. Randolph Franklin, wrf@ecse.rpi.edu, (518) 276-6077; Fax: -6261 C ECSE Dept., 6026 JEC, Rensselaer Polytechnic Inst, Troy NY, 12180 USA C More info: (1) finger -l wrf@ecse.rpi.edu C (2) http://www.ecse.rpi.edu/homepages/wrf/ C C C>>>AREA AREA 10 C ..................................................................AREA 20 C AREA 30 C FUNCTION AREA AREA 40 C AREA 50 C PURPOSE AREA 60 C TO FIND THE SIGNED AREA OF A POLYGON. AREA 70 C AREA 80 C USAGE AREA 90 C A = AREA (X, Y, N) AREA 100 C AREA 110 C DESCRIPTION OF THE PARAMETERS AREA 120 C X - N LONG VECTOR CONTAINING X-COORDINATES OF THE AREA 130 C VERTICES OF THE POLYGON. AREA 140 C Y - N LONG VECTOR CONTAINING Y-COORDINATES OF THEM. AREA 150 C N - NUMBER OF VERTICES. AREA 160 C AREA 170 C REMARKS AREA 180 C THE FIRST VERTEX MAY OPTIONALLY BE REPEATED. IF SO, N MAY AREA 190 C OPTIONALLY BE INCREASED BY 1. AREA 200 C IF THE VERTICES ARE LISTED CLOCKWISE, THE AREA RETURNED AREA 210 C WILL BE POSITIVE; IF THEY ARE LISTED COUNTERCLOCKWISE, AREA 220 C THE AREA RETURNED WILL BE NEGATIVE. AREA 230 C THUS THIS FUNCTION MAY BE USED TO DETERMINE WHICH WAY AREA 240 C THE VERTICES ARE LISTED. AREA 250 C THE INPUT POLYGON MAY BE A COMPOUND POLYGON, CONSISTING OF AREA 260 C SEVERAL DISJOINT SUBPOLYGONS. IN THIS CASE THE FIRST VERTEX AREA 270 C OF EACH SUBPOLYGON MUST BE REPEATED, AND ALL THESE FIRST AREA 280 C VERTICES MUST BE COUNTED TWICE WHEN CALCULATING N. AREA 290 C IN THIS CASE THE AREA RETURNED WILL BE THE SUM OF THE AREAS AREA 300 C OF THE SEPARATE SUBPOLYGONS. AREA 310 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 9/71. AREA 320 C AREA 330 C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED AREA 340 C NONE AREA 350 C AREA 360 C METHOD AREA 370 C THE POLYGON IS BROKEN INTO QUADRILATERALS WITH A POSSIBLE AREA 380 C EXTRA TRIANGLE, AND THE AREAS OF THESE CALCULATED. AREA 390 C AREA 400 C ..................................................................AREA 410 C AREA 420 FUNCTION AREA(X,Y,N) AREA 430 REAL X(N),Y(N) AREA 440 A=0.0 AREA 450 W=X(1) AREA 460 Z=Y(1) AREA 470 IF (N-3) 4,3,1 AREA 480 1 DO 2 I=4,N,2 AREA 490 A=A+(X(I-1)-W)*(Y(I-2)-Y(I))-(X(I-2)-X(I))*(Y(I-1)-Z) AREA 500 2 CONTINUE AREA 510 3 IF (MOD(N,2).EQ.1) A=A+(X(N)-W)*(Y(N-1)-Z)-(X(N-1)-W)*(Y(N)-Z) AREA 520 A=A/2.0 AREA 530 AREA=A AREA 540 RETURN AREA 550 4 AREA=0.0 AREA 560 RETURN AREA 570 END AREA 580 C>>>GRAV GRAV 10 C ..................................................................GRAV 20 C GRAV 30 C SUBROUTINE CGRAV GRAV 40 C GRAV 50 C PURPOSE GRAV 60 C TO FIND THE CENTROID OF A POLYGON GRAV 70 C GRAV 80 C USAGE GRAV 90 C CALL CGRAV(X,Y,N,PX,PY,A) GRAV 100 C GRAV 110 C DESCRIPTION OF PARAMETERS GRAV 120 C X - N LONG VECTOR CONTAINING X COORDINATES OF THE GRAV 130 C VERTICES OF THE POLYGON. GRAV 140 C Y - N LONG VECTOR CONTAINING Y COORDINATES OF THE GRAV 150 C VERTICES OF THE POLYGON. GRAV 160 C N - NUMBER OF THE VERTICES OF THE POLYGON GRAV 170 C PX - RETURNED X COORDINATE OF THE CENTRIOD GRAV 180 C PY - RETURNED Y COORDINATE OF THE CENTRIOD GRAV 190 C A - RETURNED AREA OF THE POLYGON GRAV 200 C GRAV 210 C REMARKS GRAV 220 C THE FIRST VERTEX MAY BE OPTIONALLY REPEATED. IF SO, GRAV 230 C N MAY BE OPTIONALLY INCREASED BY 1. GRAV 240 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 5/71 GRAV 250 C GRAV 260 C SUBROUTINES AND FUNCTIONS REQUIRED GRAV 270 C NONE GRAV 280 C GRAV 290 C METHOD GRAV 300 C THE POLYGON IS BROKEN INTO TRIANGLES AND THE MOMENT CENTRESGRAV 310 C ARE SUMMED. GRAV 320 C ..................................................................GRAV 330 C GRAV 340 SUBROUTINE CGRAV (X,Y,N,PX,PY,A) GRAV 350 REAL X(N),Y(N),R(200),S(200) GRAV 360 INTEGER O GRAV 370 C O IS OUTPUT UNIT FOR WARNINGS GRAV 380 DATA O/6/ GRAV 390 TOLER=1.E-10 GRAV 400 A=0.0 GRAV 410 PX=0.0 GRAV 420 PY=0.0 GRAV 430 XA=X(1) GRAV 440 YA=Y(1) GRAV 450 R(2)=X(2)-XA GRAV 460 S(2)=Y(2)-YA GRAV 470 IF (N-3) 7,1,2 GRAV 480 1 IF (X(3).EQ.X(1).AND.Y(3).EQ.Y(1)) GO TO 7 GRAV 490 2 DO 3 I=3,N GRAV 500 J=I-1 GRAV 510 R(I)=X(I)-XA GRAV 520 S(I)=Y(I)-YA GRAV 530 B=R(I)*S(J)-R(J)*S(I) GRAV 540 PX=PX+(R(I)+R(J))*B GRAV 550 PY=PY+(S(I)+S(J))*B GRAV 560 A=A+B GRAV 570 3 CONTINUE GRAV 580 IF (ABS(A).LT.TOLER) GO TO 4 GRAV 590 PX=PX/A/3.0+XA GRAV 600 PY=PY/A/3.0+YA GRAV 610 RETURN GRAV 620 4 R(1)=X(1) GRAV 630 R(2)=X(1) GRAV 640 S(1)=Y(1) GRAV 650 S(2)=Y(1) GRAV 660 DO 5 I=2,N GRAV 670 R(1)=AMAX1(R(1),X(I)) GRAV 680 R(2)=AMIN1(R(2),X(I)) GRAV 690 S(1)=AMAX1(S(1),Y(I)) GRAV 700 S(2)=AMIN1(S(2),Y(I)) GRAV 710 5 CONTINUE GRAV 720 PX=(R(1)+R(2))/2.0 GRAV 730 PY=(S(1)+S(2))/2.0 GRAV 740 WRITE (O,6) GRAV 750 6 FORMAT ('0CAUTION: CGRAV OF A POLYGON OF ZERO AREA REQUESTED') GRAV 760 RETURN GRAV 770 7 WRITE (O,8) N GRAV 780 8 FORMAT ('0CAUTION: CGRAV CALLED WITH',I5,' POINTS') GRAV 790 PX=(X(2)+X(1))/2.0 GRAV 800 PY=(Y(2)+Y(1))/2.0 GRAV 810 RETURN GRAV 820 END GRAV 830 C>>>PNP1 PNP1 10 C PNP1 20 C ..................................................................PNP1 30 C PNP1 40 C SUBROUTINE PNPOLY PNP1 50 C PNP1 60 C PURPOSE PNP1 70 C TO DETERMINE WHETHER A POINT IS INSIDE A POLYGON PNP1 80 C PNP1 90 C USAGE PNP1 100 C CALL PNPOLY (PX, PY, XX, YY, N, INOUT ) PNP1 110 C PNP1 120 C DESCRIPTION OF THE PARAMETERS PNP1 130 C PX - X-COORDINATE OF POINT IN QUESTION. PNP1 140 C PY - Y-COORDINATE OF POINT IN QUESTION. PNP1 150 C XX - N LONG VECTOR CONTAINING X-COORDINATES OF PNP1 160 C VERTICES OF POLYGON. PNP1 170 C YY - N LONG VECTOR CONTAING Y-COORDINATES OF PNP1 180 C VERTICES OF POLYGON. PNP1 190 C N - NUMBER OF VERTICES IN THE POLYGON. PNP1 200 C INOUT - THE SIGNAL RETURNED: PNP1 210 C -1 IF THE POINT IS OUTSIDE OF THE POLYGON, PNP1 220 C 0 IF THE POINT IS ON AN EDGE OR AT A VERTEX, PNP1 230 C 1 IF THE POINT IS INSIDE OF THE POLYGON. PNP1 240 C PNP1 250 C REMARKS PNP1 260 C THE VERTICES MAY BE LISTED CLOCKWISE OR ANTICLOCKWISE. PNP1 270 C THE FIRST MAY OPTIONALLY BE REPEATED, IF SO N MAY PNP1 280 C OPTIONALLY BE INCREASED BY 1. PNP1 290 C THE INPUT POLYGON MAY BE A COMPOUND POLYGON CONSISTING PNP1 300 C OF SEVERAL SEPARATE SUBPOLYGONS. IF SO, THE FIRST VERTEX PNP1 310 C OF EACH SUBPOLYGON MUST BE REPEATED, AND WHEN CALCULATING PNP1 320 C N, THESE FIRST VERTICES MUST BE COUNTED TWICE. PNP1 330 C INOUT IS THE ONLY PARAMETER WHOSE VALUE IS CHANGED. PNP1 340 C THE SIZE OF THE ARRAYS X AND Y MUST BE INCREASED IF N > 20. PNP1 350 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 7/70. PNP1 360 C PNP1 370 C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED PNP1 380 C NONE PNP1 390 C PNP1 400 C METHOD PNP1 410 C A VERTICAL LINE IS DRAWN THRU THE POINT IN QUESTION. IF IT PNP1 420 C CROSSES THE POLYGON AN ODD NUMBER OF TIMES, THEN THE PNP1 430 C POINT IS INSIDE OF THE POLYGON. PNP1 440 C PNP1 450 C ..................................................................PNP1 460 C PNP1 470 SUBROUTINE PNPOLY(PX,PY,XX,YY,N,INOUT) PNP1 480 REAL X(020),Y(020),XX(N),YY(N) PNP1 490 LOGICAL MX,MY,NX,NY PNP1 500 DO 1 I=1,N PNP1 510 X(I)=XX(I)-PX PNP1 520 1 Y(I)=YY(I)-PY PNP1 530 INOUT=-1 PNP1 540 DO 2 I=1,N PNP1 550 J=1+MOD(I,N) PNP1 560 MX=X(I).GE.0.0 PNP1 570 NX=X(J).GE.0.0 PNP1 580 MY=Y(I).GE.0.0 PNP1 590 NY=Y(J).GE.0.0 PNP1 600 IF(.NOT.((MY.OR.NY).AND.(MX.OR.NX)).OR.(MX.AND.NX)) GO TO 2 PNP1 610 IF(.NOT.(MY.AND.NY.AND.(MX.OR.NX).AND..NOT.(MX.AND.NX))) GO TO 3 PNP1 620 INOUT=-INOUT PNP1 630 GO TO 2 PNP1 640 3 IF((Y(I)*X(J)-X(I)*Y(J))/(X(J)-X(I))) 2,4,5 PNP1 650 4 INOUT=0 PNP1 660 RETURN PNP1 670 5 INOUT=-INOUT PNP1 680 2 CONTINUE PNP1 690 RETURN PNP1 700 END PNP1 710 C>>>PNP2 PNP2 10 C ..................................................................PNP2 20 C PNP2 30 C SUBROUTINE PNPOLY PNP2 40 C PNP2 50 C PURPOSE PNP2 60 C TO DETERMINE WHETHER A POINT IS INSIDE A POLYGON PNP2 70 C PNP2 80 C USAGE PNP2 90 C CALL PNPOLY (PX, PY, X, Y, N, INOUT ) PNP2 100 C PNP2 110 C DESCRIPTION OF THE PARAMETERS PNP2 120 C PX - X-COORDINATE OF POINT IN QUESTION. PNP2 130 C PY - Y-COORDINATE OF POINT IN QUESTION. PNP2 140 C X - N LONG VECTOR CONTAINING X-COORDINATES OF PNP2 150 C VERTICES OF POLYGON. PNP2 160 C Y - N LONG VECTOR CONTAINING Y-COORDINATES OF PNP2 170 C VERTICES OF POLYGON. PNP2 180 C N - NUMBER OF VERTICES IN THE POLYGON. PNP2 190 C INOUT - THE SIGNAL RETURNED: PNP2 200 C -1 IF THE POINT IS OUTSIDE OF THE POLYGON, PNP2 210 C 0 IF THE POINT IS ON AN EDGE OR AT A VERTEX, PNP2 220 C 1 IF THE POINT IS INSIDE OF THE POLYGON. PNP2 230 C PNP2 240 C REMARKS PNP2 250 C THE VERTICES MAY BE LISTED IN CLOCKWISE OR ANTICLOCKWISE PNP2 260 C ORDER. FOR THIS SUBROUTINE A POINT IS CONSIDERED INSIDE PNP2 270 C THE POLYGON IF IT IS LOCATED IN THE ENCLOSED AREA DEFINED PNP2 280 C BY THE LINE FORMING THE POLYGON. PNP2 290 C THE FIRST POINT MAY OPTIONALLY BE REPEATED, IF SO N MAY PNP2 300 C OPTIONALLY BE INCREASED BY 1. PNP2 310 C THE INPUT POLYGON MAY BE A COMPOUND POLYGON CONSISTING PNP2 320 C OF SEVERAL SEPARATE SUBPOLYGONS. IF SO, THE FIRST VERTEX PNP2 330 C OF EACH SUBPOLYGON MUST BE REPEATED, AND WHEN CALCULATING PNP2 340 C N, THESE FIRST VERTICES MUST BE COUNTED TWICE. PNP2 350 C INOUT IS THE ONLY PARAMETER WHOSE VALUE IS CHANGED. PNP2 360 C PNPOLY CAN HANDLE ANY NUMBER OF VERTICES IN THE POLYGON. PNP2 370 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 6/72. PNP2 380 C PNP2 390 C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED PNP2 400 C NONE PNP2 410 C PNP2 420 C METHOD PNP2 430 C A VERTICAL SEMI-INFINITE LINE IS DRAWN UP FROM THE POINT PNP2 440 C IN QUESTION. IF IT CROSSES THE POLYGON AN ODD NUMBER OF PNP2 450 C TIMES, THEN THE POINT IS INSIDE THE POLYGON. PNP2 460 C PNP2 470 C ..................................................................PNP2 480 C PNP2 490 SUBROUTINE PNPOLY (PX,PY,X,Y,N,INOUT) PNP2 500 REAL X(N),Y(N) PNP2 510 LOGICAL IX,IY,JX,JY,EOR PNP2 520 C EXCLUSIVE OR FUNCTION. PNP2 530 EOR(IX,IY)=(IX.OR.IY).AND..NOT.(IX.AND.IY) PNP2 540 INOUT=-1 PNP2 550 DO 4 I=1,N PNP2 560 XI=X(I)-PX PNP2 570 YI=Y(I)-PY PNP2 580 C CHECK WHETHER THE POINT IN QUESTION IS AT THIS VERTEX. PNP2 590 IF (XI.EQ.0.0.AND.YI.EQ.0.0) GO TO 2 PNP2 600 C J IS NEXT VERTEX NUMBER OF POLYGON. PNP2 610 J=1+MOD(I,N) PNP2 620 XJ=X(J)-PX PNP2 630 YJ=Y(J)-PY PNP2 640 C IS THIS LINE OF 0 LENGTH ? PNP2 650 IF (XI.EQ.XJ.AND.YI.EQ.YJ) GO TO 4 PNP2 660 IX=XI.GE.0.0 PNP2 670 IY=YI.GE.0.0 PNP2 680 JX=XJ.GE.0.0 PNP2 690 JY=YJ.GE.0.0 PNP2 700 C CHECK WHETHER (PX,PY) IS ON VERTICAL SIDE OF POLYGON. PNP2 710 IF (XI.EQ.0.0.AND.XJ.EQ.0.0.AND.EOR(IY,JY)) GO TO 2 PNP2 720 C CHECK WHETHER (PX,PY) IS ON HORIZONTAL SIDE OF POLYGON. PNP2 730 IF (YI.EQ.0.0.AND.YJ.EQ.0.0.AND.EOR(IX,JX)) GO TO 2 PNP2 740 C CHECK WHETHER BOTH ENDS OF THIS SIDE ARE COMPLETELY 1) TO RIGHT PNP2 750 C OF, 2) TO LEFT OF, OR 3) BELOW (PX,PY). PNP2 760 IF (.NOT.((IY.OR.JY).AND.EOR(IX,JX))) GO TO 4 PNP2 770 C DOES THIS SIDE OBVIOUSLY CROSS LINE RISING VERTICALLY FROM (PX,PY)PNP2 780 IF (.NOT.(IY.AND.JY.AND.EOR(IX,JX))) GO TO 1 PNP2 790 INOUT=-INOUT PNP2 800 GO TO 4 PNP2 810 1 IF ((YI*XJ-XI*YJ)/(XJ-XI)) 4,2,3 PNP2 820 2 INOUT=0 PNP2 830 RETURN PNP2 840 3 INOUT=-INOUT PNP2 850 4 CONTINUE PNP2 860 RETURN PNP2 870 END PNP2 880 C>>>ANTB ANTB 10 C ..................................................................ANTB 20 C ANTB 30 C SUBROUTINE ANOTB ANTB 40 C ANTB 50 C PURPOSE ANTB 60 C GIVEN 2 POLYGONS, TO FORM THEIR UNION, INTERSECTION, OR ANY ANTB 70 C ONE OF THE 16 POSSIBLE RESULTANT POLYGONS THAT MAY BE ANTB 80 C FORMED WHEN 2 POLYGONS ARE COMBINED. ANTB 90 C ANTB 100 C USAGE ANTB 110 C CALL ANOTB (AX,AY,NOA,BX,BY,NOB,ABX,ABY,NAB,ITT,LBUG) ANTB 120 C ANTB 130 C DESCRIPTION OF THE PARAMETERS ANTB 140 C AX - NOA LONG VECTOR WITH X-COORDINATES OF VERTICES OF ANTB 150 C POLYGON #1. ANTB 160 C AY - NOA LONG VECTOR WITH Y-COORDINATES OF VERTICES OF ANTB 170 C POLYGON #1. ANTB 180 C NOA - NUMBER OF VERTICES IN POLYGON #1. ANTB 190 C BX - NOB LONG VECTOR WITH X-COORDINATES OF VERTICES OF ANTB 200 C POLYGON #2. ANTB 210 C BY - NOB LONG VECTOR WITH Y-COORDINATES OF VERTICES OF ANTB 220 C POLYGON #2. ANTB 230 C NOB - NUMBER OF VERTICES IN POLYGON #2. ANTB 240 C ABX - NAB LONG VECTOR WITH X-COORDINATES OF VERTICES OF ANTB 250 C RESULTANT POLYGON. ANTB 260 C ABY - NAB LONG VECTOR WITH Y-COORDINATES OF VERTICES OF ANTB 270 C RESULTANT POLYGON. ANTB 280 C ITT CODE FOR OPERATION WANTED ON POLYGONS. LET 'A' ANTB 290 C BE POLYGON #1, 'B' BE POLYGON #2, '.' STAND FOR ANTB 300 C 'AND', '+' STAND FOR 'OR', AND '¬' STAND FOR 'NOT'ANTB 310 C CODE OPERATION ANTB 320 C 1 A.B ANTB 330 C 2 A.¬B ANTB 340 C 3 A ANTB 350 C 4 B.¬A ANTB 360 C 5 B ANTB 370 C 6 (A+B).¬(A.B) THAT IS EXCLUSIVE OR ANTB 380 C 7 A+B ANTB 390 C 8 ¬(A+B) THAT IS NOR ANTB 400 C 9 (A.B)+¬(A+B) THAT IS EXCLUSIVE NOR ANTB 410 C 10 ¬B ANTB 420 C 11 A+¬B ANTB 430 C 12 ¬A ANTB 440 C 13 B+¬A ANTB 450 C 14 ¬(A.B) THAT IS NAND ANTB 460 C 15 THE WHOLE PLANE; THE UNIVERSAL POLYGON ANTB 470 C 16 NOTHING; THE NULL POLYGON ANTB 480 C ANTB 490 C LBUG - A LOGICAL VARIABLE SUPPLIED BY THE MAIN PROGRAM. ANTB 500 C IF .TRUE. VARIOUS DEBUGGING INFO IS PRINTED. ANTB 510 C ANTB 520 C ANTB 530 C REMARKS ANTB 540 C THE INPUT POLYGONS MAY BE COMPOUND; THAT IS THEY MAY BE ANTB 550 C COMPOSED OF SEVERAL DISJOINT SUBPOLYGONS. THE FIRST ANTB 560 C VERTEX OF EACH SUBPOLYGON MUST BE REPEATED AFTER THAT ANTB 570 C POLYGON'S LAST VERTEX. THE VERTICES OF EACH SUBPOLYGON ARE ANTB 580 C TO BE LISTED CLOCKWISE. AN EXCEPTION IS IF YOU WANT A ANTB 590 C POLYGON TO HAVE A HOLE. THEN THE SUBPOLYGON REPRESENTING THEANTB 600 C HOLE WOULD BE LISTED ANTICLOCKWISE. IF THERE WERE AN ANTB 610 C ISLAND IN THE HOLE, IT WOULD BE LISTED CLOCKWISE, ETC. ANTB 620 C THE SUBPOLYGONS OF ONE POLYGON MAY BE LISTED IN ANY ORDER. ANTB 630 C THE SUBPOLYGONS OF ONE POLYGON MUST BE DISJOINT; NO EDGE ANTB 640 C OF ONE SUBPOLYGON MAY CROSS AN EDGE OF ANOTHER SUBPOLYGON ANTB 650 C BELONGING TO THE SAME POLYGON. ANTB 660 C THE SUBPOLYGONS LISTED FOR A POLYGON MUST BE CONSISTENT. ANTB 670 C FOR EXAMPLE, A POLYGON CANNOT HAVE A CLOCKWISE SUBPOLYGON ANTB 680 C INSIDE ANOTHER CLOCKWISE SUBPOLYGON, UNLESS THERE IS A ANTB 690 C COUNTERCLOCKWISE SUBPOLYGON BETWEEN THEM. ANTB 700 C EACH SUBPOLYGON MUST HAVE A NON-ZERO AREA. ANTB 710 C A POLYGON WHICH CONSISTS OF THE WHOLE PLANE EXCEPT FOR ANTB 720 C SOME REGION, IS REPRESENTED BY A COUNTERCLOCKWISE ANTB 730 C SUBPOLYGON. ALSO, THE COMPLEMENT OF A POLYGON IS OBTAINED ANTB 740 C BY MAKING ALL ITS CLOCKWISE SUBPOLYGONS ANTICLOCKWISE, AND ANTB 750 C VICE-VERSA. ANTB 760 C TO PROVIDE COMPLETE GENERALITY, THE NULL POLYGON AND THE ANTB 770 C UNIVERSAL POLYGON MUST BE CONSIDERED. THE FORMER IS ANTB 780 C REPRESENTED BY ONE VERTEX AT (0,0). THE LATTER BY ONE VERTEXANTB 790 C AT (1,1). ANTB 800 C THE OUTPUT POLYGON IS IN EXACTLY THE SAME FORMAT AS THE ANTB 810 C INPUT POLYGONS DESCRIBED ABOVE. ANTB 820 C IF THE INPUT POLYGONS ARE LARGE, VARIOUS ARRAYS IN ANOTB ANTB 830 C MUST BE INCREASED. ANTB 840 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 2/72 ANTB 850 C ANTB 860 C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED ANTB 870 C AREA - FUNCTION TO CALCULATE AREA OF A POLYGON. ANTB 880 C PNPOLY - SUBROUTINE TO TELL IF A POINT IS IN A POLYGON. ANTB 890 C ANTB 900 C METHOD ANTB 910 C SEE FLOWCHART. ANTB 920 C ANTB 930 C ..................................................................ANTB 940 C ANTB 950 SUBROUTINE ANOTB (AX,AY,NOA,BX,BY,NOB,ABX,ABY,NAB,ITT,LBUG) ANTB 960 INTEGER*2 GROUP(20,2),NEWGR(8,2) ,IF(100),NN(2),NGR(2),GO*4,GOB*4 ANTB 970 LOGICAL*1 LCR(20,20),LA,LCR2(20,20),LMEET(8,2) ANTB 980 LOGICAL LBUG ANTB 990 REAL A(2,20,2),EQ(3,20,2),CROSS(2,100),PX(10),PY(10) ANTB1000 REAL AX(NOA),AY(NOA),BX(NOB),BY(NOB),ABX(100),ABY(100),AREAS(2) ANTB1010 INTEGER*2 SAME(20,2),NOLIN(20,2),NOR(10,20,2),INO(10,20,2) ANTB1020 C EXPLANATION OF ABOVE NUMBERS: ANTB1030 C 100 IS THE MAXIMUM NUMBER OF VERTICES IN RESULT. ANTB1040 C 20 IS THE MAXIMUM NUMBER OF VERTICES IN EACH INPUT POLYGON. ANTB1050 C 10 IS THE MAXIMUM NUMBER OF INTERSECTIONS A POLYGON'S SIDE ANTB1060 C MAY HAVE WITH SIDES OF THE OTHER POLYGON. ANTB1070 C 8 IS THE MAXIMUM NUMBER OF SUBPOLYGONS IN EACH INPUT POLYGON. ANTB1080 C WHEN INCREASING THE PROGRAM, CHANGE THESE NUMBERS WHEREEVER THEY ANTB1090 C OCCUR. ANTB1100 C SUBPOLYGONS OF ONE POLYGON MUST NOT INTERSECT. ANTB1110 C THE VARIABLE O IS THE OUTPUT UNIT. ANTB1120 C NGR(I) IS THE NUMBER OF SUBPOLYGONS IN POLYGON #I. ANTB1130 C LMEET(J,I)=.TRUE. IFF SUBPOLYGON #J OF POLYGON #I INTERSECTS ANTB1140 C POLYGON #3-I. ANTB1150 C GROUP(J,I) TELLS WHICH SUBPOLYGON OF POLYGON #I, VERTEX #J OF ANTB1160 C POLYGON #I IS IN. ANTB1170 C NEWGR(J,I) GIVES THE LAST VERTEX IN POLYGON #I WHICH IS IN THAT ANTB1180 C POLYGON'S J-TH SUBPOLYGON. ANTB1190 INTEGER*2 ISTART(16)/1,-1,0,-1,0,0,-1,1,0,0,1,0,1,-1,0,0/ ANTB1200 INTEGER*2 TABYES(2,16)/-1,-1,1,-1,2*30000,-1,1,4*30000,1,1,1,1, ANTB1210 14*30000,-1,1,2*30000,1,-1,-1,-1,4*30000/ ANTB1220 C TABYES IS A DECISION TABLE USED WHEN 2 SUBPOLYGONS INTERSECT. ANTB1230 C THE ENTRIES WITH 30000 SHOULDN'T BE USED. ANTB1240 C TABNO IS THE DECISION TABLE USED FOR A DISJOINT SUBPOLYGON. ANTB1250 C LET I= NUMBER OF POLYGON TO WHICH SUBPOLYGON BELONGS. ANTB1260 C K= 1 IF THIS SUBPOLYGON IS INSIDE POLYGON#(3-I), ELSE 2. ANTB1270 C THEN TABNO(K,I,IT)=1 IF THIS SUBPOLYGON IS TO BE USED, ANTB1280 C 2 IF ITS REVERSE IS TO BE USED, ANTB1290 C 3 IF IT IS NOT TO BE USED AT ALL, ANTB1300 C 4 IF THE RESULT IS THE WHOLE PLANE AS FAR ANTB1310 C AS THIS SUBPOLYGON IS CONCERNED, ANTB1320 C 5 IF THIS ENTRY IS NOT TO BE USED. ANTB1330 INTEGER*2 TABNO(2,2,16)/1,3,1,3,3,1,2,3,4*5,2,3,3,1,4*5,2,1,2,1,4,ANTB1340 11,4,1,3,2,3,2,1,2,1,2,4*5,1,4,4,2,4*5,4,2,1,4,2,4,2,4,8*5/ ANTB1350 C IF THE RESULTANT POLYGON HAS NO VERTICES, NAB=1 IS RETURNED. ANTB1360 C IN THIS CASE, IF THE POLYGON IS NOTHING, THEN ABX(1)=ABY(1)=0.0 ANTB1370 C OR IF THE POLYGON IS THE WHOLE PLANE, ABX(1)=ABY(1)=1.0 ANTB1380 C NOLIN(I,J) : NO OF LINES THAT CROSS THE I-TH LINE OF THE J-TH ANTB1390 C POLYGON. ANTB1400 C INO(K, I, J) : K-TH LINE TO CROSS THE I-TH LINE OF THE J-TH ANTB1410 C POLYGON. ANTB1420 C CROSS(*, NOR (K,I,J) ): COORDINATES OF THE ABOVE INTERSECTION. ANTB1430 C LAST SUBSCRIPT OF AB IS NAB. ANTB1440 C EACH LINE MAY HAVE UP TO 10 INTERSECTIONS WITH OTHER LINES. ANTB1450 C IF THE POLYGONS CONSIST OF MANY SEPARATE SUBPOLYGONS, THEN THESE ANTB1460 C SUBPOLYGONS ARE EACH TREATED SEPARATELY PAIR BY PAIR. ANTB1470 C CODE FOR IT AND ITT: ANTB1480 C 1: A AND B ANTB1490 C 2: A AND NOT B ANTB1500 C 3: A ANTB1510 C 4: B AND NOT A ANTB1520 C 5: B ANTB1530 C 6: A EXCLUSIVE OR B ANTB1540 C 7: A OR B ANTB1550 C 8: NEITHER A NOR B ANTB1560 C 9: NEITHER OR BOTH A AND B ANTB1570 C 10: NOT B ANTB1580 C 11: A OR NOT B ANTB1590 C 12: NOT A ANTB1600 C 13: B OR NOT A ANTB1610 C 14: A NAND B ANTB1620 C 15: THE WHOLE PLANE ANTB1630 C 16: NOTHING ANTB1640 INTEGER O ANTB1650 C O IS OUTPUT UNIT FOR PRINTED MESSAGES ANTB1660 DATA O/6/ ANTB1670 KSA=1 ANTB1680 C IF 2 POINTS ARE CLOSER THAN TOLER, THEY ARE CONSIDERED IDENTICAL. ANTB1690 TOLER=.00001 ANTB1700 MMEET=8 ANTB1710 IT=ITT ANTB1720 LIMIT=100 ANTB1730 NAB=0 ANTB1740 DO 1 I=1,2 ANTB1750 DO 1 J=1,MMEET ANTB1760 1 LMEET(J,I)=.FALSE. ANTB1770 AREAS(1)=AREA(AX,AY,NOA) ANTB1780 AREAS(2)=AREA(BX,BY,NOB) ANTB1790 NCR=0 ANTB1800 C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<0 IF THE VERTEX BEFORE IS THE I+K-TH VERTEX. ANTB2010 C NOTE THAT IT NEVER HAPPENS THAT THE PRECEDING VERTEX IS NOT THE ANTB2020 C I-1 -TH VERTEX AT THE SAME TIME AS THE FOLLOWING VERTEX IS NOT THEANTB2030 C I+1 -TH VERTEX. ANTB2040 IG=0 ANTB2050 C FIND SUBPOLYGONS OF A. ANTB2060 DO 9 I=1,NOA ANTB2070 SAME(I,1)=0 ANTB2080 IF (I.EQ.L) GO TO 8 ANTB2090 IF (ABS(AX(I)-AX(I-1)).LT.TOLER.AND.ABS(AY(I)-AY(I-1)).LT.TOLER) GANTB2100 1O TO 9 ANTB2110 IF (ABS(AX(I)-AX(L)).GT.TOLER.OR.ABS(AY(I)-AY(L)).GT.TOLER) GO TO ANTB2120 18 ANTB2130 SAME(K,1)=NA-K ANTB2140 SAME(NA,1)=K-NA ANTB2150 K=NA+2 ANTB2160 DO 7 J=1,NOB ANTB2170 LCR(NA+1,J)=.FALSE. ANTB2180 7 CONTINUE ANTB2190 L=I+1 ANTB2200 8 NA=NA+1 ANTB2210 IF (K.EQ.NA) IG=IG+1 ANTB2220 GROUP(NA,1)=IG ANTB2230 A(1,NA,1)=AX(I) ANTB2240 A(2,NA,1)=AY(I) ANTB2250 9 CONTINUE ANTB2260 NA=NA-1 ANTB2270 IF (L.EQ.NOA+1) GO TO 11 ANTB2280 IF (LBUG) WRITE (O,10) ANTB2290 10 FORMAT (' FIRST POLYGON IS NOT CLOSED. THE PROGRAM WILL CLOSE IT.'ANTB2300 1) ANTB2310 NA=NA+1 ANTB2320 GROUP(NA,1)=IG ANTB2330 A(1,NA+1,1)=AX(L) ANTB2340 A(2,NA+1,1)=AY(L) ANTB2350 11 CONTINUE ANTB2360 NGR(1)=IG ANTB2370 L=1 ANTB2380 K=1 ANTB2390 NB=0 ANTB2400 IG=0 ANTB2410 C FIND SUBPOLYGONS OF POLYGON B. ANTB2420 DO 14 I=1,NOB ANTB2430 SAME(I,2)=0 ANTB2440 IF (I.EQ.L) GO TO 13 ANTB2450 IF (ABS(BX(I)-BX(I-1)).LT.TOLER.AND.ABS(BY(I)-BY(I-1)).LT.TOLER) GANTB2460 1O TO 14 ANTB2470 IF (ABS(BX(I)-BX(L)).GT.TOLER.OR.ABS(BY(I)-BY(L)).GT.TOLER) GO TO ANTB2480 113 ANTB2490 SAME(K,2)=NB-K ANTB2500 SAME(NB,2)=K-NB ANTB2510 K=NB+2 ANTB2520 DO 12 J=1,NOA ANTB2530 LCR(J,NB+1)=.FALSE. ANTB2540 12 CONTINUE ANTB2550 L=I+1 ANTB2560 13 NB=NB+1 ANTB2570 IF (K.EQ.NB) IG=IG+1 ANTB2580 GROUP(NB,2)=IG ANTB2590 A(1,NB,2)=BX(I) ANTB2600 A(2,NB,2)=BY(I) ANTB2610 14 CONTINUE ANTB2620 NB=NB-1 ANTB2630 IF (L.EQ.NOB+1) GO TO 16 ANTB2640 IF (LBUG) WRITE (O,15) ANTB2650 15 FORMAT (' SECOND POLYGON IS NOT CLOSED. THE PROGRAM WILL CLOSE IT'ANTB2660 1) ANTB2670 NB=NB+1 ANTB2680 GROUP(NB,2)=IG ANTB2690 A(1,NB+1,2)=BX(L) ANTB2700 A(2,NB+1,2)=BY(L) ANTB2710 16 CONTINUE ANTB2720 C NGR(I) IS THE NUMBER OF SUBPOLYGONS IN POLYGON #I. ANTB2730 NGR(2)=IG ANTB2740 NN(1)=NA ANTB2750 NN(2)=NB ANTB2760 ASSIGN 100 TO GO ANTB2770 GO TO (21,21,17,21,18,21,21,21,21,19,21,20,21,21,68,70), IT ANTB2780 C TRIVIAL CASES. ANTB2790 C THE ANSWER IS A. ANTB2800 17 I=1 ANTB2810 M=NA+1 ANTB2820 GO TO 63 ANTB2830 C THE ANSWER IS B. ANTB2840 18 I=2 ANTB2850 M=NB+1 ANTB2860 GO TO 63 ANTB2870 C THE ANSWER IS NOT B. ANTB2880 19 I=2 ANTB2890 M=NB+1 ANTB2900 GO TO 65 ANTB2910 C THE ANSWER IS NOT A. ANTB2920 20 I=1 ANTB2930 M=NA+1 ANTB2940 GO TO 65 ANTB2950 21 CONTINUE ANTB2960 C DOES EITHER POLYGON HAVE NO VERTICES.? ANTB2970 IF (NA.LE.1) GO TO 22 ANTB2980 IF (NB.GT.1) GO TO 26 ANTB2990 IF (A(1,1,2).NE.1.0) GO TO 27 ANTB3000 C B IS THE WHOLE PLANE. ANTB3010 GO TO (17,70,17,20,68,20,68,70,17,70,17,20,68,20,68,70), IT ANTB3020 22 CONTINUE ANTB3030 IF (NB.LE.1) GO TO 23 ANTB3040 C A IS THE WHOLE PLANE. ANTB3050 IF (A(1,1,1).NE.1.0) GO TO 30 ANTB3060 GO TO (18,19,68,70,18,19,68,70,18,19,68,70,18,19,68,70), IT ANTB3070 23 CONTINUE ANTB3080 IF (A(1,1,1).EQ.1.0) GO TO 24 ANTB3090 IF (A(1,1,2).NE.1.0) GO TO 29 ANTB3100 C A IS NOTHING. B IS THE WHOLE PLANE. ANTB3110 GO TO (70,70,70,68,68,68,68,70,70,70,70,68,68,68,68,70), IT ANTB3120 24 IF (A(1,1,2).EQ.1.0) GO TO 25 ANTB3130 C A IS THE WHOLE PLANE. B IS NOTHING. ANTB3140 GO TO (70,68,68,70,70,68,68,70,70,68,68,70,70,68,68,70), IT ANTB3150 C BOTH A AND B ARE THE WHOLE PLANE. ANTB3160 25 GO TO (68,70,68,70,68,70,68,70,68,70,68,70,68,70,68,70), IT ANTB3170 26 CONTINUE ANTB3180 IF (ABS(AREAS(1)).LT.TOLER) GO TO 28 ANTB3190 IF (ABS(AREAS(2)).GE.TOLER) GO TO 31 ANTB3200 27 CONTINUE ANTB3210 C B IS NOTHING. ANTB3220 GO TO (70,17,17,70,70,17,17,20,20,68,68,20,70,68,68,70), IT ANTB3230 28 IF (ABS(AREAS(2)).GE.TOLER) GO TO 30 ANTB3240 29 CONTINUE ANTB3250 C BOTH A AND B ARE NOTHING. ANTB3260 GO TO (70,70,70,70,70,70,70,68,68,68,68,68,68,68,68,70), IT ANTB3270 C A IS NOTHING. ANTB3280 30 GO TO (70,70,70,18,18,18,18,19,19,19,19,68,68,68,68,70), IT ANTB3290 31 CONTINUE ANTB3300 NI=-1 ANTB3310 C PUT VERTICES OF POLYGONS IN CROSS. ANTB3320 DO 33 I=1,2 ANTB3330 N=NN(I) ANTB3340 DO 32 J=1,N ANTB3350 INO(1,J,I)=-1 ANTB3360 INO(2,J,I)=-2 ANTB3370 NOR(1,J,I)=NCR+J ANTB3380 NOR(2,J,I)=NCR+J+1 ANTB3390 CROSS(1,NCR+J)=A(1,J,I) ANTB3400 CROSS(2,NCR+J)=A(2,J,I) ANTB3410 NOLIN(J,I)=2 ANTB3420 32 CONTINUE ANTB3430 CROSS(1,NCR+N+1)=A(1,N+1,I) ANTB3440 CROSS(2,NCR+N+1)=A(2,N+1,I) ANTB3450 NCR=NCR+N+1 ANTB3460 33 CONTINUE ANTB3470 C CALCULATE THE EDGE EQUATIONS ANTB3480 C X*EQ(1,I,K)+Y*EQ(2,I,K)=EQ(3,I,K) FOR THE I-TH SIDE OF THE K-TH ANTB3490 C POLYGON. ANTB3500 DO 37 L=1,2 ANTB3510 N=NN(L) ANTB3520 DO 34 I=1,N ANTB3530 EQ(1,I,L)=A(2,I,L)-A(2,I+1,L) ANTB3540 EQ(2,I,L)=A(1,I+1,L)-A(1,I,L) ANTB3550 EQ(3,I,L)=A(1,I,L)*EQ(1,I,L)+A(2,I,L)*EQ(2,I,L) ANTB3560 C NORMALIZE THE EQUATIONS. ANTB3570 X=SQRT(EQ(1,I,L)**2+EQ(2,I,L)**2) ANTB3580 C AN EDGE OF 0 LENGTH SHOULDN'T HAPPEN. ANTB3590 IF (X.LE.TOLER) GO TO 90 ANTB3600 DO 34 J=1,3 ANTB3610 34 EQ(J,I,L)=EQ(J,I,L)/X ANTB3620 C PUT THE EDGE EQUATIONS IN A FORM SO THAT IF A POINT NEAR AN ANTB3630 C EDGE IS SUBSTITUTED INTO THE EXPRESSION, THE RESULT IS >0 IF THE ANTB3640 C POINT IS IN THE POLYGON, AND <0 IF IT IS OUTSIDE. ANTB3650 DO 36 I=1,N ANTB3660 J=1+I ANTB3670 IF ((A(1,I,L)-A(2,I,L)+A(2,J,L))*EQ(1,I,L)+(A(2,I,L)+A(1,I,L)-A(1,ANTB3680 1J,L))*EQ(2,I,L).GT.EQ(3,I,L)) GO TO 36 ANTB3690 DO 35 K=1,3 ANTB3700 EQ(K,I,L)=-EQ(K,I,L) ANTB3710 35 CONTINUE ANTB3720 36 CONTINUE ANTB3730 37 CONTINUE ANTB3740 C THE FOLLOWING DETERMINES WHICH SIDES OF THE 2 POLYGONS CROSS. ANTB3750 C LCR(I,J)=.TRUE. IF THE I-TH SIDE OF THE 1-ST POLYGON CROSSES THE ANTB3760 C J-TH SIDE OF THE 2-ND POLYGON. ANTB3770 DO 39 I=1,NA ANTB3780 DO 39 J=1,NB ANTB3790 C K TELLS WHETHER EDGE 1,I CROSSES PATH OF EDGE 2,J. ANTB3800 K=(MIN0(0,INT(SIGN(1.0,A(1,I,1)*EQ(1,J,2)+A(2,I,1)*EQ(2,J,2)-EQ(3,ANTB3810 1J,2))))+MIN0(0,INT(SIGN(1.0,A(1,I+1,1)*EQ(1,J,2)+A(2,I+1,1)*EQ(2,JANTB3820 2,2)-EQ(3,J,2)))))+1 ANTB3830 C L TELLS WHETHER EDGE 2,J CROSSES PATH OF EDGE 1,I. ANTB3840 L=MIN0(0,INT(SIGN(1.0,A(1,J,2)*EQ(1,I,1)+A(2,J,2)*EQ(2,I,1)-EQ(3,IANTB3850 1,1))))+MIN0(0,INT(SIGN(1.0,A(1,J+1,2)*EQ(1,I,1)+A(2,J+1,2)*EQ(2,I,ANTB3860 21)-EQ(3,I,1))))+1 ANTB3870 IF (.NOT.LCR(I,J)) GO TO 39 ANTB3880 IF (K.NE.0.OR.L.NE.0) GO TO 38 ANTB3890 NI=I ANTB3900 NJ=J ANTB3910 Y=EQ(1,I,1)*EQ(2,J,2)-EQ(2,I,1)*EQ(1,J,2) ANTB3920 X=(EQ(3,I,1)*EQ(2,J,2)-EQ(2,I,1)*EQ(3,J,2))/Y ANTB3930 Y=(EQ(1,I,1)*EQ(3,J,2)-EQ(3,I,1)*EQ(1,J,2))/Y ANTB3940 K=NOLIN(I,1)+1 ANTB3950 NOLIN(I,1)=K ANTB3960 INO(K,I,1)=J ANTB3970 NCR=NCR+1 ANTB3980 NOR(K,I,1)=NCR ANTB3990 CROSS(1,NCR)=X ANTB4000 CROSS(2,NCR)=Y ANTB4010 K=NOLIN(J,2)+1 ANTB4020 NOLIN(J,2)=K ANTB4030 INO(K,J,2)=I ANTB4040 NOR(K,J,2)=NCR ANTB4050 LMEET(GROUP(I,1),1)=.TRUE. ANTB4060 LMEET(GROUP(J,2),2)=.TRUE. ANTB4070 GO TO 39 ANTB4080 38 LCR(I,J)=.FALSE. ANTB4090 39 CONTINUE ANTB4100 C NI=-1 IF THE 2 POLYGONS DON'T INTERSECT AT ALL. ANTB4110 C FILL ARRAY NEWGR. ANTB4120 C NEWGR(J,I)=LAST VERTEX IN SUBPOLYGON #J OF POLYGON #I. ANTB4130 DO 41 I=1,2 ANTB4140 N=NN(I)+1 ANTB4150 K=1 ANTB4160 DO 40 J=1,N ANTB4170 IF (GROUP(J,I).EQ.K) GO TO 40 ANTB4180 NEWGR(K,I)=J-1 ANTB4190 K=K+1 ANTB4200 40 CONTINUE ANTB4210 NEWGR(K,I)=N ANTB4220 41 CONTINUE ANTB4230 IF (.NOT.LBUG) GO TO 54 ANTB4240 WRITE (O,42) (I,I=1,NB) ANTB4250 42 FORMAT (' LCR',40I3) ANTB4260 DO 43 J=1,NA ANTB4270 WRITE (O,44) J,(LCR(J,I),I=1,NB) ANTB4280 43 CONTINUE ANTB4290 44 FORMAT (I5,(T5,40(1XL2))) ANTB4300 DO 48 L=1,2 ANTB4310 M=NN(L) ANTB4320 WRITE (O,45) L,(SAME(J,L),J=1,M) ANTB4330 45 FORMAT (' POLYGON #',I3/' SAME= ',(20I5)) ANTB4340 WRITE (O,46) (NOLIN(J,L),J=1,M) ANTB4350 46 FORMAT (' NOLIN=',(20I5)) ANTB4360 WRITE (O,47) ((INO(J,K,L),J=1,10),(NOR(J,K,L),J=1,10),K=1,M) ANTB4370 47 FORMAT (' INO:',60X,'NOR:'/(10I5,15X10I5)) ANTB4380 48 CONTINUE ANTB4390 WRITE (O,49) NCR,(J,(CROSS(I,J),I=1,2),J=1,NCR) ANTB4400 49 FORMAT (' NCR,CROSS:',I5/(5(I10,2F5.1))) ANTB4410 DO 50 L=1,2 ANTB4420 N=NN(L) ANTB4430 WRITE (O,51) L,(J,(EQ(I,J,L),I=1,3),(A(I,J,L),I=1,2),GROUP(J,L),J=ANTB4440 11,N) ANTB4450 50 CONTINUE ANTB4460 51 FORMAT ('0EDGE EQUATIONS FOR POLYGON',I3,T70,'VERTEX',T100,'GROUP'ANTB4470 1/(' SIDE',I3,':',G15.5,'*X+',G15.5,'*Y=',G15.5,T65,2G15.5,T100,I3)ANTB4480 2) ANTB4490 DO 52 I=1,2 ANTB4500 L=NGR(I) ANTB4510 WRITE (O,53) I,L,(J,LMEET(J,I),NEWGR(J,I),J=1,L) ANTB4520 52 CONTINUE ANTB4530 53 FORMAT ('0POLYGON #',I1,' HAS',I3,' GROUPS.',8(T30,I2,':',L5,I5/))ANTB4540 C INO(*,I,J) CONTAINS THE LINE NUMBERS OF THE (3-J)-TH POLYGON ANTB4550 C THAT INTERSECT THE I-TH LINE OF THE J-TH POLYGON. ANTB4560 C THESE INTERSECTIONS MUST NOW BE SORTED. ANTB4570 54 CONTINUE ANTB4580 C SOME SUBPOLYGONS MAY BE DISJOINT, THEY DO NOT INTERSECT WITH THE ANTB4590 C OTHER POLYGON. TREAT THEM FIRST. ANTB4600 ASSIGN 61 TO GO ANTB4610 ASSIGN 100 TO GOB ANTB4620 C *****THE FOLLOWING STATEMENTS SIMULATE 2 DO LOOPS SINCE WATFIV ANTB4630 C DOES NOT RECKOGNIZE THE EXTENDED RANGE OF A DO LOOP. ANTB4640 C DO 57 I=1,2 ANTB4650 I=0 ANTB4660 55 I=I+1 ANTB4670 IF (I.GT.2) GO TO 62 ANTB4680 N=NGR(I) ANTB4690 C DO 57 J=1,N ANTB4700 J=0 ANTB4710 56 J=J+1 ANTB4720 IF (J.GT.N) GO TO 55 ANTB4730 IF (LMEET(J,I)) GO TO 61 ANTB4740 KSA=1 ANTB4750 IF (J.GT.1) KSA=NEWGR(J-1,I)+1 ANTB4760 M=NEWGR(J,I)+1-KSA ANTB4770 C KSA IS THE FIRST VERTEX OF SUBPOLYGON #J OF POLYGON #I. ANTB4780 IF (I.EQ.2) GO TO 57 ANTB4790 C I AM TRYING TO SELECT A POINT ON POLYGON A WHICH IS FAR ENOUGH ANTB4800 C FROM AN EDGE OF POLYGON B TO AVOID ROUND OFF ERROR. THIS RANDOM ANTB4810 C POINT SHOULD BE OK. ANTB4820 TEMP1=AX(KSA)*0.1554304+AX(KSA+1)*0.8445696 ANTB4830 TEMP2=AY(KSA)*0.1554304+AY(KSA+1)*0.8445696 ANTB4840 CALL PNPOLY (TEMP1,TEMP2,BX,BY,NOB,K) ANTB4850 GO TO 58 ANTB4860 57 TEMP1=BX(KSA)*0.1554304+BX(KSA+1)*0.8445696 ANTB4870 TEMP2=BY(KSA)*0.1554304+BY(KSA+1)*0.8445696 ANTB4880 CALL PNPOLY (TEMP1,TEMP2,AX,AY,NOA,K) ANTB4890 58 K=TABNO((3-SIGN(1.0,K*AREAS(3-I)))/2,I,IT) ANTB4900 GO TO (63,65,60,59,90), K ANTB4910 59 ASSIGN 67 TO GOB ANTB4920 GO TO 61 ANTB4930 60 ASSIGN 69 TO GOB ANTB4940 61 CONTINUE ANTB4950 GO TO 56 ANTB4960 62 CONTINUE ANTB4970 IF (NI.GT.0) GO TO 71 ANTB4980 GO TO GOB, (100,67,69) ANTB4990 C INCLUDE THIS POLYGON. ANTB5000 63 DO 64 L=1,M ANTB5010 NAB=NAB+1 ANTB5020 ABX(NAB)=A(1,KSA-1+L,I) ANTB5030 ABY(NAB)=A(2,KSA-1+L,I) ANTB5040 64 CONTINUE ANTB5050 GO TO GO, (61,100) ANTB5060 C INCLUDE THE COMPLEMENT OF THIS POLYGON. ANTB5070 65 DO 66 L=1,M ANTB5080 NAB=NAB+1 ANTB5090 ABX(NAB)=A(1,M+KSA-L,I) ANTB5100 ABY(NAB)=A(2,M+KSA-L,I) ANTB5110 66 CONTINUE ANTB5120 GO TO GO, (61,100) ANTB5130 C THE RESULT IS THE WHOLE PLANE. ANTB5140 67 IF (NAB.GT.1) GO TO 100 ANTB5150 68 ABX(1)=1.0 ANTB5160 ABY(1)=1.0 ANTB5170 NAB=1 ANTB5180 GO TO 100 ANTB5190 C THE RESULT IS NOTHING. ANTB5200 69 IF (NAB.GT.1) GO TO 100 ANTB5210 70 ABX(1)=0.0 ANTB5220 ABY(1)=0.0 ANTB5230 NAB=1 ANTB5240 GO TO 100 ANTB5250 C THE 2 POLYGONS INTERSECT. $$$$$ ANTB5260 71 CONTINUE ANTB5270 ASSIGN 100 TO GO ANTB5280 DO 77 I=1,2 ANTB5290 K=NN(I) ANTB5300 DO 76 J=1,K ANTB5310 L=NOLIN(J,I) ANTB5320 DO 72 M=1,L ANTB5330 II=NOR(M,J,I) ANTB5340 PX(M)=CROSS(1,II) ANTB5350 PY(M)=CROSS(2,II) ANTB5360 72 CONTINUE ANTB5370 N=L-1 ANTB5380 LA=PX(1).LT.PX(2)-TOLER.OR.ABS(PX(1)-PX(2)).LT.TOLER.AND.PY(1).LE.ANTB5390 1PY(2) ANTB5400 DO 75 II=1,N ANTB5410 JJ=II+1 ANTB5420 DO 75 IJ=JJ,L ANTB5430 IF (.NOT.LA) GO TO 73 ANTB5440 IF (PX(II).LT.PX(IJ)-TOLER.OR.ABS(PX(II)-PX(IJ)).LT.TOLER.AND.PY(IANTB5450 1I).LE.PY(IJ)) GO TO 75 ANTB5460 GO TO 74 ANTB5470 73 IF (PX(IJ).LT.PX(II)-TOLER.OR.ABS(PX(II)-PX(IJ)).LT.TOLER.AND.PY(IANTB5480 1J).LE.PY(II)) GO TO 75 ANTB5490 74 DEN=PX(II) ANTB5500 PX(II)=PX(IJ) ANTB5510 PX(IJ)=DEN ANTB5520 DEN=PY(II) ANTB5530 PY(II)=PY(IJ) ANTB5540 PY(IJ)=DEN ANTB5550 M=INO(II,J,I) ANTB5560 INO(II,J,I)=INO(IJ,J,I) ANTB5570 INO(IJ,J,I)=M ANTB5580 M=NOR(II,J,I) ANTB5590 NOR(II,J,I)=NOR(IJ,J,I) ANTB5600 NOR(IJ,J,I)=M ANTB5610 75 CONTINUE ANTB5620 76 CONTINUE ANTB5630 77 CONTINUE ANTB5640 IF (.NOT.LBUG) GO TO 79 ANTB5650 DO 78 L=1,2 ANTB5660 M=NN(L) ANTB5670 WRITE (O,45) L,(SAME(J,L),J=1,M) ANTB5680 WRITE (O,46) (NOLIN(J,L),J=1,M) ANTB5690 WRITE (O,47) ((INO(J,K,L),J=1,10),(NOR(J,K,L),J=1,10),K=1,M) ANTB5700 78 CONTINUE ANTB5710 79 CONTINUE ANTB5720 IF (IT.EQ.6) GO TO 80 ANTB5730 IF (IT.NE.9) GO TO 83 ANTB5740 IT=1 ANTB5750 ASSIGN 104 TO GO ANTB5760 GO TO 81 ANTB5770 80 IT=2 ANTB5780 ASSIGN 105 TO GO ANTB5790 81 CONTINUE ANTB5800 DO 82 I=1,NA ANTB5810 DO 82 J=1,NB ANTB5820 82 LCR2(I,J)=LCR(I,J) ANTB5830 NIB=NI ANTB5840 NJB=NJ ANTB5850 83 L=1 ANTB5860 C L IS THE POLYGON WHOSE EDGE I AM ON. ANTB5870 C NI,NJ ARE THE 2 INTERSECTING SIDES. ANTB5880 C SINCE THE DIRECTION OF THE CYCLE IS IMPORTANT, THE FOLLOWING GETS ANTB5890 C ME STARTED IN THE RIGHT DIRECTION. ANTB5900 IF (INT(SIGN(1.0,(2.0*A(1,NI,1)-A(1,NI+1,1))*EQ(1,NJ,2)+(2.0*A(2,NANTB5910 1I,1)-A(2,NI+1,1))*EQ(2,NJ,2)-EQ(3,NJ,2)))*ISTART(IT)) 84,90,88 ANTB5920 84 L=2 ANTB5930 K=NJ ANTB5940 NJ=NI ANTB5950 NI=K ANTB5960 C THE FOLLOWING FINDS 1 CLOSED CYCLE OF VERTICES FOR THE NEW ANTB5970 C POLYGON. ANTB5980 85 IF (NJ.GT.0) GO TO 88 ANTB5990 C THIS NEW POINT IS THE CORNER OF AN OLD POLYGON. ANTB6000 IF (NJ.NE.-2) GO TO 86 ANTB6010 K=SAME(NI,L) ANTB6020 IF (K.GE.0) K=1 ANTB6030 NI=NI+K ANTB6040 86 CONTINUE ANTB6050 NAB=NAB+1 ANTB6060 IF (NAB.EQ.LIMIT) GO TO 97 ANTB6070 IF (NI.LE.0) GO TO 97 ANTB6080 ABX(NAB)=A(1,NI,L) ANTB6090 ABY(NAB)=A(2,NI,L) ANTB6100 IF (NJ.EQ.-1) GO TO 87 ANTB6110 NJ=INO(2,NI,L) ANTB6120 GO TO 85 ANTB6130 87 CONTINUE ANTB6140 K=SAME(NI,L) ANTB6150 IF (K.LE.0) K=-1 ANTB6160 NI=NI+K ANTB6170 NJ=INO(NOLIN(NI,L)-1,NI,L) ANTB6180 GO TO 85 ANTB6190 C THIS NEW POINT IS THE INTERSECTION OF 2 OLD EDGES. ANTB6200 88 L=3-L ANTB6210 N=NOLIN(NJ,L) ANTB6220 DO 89 I=1,N ANTB6230 IF (INO(I,NJ,L).EQ.NI) GO TO 92 ANTB6240 89 CONTINUE ANTB6250 C SOMETHING IS WRONG. ANTB6260 90 WRITE (O,91) I,J,K,L,M,N,II,IJ,JI,JJ,NA,NB,NI,NJ,NN,NO,MM,IT,ITT,KANTB6270 1SA ANTB6280 91 FORMAT ('-********************HELP SUBROUTINE ANOTB IN TROUBLE. ANTB6290 1 ',/,20X,'IT WAS NOT THOUGHT POSSIBLE THAT THIS ERROR COULD ',/,20ANTB6300 2X'OCCUR. PLEASE SEND YOUR LISTING TO R FRANKLIN, U OF O.',/,(14I7ANTB6310 3)) ANTB6320 RETURN ANTB6330 92 M=NOR(I,NJ,L) ANTB6340 NAB=NAB+1 ANTB6350 IF (NAB.EQ.LIMIT) GO TO 97 ANTB6360 ABX(NAB)=CROSS(1,M) ANTB6370 ABY(NAB)=CROSS(2,M) ANTB6380 IF (L.EQ.1) GO TO 93 ANTB6390 IF (.NOT.LCR(NI,NJ)) GO TO 95 ANTB6400 LCR(NI,NJ)=.FALSE. ANTB6410 GO TO 94 ANTB6420 93 IF (.NOT.LCR(NJ,NI)) GO TO 95 ANTB6430 LCR(NJ,NI)=.FALSE. ANTB6440 94 J=NOR(1,NJ,L) ANTB6450 C WHICH WAY TO TURN? ANTB6460 M=SIGN(1.0,CROSS(1,J)*EQ(1,NI,3-L)+CROSS(2,J)*EQ(2,NI,3-L)-EQ(3,NIANTB6470 1,3-L)) ANTB6480 IF (M.EQ.0) GO TO 90 ANTB6490 NI=NJ ANTB6500 NJ=INO(I+M*TABYES(L,IT),NJ,L) ANTB6510 GO TO 85 ANTB6520 C FIND THE START OF ANOTHER CYCLE OF VERTICES. ANTB6530 95 CONTINUE ANTB6540 DO 96 NI=1,NA ANTB6550 DO 96 NJ=1,NB ANTB6560 IF (LCR(NI,NJ)) GO TO 83 ANTB6570 96 CONTINUE ANTB6580 C NEW POLYGON HAS BEEN FOUND COMPLETELY ANTB6590 GO TO GO, (100,104,105) ANTB6600 97 CONTINUE ANTB6610 C THE NEW POLYGON IS TOO LARGE FOR ITS ARRAY. ANTB6620 WRITE (3,98) NAB,I,J,K,L,M,N,NA,NB,NI,NJ,NN,NOA,MOA,NOB,MOB,LIMIT ANTB6630 98 FORMAT (' NEW POLYGON FOUND BY ANOTB TOO LARGE FOR ITS ARRAY'/(10IANTB6640 110)) ANTB6650 WRITE (3,99) LCR,LA ANTB6660 99 FORMAT (50L2) ANTB6670 100 CONTINUE ANTB6680 C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>CHTH CHTH 10 C ..................................................................CHTH 20 C CHTH 30 C SUBROUTINE CRHTCH CHTH 40 C CHTH 50 C PURPOSE CHTH 60 C TO PLOT CROSS-HATCH LINES AT ANY ANGLE AND SPACING INSIDE CHTH 70 C A GIVEN POLYGON WHICH MAY BE COMPOUND. CHTH 80 C CHTH 90 C USAGE CHTH 100 C CALL CRHTCH (X, Y, N, TN, D) CHTH 110 C CHTH 120 C DESCRIPTION OF THE PARAMETERS CHTH 130 C X - A VECTOR, N LONG, WHICH HAS THE X-COORDINATES OF CHTH 140 C THE VERTICES OF THE POLYGON. CHTH 150 C Y - A VECTOR, N LONG, WITH THE Y-COORDINATES. CHTH 160 C N - THE NUMBER OF VERTICES OF THE POLYGON CHTH 170 C TN - THE SLOPE OF THE LINES TO BE DRAWN. IF TN >100. CHTH 180 C THEN THE LINES ARE DRAWN VERTICALLY. CHTH 190 C D - THE DISTANCE BETWEEN THE CROSS-HATCH LINES. CHTH 200 C CHTH 210 C REMARKS CHTH 220 C IF THE POLYGON IS SIMPLE (I.E. ONE CONTINUOUS AREA NOT CHTH 230 C BROKEN UP INTO SEPARATE SUBPOLYGONS), THE FIRST CHTH 240 C VERTEX MAY OPTIONALLY BE REPEATED AT THE END OF X AND Y. CHTH 250 C IF THE POLYGON IS COMPOUND, THE SEPARATE CHTH 260 C SUBPOLYGONS MUST EACH HAVE THEIR FIRST VERTICES REPEATED. CHTH 270 C TO CALCULATE N, COUNT EACH REPEATED VERTEX TWICE. CHTH 280 C IF N > 200, THEN THE SIZE OF PX, PY, AND ORD MUST BE CHTH 290 C INCREASED IN THE SUBROUTINE. CHTH 300 C THE POLYGON IS HATCHED ON ITS INTERIOR WHETHER IT IS CHTH 310 C CLOCKWISE OR COUNTERCLOCKWISE. CHTH 320 C ONE HATCH LINE WILL GO THRU VERTEX #1 OF THE POLYGON. THIS CHTH 330 C MEANS THAT IF THE SAME POLYGON IS CROSS-HATCHED 3 OR MORE CHTH 340 C TIMES WITH D THE SAME BUT TN DIFFERENT, THEN THE CHTH 350 C CROSS-HATCH LINES AT DIFFERENT ANGLES WILL BE CONCURRANT. CHTH 360 C IF NO HATCHES ARE POSSIBLE BECAUSE D IS EITHER TOO LARGE CHTH 370 C OR LESS THAN 0, THEN A MESSAGE IS PRINTED ON UNIT #O. CHTH 380 C TO ACCOMODATE PLOTTERS WITHOUT A FAST SLEWING MODE CHTH 390 C THE PROGRAM DRAWS EACH LINE IN THE OPPOSITE CHTH 400 C DIRECTION TO THE PREVIOUS LINE. CHTH 410 C CRHTCH DOES NOT ALTER ANY OF ITS ARGUMENTS. CHTH 420 C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 4/72 CHTH 430 C CHTH 440 C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED CHTH 450 C LIN - PLOTTER SUBROUTINE TO PLOT A LINE. CHTH 460 C CHTH 470 C METHOD CHTH 480 C SEE FLOWCHART CHTH 490 C CHTH 500 C ..................................................................CHTH 510 C CHTH 520 SUBROUTINE CRHTCH (X,Y,N,TN,D) CHTH 530 REAL X(1),Y(1),XH(2),YH(2),ORD(200),LO,HI,PX(200),PY(200) CHTH 540 INTEGER O CHTH 550 C OUTPUT UNIT FOR PRINTED MESSAGES CHTH 560 DATA O/6/ CHTH 570 IF (TN.GE.100) GO TO 1 CHTH 580 C=1.0/SQRT(1.0+TN*TN) CHTH 590 S=TN*C CHTH 600 GO TO 2 CHTH 610 1 C=0.0 CHTH 620 S=1.0 CHTH 630 2 LO=.1E50 CHTH 640 HI=-.1E50 CHTH 650 C ROTATE POLYGON SO THAT HATCH LINES WILL BE HORIZONTAL. CHTH 660 DO 3 I=1,N CHTH 670 PX(I)=X(I)*C+Y(I)*S CHTH 680 P=-X(I)*S+Y(I)*C CHTH 690 C LO, HI ARE THE LOWEST, HIGHEST Y-COORDINATES OF POLYGON. CHTH 700 IF (P.LT.LO) LO=P CHTH 710 IF (P.GT.HI) HI=P CHTH 720 PY(I)=P CHTH 730 3 CONTINUE CHTH 740 P=AMOD(PY(1)-LO,D)+LO CHTH 750 IF (P.GT.LO) P=P-D CHTH 760 K=(HI-P)/D CHTH 770 IF (K.GT.0.AND.D.GT.0.0) GO TO 7 CHTH 780 WRITE (O,4) LO,HI,D CHTH 790 4 FORMAT (' NO CROSS HATCHES POSSIBLE',3G20.10) CHTH 800 WRITE (O,5) (PX(I),PY(I),I=1,N) CHTH 810 5 FORMAT (4(3X2G15.5)) CHTH 820 WRITE (O,6) K,N,P,LO,HI,D,TN CHTH 830 6 FORMAT (' K=',I10,' N=',I10,' P=',G15.5,' LO=',G15.5,' HI=',G1CHTH 840 15.5,' D=',F10.4,' TN=',F10.4) CHTH 850 RETURN CHTH 860 7 CONTINUE CHTH 870 C K IS THE NUMBER OF HATCH LINES TO BE DRAWN. CHTH 880 DO 15 I=1,K CHTH 890 M=0 CHTH 900 P=P+D CHTH 910 DO 8 L=1,N CHTH 920 J=1+MOD(L,N) CHTH 930 IF (PY(L).LE.P.AND.PY(J).LE.P) GO TO 8 CHTH 940 IF (PY(L).GT.P.AND.PY(J).GT.P) GO TO 8 CHTH 950 M=M+1 CHTH 960 ORD(M)=PX(L)+(PX(J)-PX(L))*(P-PY(L))/(PY(J)-PY(L)) CHTH 970 8 CONTINUE CHTH 980 IF (MOD(M,2).NE.0) WRITE (O,9) CHTH 990 9 FORMAT (' ERROR IN CRHTCH DUE TO ROUND-OFF ERROR.') CHTH1000 IF (M.EQ.0) GO TO 15 CHTH1010 C SORT ORD. CHTH1020 DO 10 J=2,M CHTH1030 NS=M+2-J CHTH1040 DO 10 L=2,NS CHTH1050 IF (ORD(L-1).LE.ORD(L)) GO TO 10 CHTH1060 Q=ORD(L) CHTH1070 ORD(L)=ORD(L-1) CHTH1080 ORD(L-1)=Q CHTH1090 10 CONTINUE CHTH1100 IF (MOD(I,2).EQ.0) GO TO 11 CHTH1110 L=1 CHTH1120 NS=0 CHTH1130 GO TO 12 CHTH1140 11 L=-1 CHTH1150 NS=M+1 CHTH1160 12 DO 14 J=2,M,2 CHTH1170 DO 13 JJ=1,2 CHTH1180 NS=NS+L CHTH1190 XH(JJ)=ORD(NS)*C-P*S CHTH1200 YH(JJ)=ORD(NS)*S+P*C CHTH1210 13 CONTINUE CHTH1220 CALL LIN (XH,YH,2) CHTH1230 14 CONTINUE CHTH1240 15 CONTINUE CHTH1250 RETURN CHTH1260 END CHTH1270