C C $$$$$$$$$$$$$$$$$$$$$ C * * C * PACK * C * * C $$$$$$$$$$$$$$$$$$$$$ C C FOR A FULL DESCRIPTION,SEE T.R. 23 C C PACK REDUCES AND OUTPUTS DATA C PROGRAMME IS DIMENSIONED TO TAKE 5000 DATA POINTS C N IS NUMBER OF POINTS C X,Y COORDS OF N POINTS C DIMENSION X(5000),Y(5000),DX(5000),DY(5000),SEG(20),TAA(2) DIMENSION XX(5000),YY(5000) COMMON/INOUT/ITYP,ITAP,ITAPO COMMON/FEAT/ICDE,ISUBCD,J,ISCO,ISEG COMMON/DAVE/XXD(5000),YYD(5000),IJI IJI=0 IFLAG=1 XMAXD=0. YMAXD=0. XMIND=1000000. YMIND=1000000. C DO 657 JJ=1,20 657 SEG(JJ)=0 C C DELTA IS DIGITISER INCREMENT IN MICRONS I.E. LEAST COUNT C C ERR IS THE ALLOWABLE PLOTTING ERROR AT REDUCED SCALE IN MICRONS C C ISD IS THE DENOMINATOR OF THE SCALE OF THE INPUT DATA C C IOSD IS THE DENOMINATOR OF THE SCALE OF THE OUTPUT DATA C READ(5,104)ERR,DELTA READ(5,100)N READ(5,374)ISD,IOSD 374 FORMAT(2I10) DRAT=IOSD/ISD C C DRAT IS THE RATIO OF INPUT /OUTPUT SCALES C ERR=ERR*DRAT C C ERROR IS NOW THE EQUIVALENT ERROR AT INPUT SCALE C ND=N 100 FORMAT(I4) NNNN=N/5 JD=1 DO 101 JJ=1,NNNN C C READING X AND Y COORDINATES,5 SETS AT A TIME C READ(5,102) X(JD),Y(JD),X(JD+1),Y(JD+1),X(JD+2),Y(JD+2), *X(JD+3),Y(JD+3),X(JD+4),Y(JD+4) YY(JD)=X(JD)/1000. C C SCALING X AND Y COORDINATES,PURELY FOR UNB PLOTTING PURPOSES C XX(JD)=Y(JD)/1000. YY(JD+1)=X(JD+1)/1000. YY(JD+2)=X(JD+2)/1000. YY(JD+3)=X(JD+3)/1000. YY(JD+4)=X(JD+4)/1000. XX(JD+1)=Y(JD+1)/1000. XX(JD+2)=Y(JD+2)/1000. XX(JD+3)=Y(JD+3)/1000. XX(JD+4)=Y(JD+4)/1000. JD=JD+5 101 CONTINUE C C SEARCHING FOR MAXIMUM AND MINIMUM COORDINATES,PURELY FOR UNB PLOTTING C PURPOSES DO 133 JJ=1,N IF(X(JJ).GT.XMAXD)XMAXD=X(JJ) IF(X(JJ).LT.XMIND)XMIND=X(JJ) IF(Y(JJ).GT.YMAXD)YMAXD=Y(JJ) IF(Y(JJ).LT.YMIND)YMIND=Y(JJ) 133 CONTINUE XMAXD=XMAXD/1000. YMAXD=YMAXD/1000. XMIND=XMIND/1000. YMIND=YMIND/1000. XMING=XMIND-2. YMING=YMIND-2. XMAXG=XMAXD+2. YMAXG=YMAXD+2. 999 FORMAT(' ',2F20.5) 183 CONTINUE 104 FORMAT(2F4.0) 102 FORMAT(5(F7.0,1X,F7.0,1X)) WRITE(6,167) 167 FORMAT('1') C C THE PLOTTING SEQUENCE USED ON UNB COMPUTER PLOTTING SYSTEM C C IFLAG=1 GIVES PLOT OF ORIGINAL COORDINATES C IFLAG = 2 GIVES PLOT OF PACKED COORDINATES C CALL DEVICE(611) CALL AREA(2.,2.) CALL GRID(YMING,XMING,YMAXG,XMAXG,1.,1.) IF(IFLAG.EQ.1)CALL PRNTCH('$') IF(IFLAG.NE.1)CALL PRNTCH('*') IF(IFLAG.NE.1)N=IJI CALL SETPLT(YMIND,XMIND,YMAXD,XMAXD) IF(IFLAG.NE.1)CALL NOWPLT(0.,XX(1),YY(1)) IF(IFLAG.NE.1)CALL NOWPLT(1.,XXD(1),YYD(1)) IF(IFLAG.EQ.1)CALL NOWPLT(0.,XX(1),YY(1)) DO 135 IJ=2,N IF(IFLAG.NE.1)CALL NOWPLT(1.,XXD(IJ),YYD(IJ)) IF(IFLAG.EQ.1)CALL NOWPLT(1.,XX(IJ),YY(IJ)) 135 CONTINUE IF(IFLAG.NE.1)CALL NOWPLT(1.,XX(ND),YY(ND)) CALL ENDPLT IF(IFLAG.NE.1)GO TO 7968 C C C DEVICE ASSIGNMENTS AT UNB -- 6 IS THE LINE PRINTER C 5 IS THE CARD READER C 7969 ITYPE=6 ITAP=5 ITAPO=6 C C CODES MOR INDIVIDUAL MAP ELEMENTS. THE PRESENT SYSTEM ASSUMES ALL COORDS C WILL BE OF THE SAME ELEMENT C ICDE=1 ISUBCD=1 ISCO=1 ISCI=1 ISEG=1 ISCOF=ISCO-1 C C ISCOF FORMER NUMBER OF OUTPUT POINTS C EPS=ERR C C EPS IS NUMBER OF DIGITIZER STEPS C EPSI=1.5*EPS EPSE=8.0*EPS C2=13.615*EPS*EPS C1=8.000894*EPS C3=2.82*EPS C C COEFFICIENTS OF HYPERBOLA C NSEG=0 IB=1 XP=X(1) YP=Y(1) 1 X1=XP Y1=YP J=1 IF(IB.EQ.N)IFLAG=2 IF(IB.EQ.N)GO TO 183 DO 2 I=1,2 DXX=X(IB+I)-X(IB+I-1) DYY=Y(IB+I)-Y(IB+I-1) DS=1.0/SQRT(DXX*DXX+DYY*DYY) 2 TAA(I)=DS*DXX C C TAA IS COS C TA=0.5*(TAA(1)+TAA(2)) C C TA AVERAGE OF 2 COSINES C TB=SQRT(1-TA*TA) C C TB CORRESPONDING AVERAGE SINE C DO 3 NN=IB,N IF(ABS(TA*(Y(NN)-YP)-TB*(X(NN)-XP)).GT.EPSE) GO TO 4 C C ARE POINTS IN TUBE EPSE WIDE C 3 CONTINUE NN=N 4 L=1 ST=0.2 5 DXX=X(NN)-XP DYY=Y(NN)-YP D=1.0/SQRT(DXX*DXX+DYY*DYY) C C D IS 1/SEGMENT LENGTH C TX=D*DYY C C TX IS SIN THETA C TY=D*DXX C C TY IS COS THETA C DO 6 I=IB,NN DNN=I-IB+1 IF(ABS(TX*(X(I)-XP)-TY*(Y(I)-YP)).GT.EPS) GO TO 16 C C CHECK RIGOROUSLY IF POINTS IN TUBE EPS WIDE C 6 CONTINUE IF(NN.GE.N) GO TO 30 IF(L.EQ.1) GO TO 18 IF(ABS(ST*DNN).LE.1.0) GO TO 30 19 ST=-0.5*ST L=-L 18 IF(ABS(ST*DNN).LE.1) GO TO 17 NN=NN+IFIX(ST*DNN) IF(NN.LE.N)GO TO 14 NN=N GO TO 5 17 NN=NN+L 14 IF(NN.GE.IB) GO TO 5 NN=IB GO TO 5 16 IF(L.EQ.1) GO TO 19 GO TO 18 30 SEG(J)=1.0/D IF(NN.LT.N) GO TO 7 C C OPTIONAL PRINTOUT OF ORIGINAL COORDINATES C DO 962 I=1,N 963 FORMAT(' ',50X,2F12.5) 962 CONTINUE C CALL REDOUT(J,X1,Y1,SEG,X(N),Y(N),EPS) C C UPLOT CALLED TO REDUCE PACKED DATA C CALL UPLOT(J,X1,Y1,SEG,X(N),Y(N),EPS) IF(IB.EQ.N)IFLAG=2 IF(IB.EQ.N)GO TO 183 ISCOF=ISCO-ISCOF C C ISCOF NOW NUMBER OF OUTPUT POINTS IN THIS BLOCK C ITYP=6 C C PRESENT SYSTEM IS SET UP TO DOUBLE ERROR EPSILON AND RE-ITERATE C IF ONLY ONE PACKED SET OF DATA IS WANTED(I.E. TO THE ORIGINAL ERROR SPECS C READ IN),REPLACE THE TWO FOLLOWING CARDS WITH ONE CARD ---STOP C IFLAG=2 GO TO 183 7 IB=NN+1 XPP=XP YPP=YP XP=X(NN) YP=Y(NN) 8 TI=D*(XP-XPP) C C COS OF LINE SEGMENT C T2=D*(YP-YPP) C C SIN OF LINE SEGMENT C DO 10 I=IB,N C C CHECK POINTS BEYOND SEGMENT C DX(I)=TI*(X(I)-XP)+T2*(Y(I)-YP) C C DX DISPLACEMENT IN LINE DIRECTION C DY(I)=TI*(Y(I)-YP)-T2*(X(I)-XP) C C DY DISPLACEMENT IN DIRECTION NORMAL TO LINE C IF(DX(I).LE.0.0) GO TO 9 C C IF LINE TURNS BACK PAST FIRST POINT THEN CORNER POINT C IF(ABS(DY(I))*(DX(I)+C3)-(C1*DX(I)+C2).GT.0.0) GO TO 11 C C CHECK IF DY STILL WITHIN HYPERBOLA C 10 CONTINUE GO TO 9 11 C1B=SIGN(C1,DY(I)) C2B=SIGN(C2,DY(I)) C C C1B,C2B SIGNS SAME AS DY C DDX=DX(I)-DX(I-1) C C DDX DIFFERENCE BETWEEN LAST 2DX'S C IF(ABS(DDX).GT.2.0E-8) GO TO 12 C C IS HYPERBOLA CUT BY NEARLY NORMAL LINE TO AXIS C DXI=DX(I) GO TO 13 12 AL=(DY(I)-DY(I-1))/DDX C C AL IS SLOPE OF LAST TWO POINTS C Q=DY(I)-AL*DX(I) C C Q IS Y INTERCEPT OF STRAIGHT LINE BETWEEN LAST 2 POINTS C CCC=C1B*C3-C2B DXBB=DX(I) DXI=DX(I-1) DB1=(C1B*DXBB+C2B)/(DXBB+C3) DEN1=DXBB+C3 DEN2=C3*DEN1 31 DXIB=DXI ALB=CCC/(DXIB*DEN1+DEN2) QB=QB1-ALB*DXBB DXI=(Q-QB)/(ALB-AL) IF(ABS(DXI-DXIB).GT.0.5*DELTA) GO TO 31 13 DYI=(C1B*DXI+C2B)/(DXI+C3) SEM=SQRT(DXI*DXI+DYI*DYI) TX=DYI/SEM TY=DXI/SEM II=I-1 DO 15 L=IB,II IF(ABS(TX*DX(L)-TY*DY(L)).GT.EPS) GO TO 9 15 CONTINUE IB=1 J=J+1 SEG(J)=SIGN(SEM,DYI) XPP=XP YPP=YP D=1.0/SEM DDXI=DXI-DX(I) DDYI=DYI-DY(I) XP=X(I)+T1*DDXI-T2*DDYI YP=Y(I)+T1*DDYI+T2*DDXI GO TO 8 9 CALL REDOUT(J,X1,Y1,SEG,XP,YP,EPS) C C C UPLOT CALLED TO REDUCE PACKED DATA C CALL UPLOT(J,X1,Y1,SEG,XP,YP,EPS) IF(IB.EQ.N)IFLAG=2 IF(IB.EQ.N)GO TO 183 C C OUTPUT BLOCK OF REDUCED DATA C GO TO 1 7968 IFLAG=1 C C CALCULATION OF PACKING FACTOR,AND PRINTOUT OF ERRORS AND SCALES C FACTP=FLOAT(ND)/(FLOAT(IJI)+2.0) WRITE(6,5432)FACTP 5432 FORMAT(' ',20X,'PACKING FACTOR = ',F20.2) XERR=ERR/DRAT WRITE(6,5434)XERR 5434 FORMAT(' ',20X,'ERROR AT REDUCED SCALE IN MICRONS = ',F20.10) 5433 FORMAT(' ',20X,' ERROR AT ORIGINAL SCALE = ',F20.10) WRITE(6,5433)ERR WRITE(6,5435)ISD,IOSD 5435 FORMAT(' ','INPUT DENOMINATOR = ',I10,' OUTPUT DENOMINATOR = ', *I10) C C THIS SECTION DOUBLES THE ERROR TUBE AND REPEATS THE PACKING PROCEDURE C ERR=ERR*2. IJI=0 IB=N-1 N=ND IF(ERR.GT.10000.)STOP C MAXIMUM ERROR ALLOWED = 10000 MICROMETERS C GO TO 7969 40 FORMAT(' ',2X,6(5X,I4),5X,F8.2) END SUBROUTINE REDOUT(J,X1,Y1,SEG,XP,YP,EPS) C C REDOUT OUTPUTS RESULTS OF DATRED ONTO ITAPO AFTER REDUCTION C J NUMBER OF SEGMENTS IN THIS RECORD C X1,Y1 INITIAL POINT IN THIS RECORD C SEG ARRAY OF SEGMENT LENGTHS... J OF THEM C XP,YP FINAL POINT IN THIS RECORD C EPS TOLERANCE SPECIFIED FOR THIS REDUCTION C DIMENSION SEG(J),ISSEG(20) C EQUIVALENCE OUTPUTS TO INTEGERS COMMON/INOUT/ITP,ITAP,ITAPO COMMON/FEAT/ICDE,ISUBCD,ISCI,ISCO,ISEG ITAPO=6 C C IF INTERMEDIATE STORAGE OF DATA REQUIRED,CHANGE UNIT NUMBER C C C ICDE,ISUBCD ARE CODE AND SUBCODE OF FEATURE C J NUMBER OF SEGMENTS IN THIS RECORD C ISCO OUTPUT RECORD NUMBER C IF(ISCO .GT. 1)GO TO 1 ISEG=0 C C ISEG NUMBER OF SEGMENT LENGTHS OUTPUT PER FEATURE C IX1=X1+.5 IY1=Y1+.5 C C INITIAL POINT MADE AN INTEGER C WRITE(ITAPO,40)IX1,IY1,ICDE,ISUBCD IEPS=EPS+.5 C C SET ERROR AS INTEGER C WRITE(ITAPO,41)IEPS C C IF FIRST RECORD OF FEATURE WRITE FIRST POINT X1,Y1 AND EPS C 1 IXP=XP+.5 IYP=YP+.5 C C MAKE INTEGERS OF FINAL POINTS C WRITE(ITAPO,40)IXP,IYP,ISCO,J C C WRITE FINAL POINT XP,YP RECORD ISCO AND NO SEG LENGTHS J C ISCO=ISCO+1 IF(J .EQ. 1) RETURN DO 100 I=1,J ISSEG(I)=SEG(I)+.5 100 CONTINUE C C CHANGE SEGMENT LENGTHS TO INTEGERS C ISEG=ISEG+J WRITE(ITAPO,42)(ISSEG(I),I=1,J) C C WRITE J SEGMENT LENGTHS C RETURN 40 FORMAT(' ',2I6,2I4) 41 FORMAT(' ',I6) 42 FORMAT(' ',20X,'SEGMENT LENGTHS ARE',20X//20I6) END SUBROUTINE UPLOT(M,XI,YI,S,XN,YN,E) C C SUBROUTINE TO UNPACK PACKED DATA C DIMENSION S(20),GX(22),GY(22) COMMON/INOUT/ITYP,ITAP,ITAPO COMMON/DAVE/XXD(5000),YYD(5000),IJI DATA STI/1.0000/ C C IF A CONVERSION OF UNITS IS REQUIRED,CHANGE STI C IF(M .LE. 0) RETURN C C M LE 0 FOR DUMMY CALL C IF(M.GE.2) GO TO 200 C C FOR SINGLE SEGMENT NEED ONLY PLOT END POINTS C EX=XN*STI EY=YN*STI WRITE(6,199)EX,EY IJI=IJI+1 XXD(IJI)=EY/1000. YYD(IJI)=EX/1000. C C SCALING FOR UNB PLOTTING PURPOSES ONLY C IFLAG=2 199 FORMAT(' ',4X,'UPLOT COORDINATES ARE ',2F20.10) RETURN C C CALCULATION OF UNPACKING C 200 GX(1)=0.0 GY(1)=0.0 GX(2)=S(1) GY(2)=0.0 C01=8.000894*E C02=13.615*E*E C03=2.82*E C C COEFFS OF HYPERBOLA C A0=C02*C02 A1=2*C01*C02 A2=C01*C01 V=E/SQRT(FLOAT(M)) DO 20 I=2,M SS=S(I)*S(I) SI=ABS(S(I)) B0=C03*C03*SI B1=2*C03*SI+C03*C03 B2=2*C03+SI X=ABS(S(I)) 21 P2=((A2*X)+A1)*X+A0 P3=((X+B2)*X+B1)*X+B0 X=SI-P2/P3 Y=(C01*X+C02)/(X+C03) IF(ABS(X*X+Y*Y-SS) .GT. ABS(S(I)*V)) GO TO 21 Y=SQRT(SS-X*X) Y=SIGN(Y,S(I)) T1=(GX(I)-GX(I-1))/ABS(S(I-1)) T2=(GY(I)-GY(I-1))/ABS(S(I-1)) C T1,T2 COS, SIN OF LINE SEGMENT ANGLE REL TO AXIS GX(I+1)=GX(I)+T1*X-T2*Y 20 GY(I+1)=GY(I)+T2*X+T1*Y C C GX,GY REL COORDS OF END PT OF LINE SEG C 24 MM=M+1 D=1.0/(GX(MM)*GX(MM)+GY(MM)*GY(MM)) C1=(GY(MM)*(YN-YI)+GX(MM)*(XN-XI))*D C2=(-GY(MM)*(XN-XI)+GX(MM)*(YN-YI))*D C C EVERYTHING EXPRESSED IN DIGITIZER STEPS C DO 41 I=2,MM EX=(XI+C1*GX(I)-C2*GY(I))*STI EY=(YI+C2*GX(I)+C1*GY(I))*STI C C COMPUTE END POINTS FOR EACH LINE SEGMENT C WRITE(6,199)EX,EY IJI=IJI+1 XXD(IJI)=EY/1000. YYD(IJI)=EX/1000. C C SCALING FOR UNB PLOTTING PURPOSES ONLY C IFLAG=2 C DRAW LINE SEGMENTS 41 CONTINUE RETURN END