perm filename BUG.SAI[S,AIL]1 blob sn#028537 filedate 1973-03-12 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00013 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001	   VALID 00019 PAGES
C00003 00002	BEGIN "CURVE" 
C00011 00003	⊃ Following are the macros to unpack,packed info in arrays.
C00013 00004	⊃ OUTER BLOCK STARTS HERE 
C00015 00005	SIMPLE INTEGER  PROCEDURE INT3DBD(SHORT INTEGER BD,LOW,HIGHREAL ARRAY PLANE 
C00021 00006	SIMPLE PROCEDURE BESTPLANE(INTEGER IFRST,ILASTREAL ARRAY X,Y,Z,U,V,GPLANE,COORD)
C00029 00007	SIMPLE PROCEDURE LINEDPY(REAL ARRAY X,Y,Z,ANSINTEGER IFIRST,ILAST)
C00034 00008	PARABOLA
C00037 00009	INTEGER PROCEDURE AXIS (SHORT INTEGER  GRPREFERENCE INTEGER BD1,BD2,LOW1F,LOW2F,
C00051 00010	⊃ PART OF AXIS 
C00062 00011	PROCEDURE MK_PIECES
C00067 00012	⊃ MAIN PART OF OUTER 
C00068 00013	MAIN PROGRAM 
C00070 ENDMK
C⊗;
BEGIN "CURVE" 

COMMENT MODIFICATIONS OF CURV17.SAI,TERMINATION CONDITIONS PUT IN THE AXIS
	FINDER ;

REQUIRE "{}{}" DELIMITERS ;
DEFINE CRLF ={'15&'12},LF={'12},LF =  {'12},
	⊃={COMMENT},
	ITYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVS(FOO)&" ")END},
	RTYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVG(FOO)&" ")END},
	OTYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVOS(FOO)&" ")END},
	ATYPE(FOO)={BEGIN OUT(OUTFILE,"FOO"&" "&CVS(FOO)&" ")END},
	BTYPE(FOO)={BEGIN OUT(OUTFILE,"FOO"&" "&CVG(FOO)&" ")END},
	IOUT(CH,FOO)={BEGIN OUT(CH,CVS(FOO))END},
	ROUT(CH,FOO)={BEGIN OUT(CH,CVG(FOO))END},
	SIARRAY = {     INTEGER ARRAY},
	SRARRAY = {     REAL ARRAY}  ;
DEFINE MAXGROUPS = {25}, MAXBDS = {50} ;
INTEGER DISFINAL,DISPOINTS,DISCONT,FMODE,LMODE,
	FIRSTONLY,RUNBEFORE;
REAL PDEN;
INTEGER PPOG,LPOG,FINALPOG,
	COUNT,BRCHAR,EOF,FLAG,NPOINTS,NLINES,NNEW,NSEG,INDEX,WORDS,
	UPSIGN,RSIGN,SCANSIGN,DEBUG,DISP ;
DEFINE TTY={1},DSKIN={2},DSKOUT={3},PTFILE={4},PLOTFILE={5},OUTFILE= {6},BDFILE={7};
DEFINE P={'17};
DEFINE LOOP={←1 STEP 1 UNTIL};

REQUIRE -1 NEW_ITEMS ;

REQUIRE "DPYSUB.HDR[1,PDQ]" SOURCE_FILE;

COMMENT FOLLOWING DECLARATIONS HAVE COME FROM DPY3D.SAI[1,GJA];
DEFINE CALFILE = {1} , DSKIN = {4} ;
DEFINE PI = {3.1415926},BIG = {1.0@30} ;
REQUIRE "SAITRG[SYS,BGB]" SOURCE_FILE;
REQUIRE "MATRIX.HDR[L,RKN]" SOURCE_FILE ;

STRING FILNAM;
INTERNAL INTEGER MAXPOINTS,MAXLINES,MAXSEGS;
INTERNAL INTEGER ARRAY DPYBUF[1:1000];
INTERNAL INTEGER POG,WHICHFRAME;
INTERNAL REAL DPYMULT,XOFFSET,YOFFSET ;
INTERNAL REAL ARRAY SHIFT,UNSHIFIT,MA,MB[1:4,1:4];
INTERNAL REAL ARRAY ZEROPLANE[1:2,1:4] ;
INTERNAL REAL ARRAY LASLOC,LASAXIS[1:4] ;
INTERNAL INTEGER ARRAY ORIENT,SSIGN[1:2];
INTERNAL REAL ARRAY CAMTRANS[1:4,1:4],CAMVECT[1:8];
INTERNAL REAL STEPANGLE ;
INTERNAL REAL SEGTHRESH;


DEFINE NOTFOUND = {-10000.0},ORTHOGONAL = {-20000.0} ,
       FLAT = {-30000.0}, SPHERE = {-40000.0} ,OUTRANGE= {-50000.0} ;


⊃ INTEGER ARRAY PTS[1:MAXPOINTS]; ⊃ NOT USED ;
⊃ INTEGER ARRAY MID,BASSIZ[1:MAXSEGS]; ⊃ PTS AND BASSIZ NOT USED ANY MORE ;
⊃  PTS has coords,packed two to a word , OFFSET + 200.
  BASSIZ - Packed array.Left half contains BASE(accessed by a macro).
	   BASE points to first location in PTS for the segment.
	   Right half contains SIZSEG,the no of words of PTS taken up
	   by this segment.
  MID - Midpoint of segment, X in left half,Y in the right half.
	Stored values are just the sum of XBEG and YBEG.To get
	the right value,divide by two and subtract the offset 200.
  LINNO - line no of the segment.
BEGEND- Coordinates of the first and last point of the segment.
	Info pac
⊂!eft half contains beginning.First nine bits
	contain x coord,next nine y coord.OFFSET 200.
FLSEG-	Left half contains serial no of first segment in this line.
	Right half contains last segment.
GRPNO-	Group that this seg belongs to.-1 indicates does not belong to any.
	0 indicates seg not present.(larger than SEGS2).
BASGRP- Tells where the first segment in this group is located in GRPLST.
	(Gives index in GRPLST).
GRPLST- Linear list of segments.A group is contiguos no of members.No is 
	availaible in GRPSIZ.Starting point in BASGRP.Note some of the 
	entries may not belong to any group.
LINE1,SEGS1- Total no of the lines and segments in first scan.
LINE2,SEGS2- Total no of the lines and segments in both scans. ;

⊃ REQUIRE  "SOBMAT[SYS,HE]" LOAD_MODULE ;
⊃ EXTERNAL PROCEDURE QSOLVE (REAL ARRAY A,X,B);


REQUIRE "FILIN[1,RKN]" LOAD_MODULE ;
EXTERNAL SIMPLE PROCEDURE FILIN(REFERENCE STRING FILE);
EXTERNAL SIMPLE PROCEDURE GET_ARBDS; ⊃ PROCEDURES TO DEFINE MAX BOUNDS FOR PTFILE;

⊃ Following are the macros to unpack,packed info in arrays.;
	
DEFINE 

	BASE(ISEG) =  {(BASIZ[ISEG] LSH -18) } ,

	SIZSEG(ISEG) = { (BASIZ[ISEG] LAND '777777 ) },

	FIRSTSEG(LINE) = { (FLSEG[LINE] LSH - 18) } ,

	LASTSEG(LINE) = { (FLSEG[LINE] LAND '777777 ) } ,
	XBEG(SEG)= {( (BEGEND[SEG] LSH -27) - 200) },
	YBEG(SEG)= {(((BEGEND[SEG] LSH -18) LAND '777) - 200 )},
	XEND(SEG)= {(((BEGEND[SEG] LAND '777777) LSH -9)- 200) },
	YEND(SEG)= {( (BEGEND[SEG] LAND  '777) - 200 )} ,
	MIDX(SEG)= {( MID[SEG] LSH -18 )} ,
	MIDY(SEG)= {( MID[SEG] LAND '777777 ) } ,

	BDNOBEG(SEG) = { (BDNO[SEG] LSH -18 )} ,
	BDNOEND(SEG) = { (BDNO[SEG] LAND '777777 )} ,
	BDLOCBEG(SEG) = { (BDLOC[SEG] LSH -18 )} ,
	BDLOCEND(SEG) = { (BDLOC[SEG] LAND '777777 )} ,
	GET_BDNO(SEG) = { (IF SEG>0 THEN BDNOBEG(SEG) ELSE BDNOEND(-SEG)) },
	GET_BDLOC(SEG)= { (IF SEG>0 THEN BDLOCBEG(SEG) ELSE BDLOCEND(-SEG)) },
	BDSEG(LOC)= { ( (BDARRY[LOC] LAND '777777) - '100000 ) } ,
	BDMARK(LOC) = { (BDARRY[LOC] LSH -18) } ;
	
	
		
⊃ OUTER BLOCK STARTS HERE ;

PROCEDURE OUTER ;
BEGIN "OUTER" 

INTERNAL INTEGER ARRAY LINNO,BEGEND,GRPNO,GRPLST[1:MAXSEGS];
INTERNAL INTEGER ARRAY FLSEG,MIRANG[1:MAXLINES] ;
INTERNAL REAL ARRAY COLLIN[1:MAXLINES,1:4,1:3];
INTERNAL INTEGER ARRAY GRPSIZ,BASGRP[1:MAXGROUPS];
INTERNAL REAL ARRAY WIDTH[1:MAXSEGS];
INTERNAL INTEGER LINE1,LINE2,SEGS1,SEGS2;
INTERNAL INTEGER ARRAY BASBD,SIZBD[1:MAXBDS],BDARRY[1:2*MAXSEGS];
INTERNAL INTEGER ARRAY BDNO,BDLOC[1:MAXSEGS];
INTERNAL INTEGER NUMBDS; ⊃ Actual no of boundary segments found ;

REQUIRE "CUR19A[1,RKN]" LOAD_MODULE ;
EXTERNAL PROCEDURE READPTS(STRING FILE);
EXTERNAL SIMPLE PROCEDURE READCAL ;
EXTERNAL PROCEDURE BDIN ;
EXTERNAL SIMPLE PROCEDURE DPYSEG(INTEGER ISEG);
EXTERNAL PROCEDURE DISPLAY ;
EXTERNAL SIMPLE PROCEDURE GROUP ;
EXTERNAL SIMPLE PROCEDURE GETXYBD(INTEGER BD,INDX;REFERENCE REAL X,Y);
EXTERNAL SIMPLE PROCEDURE THREED(INTEGER LINE;REAL U,V;REFERENCE REAL X,Y,Z);
EXTERNAL SIMPLE PROCEDURE TWOD(REAL X,Y,Z;REFERENCE REAL U,V);
EXTERNAL PROCEDURE GET_BREAKS ;




SIMPLE INTEGER  PROCEDURE INT3DBD(SHORT INTEGER BD,LOW,HIGH;REAL ARRAY PLANE ;
			REFERENCE REAL XINT,YINT,ZINT;REFERENCE INTEGER BDINDX);
BEGIN "INT3DBD"
⊃ Finds intersection of given plane with the specified boundary segment.;
⊃ This one for 3d intersections ;

SHORT REAL MD,FD,LD,FX,FY,FZ,LX,LY,LZ,MX,MY,MZ,A,B,C,D,DIFF,UINT,VINT;
SHORT INTEGER  FU,FV,LU,LV,MU,MV,MID,DEBUG,BASB,FIRST,LAST,LINE,SEG1,SEG2,SIZB;
BOOLEAN INTDEB ;
LABEL L1,L2,L3;
DEFINE NEAR_LOW = {0}, NEAR_HIGH = {1}, FOUND = {-1} ;

⊃ Finds the coords of point in the current segment.You set BASSEG before
   calling it.;
SIMPLE  PROCEDURE COORD3(SHORT INTEGER INDX ; REFERENCE REAL X,Y,Z);
BEGIN "COORD3"
INTEGER SEG;
SEG← BDSEG(BASB+((INDX-1+SIZB) MOD SIZB)) ;
	⊃ The term INDX-1+SIZB is to account for 
	  a. negative INDX by adding SIZB to it.
	  b.-1 is so that INDX= SIZB does not give zero. ;
IF SEG> 0 THEN
	THREED(LINE← LINNO[SEG],XBEG(SEG),YBEG(SEG),X,Y,Z)
     ELSE 
	THREED(LINE← LINNO[-SEG],XEND(-SEG),YEND(-SEG),X,Y,Z)

END "COORD3";

IF HIGH= LOW THEN RETURN(0);
A← PLANE[1]; B←PLANE[2]; C← PLANE[3]; D← PLANE[4];
BASB← BASBD[BD];
SIZB← SIZBD[BD];

FIRST← LOW;LAST← HIGH;

COORD3(FIRST,FX,FY,FZ);
COORD3(LAST,LX,LY,LZ);

FD← A* FX + B*FY +C*FZ + D ;
LD← A* LX + B*LY +C*LZ + D ;

IF INTDEB THEN BEGIN
	RTYPE(FX);RTYPE(FY);RTYPE(FZ);RTYPE(LX);RTYPE(LY);RTYPE(LZ);
	RTYPE(FD);RTYPE(LD);OUTSTR(CRLF);END;

IF FD * LD ≤0 THEN BEGIN
	⊃ BINARY SEARCH WILL WORK HERE;
L1:	WHILE FIRST< LAST -1 DO BEGIN
		DIFF← LAST - FIRST ; ⊃ So that DIFF may be FSCed by compiler 
				rather than FLOATed. ;
		MID← FIRST + ABS(DIFF*FD/(FD - LD)) ;⊃ Note the minus sign 
		   in (FD - LD).Eliminates taking absolute values since 
		   FD and LD are both known to be of opposite sign.;
		IF MID = FIRST THEN MID←MID +1;
		IF MID = LAST THEN MID←LAST-1;
		COORD3(MID,MX,MY,MZ);
		MD← A* MX + B*MY + C*MZ + D ;
		IF FD * MD > 0 THEN BEGIN
			FIRST←MID;FD←MD;END
		ELSE BEGIN
			LAST←MID;LD←MD;END;
		IF INTDEB THEN BEGIN
			RTYPE(FD);RTYPE(LD);RTYPE(MD);OUTSTR(CRLF);
			RTYPE(FIRST);RTYPE(LAST);RTYPE(MID);OUTSTR(CRLF);
			END;		
	    END; ⊃ WHILE FIRST<LAST ;
	    END 
   ELSE BEGIN
      IF (HIGH-LOW)>6 THEN RETURN(IF ABS(FD)< ABS(LD) THEN NEAR_LOW ELSE NEAR_HIGH);
	⊃ Both ends on the same side of the line.No intersection.;
	⊃ SEARCH STEP BY STEP,IF LESS THEN FIVE POINTS ;
L2:	FOR LAST← LOW +1 STEP 1 UNTIL HIGH DO BEGIN
		COORD3(LAST,LX,LY,LZ);
		LD← A* LX + B* LY + C*LZ +D;
		IF FD*LD≤0 THEN DONE ELSE BEGIN
			IF LAST=HIGH THEN RETURN(0);
			⊃ ********* NOTE IT DOES NOT RETURN WHICH END IT'S NEAR ;
			FIRST← LAST; FD← LD ; 
			END;
		IF INTDEB THEN BEGIN
			RTYPE(FD);RTYPE(LD);OUTSTR(CRLF);
			RTYPE(FIRST);RTYPE(LAST);OUTSTR(CRLF);
			END;		
		END; ⊃ FOR LAST← ;
	END; ⊃ ELSE BEGIN ;

L3: COORD3(FIRST,FX,FY,FZ);
    COORD3(LAST,LX,LY,LZ);


FD←ABS(FD);LD←ABS(LD);
	

⊃ Find intersection by linear interpolation between the two points.;
DIFF← LX - FX ;
XINT←( FX+ (DIFF* FD)/(FD + LD) );
DIFF← LY - FY ;
YINT←( FY+ (DIFF* FD)/(FD + LD) );
DIFF← LZ - FZ ;
ZINT←( FZ+ (DIFF* FD)/(FD + LD) );
BDINDX← IF ABS(FD)<ABS(LD) THEN FIRST ELSE LAST ;

IF INTDEB THEN BEGIN
	RTYPE(XINT); RTYPE(YINT);RTYPE(ZINT);
	OUTSTR(CRLF);END;
 IF DISP THEN BEGIN
	TWOD(XINT,YINT,ZINT,UINT,VINT);
	APOINT(UINT*DPYMULT +XOFFSET,VINT*DPYMULT+YOFFSET);
	DPYOUT(1) ;
	END;
RETURN(FOUND);

END "INT3DBD";
SIMPLE PROCEDURE BESTPLANE(INTEGER IFRST,ILAST;REAL ARRAY X,Y,Z,U,V,GPLANE,COORD);
BEGIN "BESTPLANE"
COMMENT FIND THE BEST PLANE FITTING THE POINTS IN X,Y,Z[IFRST:ILAST]
    AND RELOCATE THOSE POINTS TO THE PLANE, PLACING THE RELOCATED POINTS IN
    U,V[IFRST:ILAST].
    SEE NOTEBOOK, VOLUME 11, PAGES 90-92;
INTEGER I,J,N,I1,I2,I3; REAL X1,Y1,Z1,D ;
OWN REAL ARRAY PSUMS[1:4,1:4];

MZERO(PSUMS);
FOR I←IFRST STEP 1 UNTIL ILAST DO BEGIN "BUILD PSUMS MATRIX"
    PSUMS[1,4]←PSUMS[1,4]+(X1←X[I]);
    PSUMS[2,4]←PSUMS[2,4]+(Y1←Y[I]);
    PSUMS[3,4]←PSUMS[3,4]+(Z1←Z[I]);
    PSUMS[1,1]←PSUMS[1,1]+X1*X1;
    PSUMS[1,2]←PSUMS[1,2]+X1*Y1;
    PSUMS[1,3]←PSUMS[1,3]+X1*Z1;
    PSUMS[2,2]←PSUMS[2,2]+Y1*Y1;
    PSUMS[2,3]←PSUMS[2,3]+Y1*Z1;
    PSUMS[3,3]←PSUMS[3,3]+Z1*Z1;
    END;
PSUMS[4,4]←N←ILAST-IFRST+1;
COMMENT TRANSLATE CENTER OF GRAVITY TO ORIGIN OF COORDINATES;
FOR I LOOP 3 DO BEGIN
    PSUMS[I,I]←PSUMS[I,I]-PSUMS[I,4]↑2/N;
    FOR J←I+1 STEP 1 UNTIL 3 DO
	PSUMS[J,I]←PSUMS[I,J]←PSUMS[I,J]-PSUMS[I,4]*PSUMS[J,4]/N;
    END;
COMMENT GUESS AT ORIENTATION FOR HOMOGENEOUS SOLUTION;
I1←IF PSUMS[1,1]<PSUMS[2,2] THEN 1 ELSE 2;
I1←IF PSUMS[I1,I1]<PSUMS[3,3] THEN I1 ELSE 3;
I2←I1 MOD 3 + 1;
I3←I2 MOD 3 + 1;
GPLANE[I1]←1.0;
COMMENT HOMOGENEOUS SOLUTION;
D←PSUMS[I2,I2]*PSUMS[I3,I3] - PSUMS[I2,I3]↑2;
GPLANE[I2]←-(PSUMS[I1,I2]*PSUMS[I3,I3] - PSUMS[I1,I3]*PSUMS[I2,I3])/D;
GPLANE[I3]←-(PSUMS[I1,I3]*PSUMS[I2,I2] - PSUMS[I1,I2]*PSUMS[I2,I3])/D;
COMMENT NORMALIZE;
D←SQRT(GPLANE[1]↑2 + GPLANE[2]↑2 + GPLANE[3]↑2);
FOR I LOOP 3 DO GPLANE[I]←GPLANE[I]/D;
GPLANE[4]←-(GPLANE[1]*PSUMS[1,4] + GPLANE[2]*PSUMS[2,4]
    +GPLANE[3]*PSUMS[3,4])/N;

COMMENT FORM COORDINATION MATRIX;
MCOORD(COORD,GPLANE);
COMMENT GET PROJECTION ON GPLANE ACCORDING TO INVERSE
    COORDINATION MATRIX;
FOR I←IFRST STEP 1 UNTIL ILAST DO BEGIN "PROJECT"
    U[I]←COORD[1,1]*X[I] COMMENT COORD[2,1] IS ZERO;
	+ COORD[3,1]*Z[I];
    V[I]←COORD[1,2]*X[I] + COORD[2,2]*Y[I]
	+ COORD[3,2]*Z[I];
    COMMENT THE SEGNO MATRIX WILL BE GARBAGE;
    END;
END "BESTPLANE";

SIMPLE REAL  PROCEDURE LINEFIT(INTEGER IFIRST,ILAST;REAL ARRAY U,V,ANS);
    BEGIN "LINEFIT"
COMMENT COMPUTE LOW ORDER SUMS
	PERFORM LEAST SQUARES  ;

REAL X,Y,Z,ERROR,SUMX,SUMXX,SUMY,SUMXY,SUMYY;
INTEGER I,J,DX,DY,BADPOINT,CONSEC,MAXOUT,N ;
REAL A,B,C,LAMBDA,TEMP,BADVAL,LAST;

N← ILAST-IFIRST + 1 ;


COMMENT COMPUTE LOW ORDER SUMS;
SUMX←SUMXX←SUMY←SUMXY←SUMYY←0;
FOR I←IFIRST STEP 1 UNTIL ILAST DO BEGIN
	SUMX←SUMX+(X←U[I]);
	SUMY←SUMY+(Y←V[I]);
	SUMXX←SUMXX+X*X;
	SUMYY←SUMYY+Y*Y;
	SUMXY←SUMXY+X*Y END;

COMMENT LEAST SQUARES COMPUTATION;
A←N*SUMXX-SUMX*SUMX;
B←N*SUMXY-SUMX*SUMY;
C←N*SUMYY-SUMY*SUMY;
IF B=0 THEN BEGIN
	IF A>C THEN X←0 ELSE X←1.0;
	Y←1.0-X END
    ELSE BEGIN
	LAMBDA←(A+C-((A-C)↑2+4*B*B)↑.5)/2;
	X←-B;
	Y←A-LAMBDA;
	A←(X↑2+Y↑2)↑.5;
	X←X/A;Y←Y/A END;
ANS[1]←X;
ANS[2]←Y;
ANS[3]←Z←-(SUMX*X+SUMY*Y)/N;

COMMENT OVERALL FIT;
A←SUMXX*X+SUMXY*Y+SUMX*Z;
B←SUMXY*X+SUMYY*Y+SUMY*Z;
C←SUMX*X+SUMY*Y+N*Z;
ERROR←((A*X+B*Y+C*Z)/N)↑.5;
RETURN(ERROR);

END "LINEFIT" ;


PROCEDURE LINE3D(INTEGER IFIRST,ILAST;REAL ARRAY X,Y,Z,ANS);
BEGIN "LINE3D" 

INTEGER I,J,K ;
REAL DEN ;
REAL ARRAY U,V[IFIRST:ILAST];
OWN REAL ARRAY GPLANE[1:4],COORD[1:4,1:4],XL,YL,ZL,UL,VL[1:2],ANS1[1:3];


BESTPLANE(IFIRST,ILAST,X,Y,Z,U,V,GPLANE,COORD);
LINEFIT(IFIRST,ILAST,U,V,ANS1);

⊃ NOW WE HAVE THE TWOD LINE,CONVERT TO THREE-D PARAMETERS ;

IF ABS(ANS1[1])>0.5 THEN BEGIN
	VL[1]← 1.0 ; VL[2]← 10.0 ;
	FOR I LOOP 2 DO
		UL[I]← -(ANS1[2]*VL[I] + ANS1[3])/ANS1[1] ;
	END
	ELSE BEGIN 
	   UL[1]← 1.0 ; UL[2]← 10.0 ;
  	   FOR I LOOP 2 DO
		VL[I]← -(ANS1[1]* UL[I] + ANS1[3])/ANS1[2] ;
	END;	

FOR I LOOP 2 DO BEGIN
	XL[I]← COORD[1,1]* UL[I] + COORD[1,2]* VL[I]  - COORD[1,3]*GPLANE[4];
					⊃ THE W COORDINATE FOR ALL POINTS IS THE DIST. 
					  OF THE PLANE FROM THE ORIGIN (- GPLANE[4]);
	YL[I]← COORD[2,1]* UL[I] + COORD[2,2]* VL[I]  - COORD[2,3]*GPLANE[4];
	ZL[I]← COORD[3,1]* UL[I] + COORD[3,2]* VL[I]  - COORD[3,3]*GPLANE[4];
	END;


ANS[1]← XL[2]- XL[1] ; ⊃ THE THREE DIRECTION COSINES ;
ANS[2]← YL[2]- YL[1] ;
ANS[3]← ZL[2]- ZL[1] ;

⊃ NORMALIZE THE COEFFICIENTS ;

DEN← SQRT(ANS[1]↑2 + ANS[2]↑2 + ANS[3]↑2 );


FOR I LOOP 3 DO ANS[I]← ANS[I]/DEN ;
ANS[4]← 1 ;
ANS[5]← XL[1] ; ⊃ THE CO-ORDINATES OF ONE POINT ON THE LINE ;
ANS[6]← YL[1] ;
ANS[7]← ZL[1] ;

ANS[8]← IFIRST ;ANS[9]← ILAST ; ⊃ STORE THE INDICES ;

END "LINE3D" ;

NE ;

I1← IF ABS(LINE[1])> ABS(LINE[2]) THEN 1 ELSE 2 ;
I1← IF ABS(LINE[I1])> ABS(LINE[3]) THEN I1 ELSE 3 ;

N1← LINE[8] ; N2← LINE[9]; ⊃ INDICES OF END POINTS ;

F[1]← X[N1];F[2]← Y[N1]; F[3]← Z[N1];
L[1]← X[N2];L[2]← Y[N2]; L[3]← Z[N2];

SIGN← IF (L[I1]-F[I1])> 0 THEN 1.0 ELSE -1.0 ;

PREVN← N - STEPSIZ ;
LP[1]← X[PREVN]; LP[2]← Y[PREVN]; LP[3]← Z[PREVN] ; ⊃ COORDINATES OF THE LAST POINT ;
INCR← ABS(STEPMULT * RAD[PREVN] *STEPSIZ) MAX MINSTEP ;
IF STEPSIZ<0 THEN INCR← -INCR ;

XYZ[I1]← LP[I1] + INCR * SIGN ;


I2← I1 ;
FOR I LOOP 2 DO BEGIN
	I2← I2 MOD 3 + 1 ;
	XYZ[I2]← LINE[I2+4] + LINE[I2]*(XYZ[I1] - LINE[I1+4])/LINE[I1] ;
	END;

XNP← XYZ[1] ;YNP← XYZ[2]; ZNP← XYZ[3] ;


END "LINEEXP" ;



COMMENT PARABOLA;
REAL SIMPLE PROCEDURE APARABOLA(SHORT INTEGER INDXL,INDXH;REAL ARRAY RAD;
			REFERENCE REAL A,B,C;REFERENCE INTEGER OFFSET);
BEGIN "APARABOLA"
COMMENT FINDS THE LEAST SQUARES FIT OF THE WIDTHS OF SEGMENTS IN
GROUP TO THE EQUATION    WIDTH = A*(X-OFFSET)↑2 + B*(X-OFFSET) + C
	SEE NOTEBOOK, VOL 11, PAGES 89-90;

REAL N,N2,N3,N4,N5,N6,N7,N8,N9,D,Y,XY,XXY,YY,E,F,G;
INTEGER I,J;

OFFSET← INDXL -1 ;

COMMENT GET N MATRIX;
N←INDXH - INDXL +1 ;

IF N≤2 THEN BEGIN
	⊃ Not enough points to fit a parabola,so fit a straight line ;
	A← 0 ; C← RAD[INDXL];
	B← (RAD[INDXH] - C);
	RETURN(0); ⊃ No error ;
	END;


N2←N*N; N3←N2*N; N4←N3*N; N5←N4*N; N6←N5*N; N7←N6*N; N8←N7*N; N9←N8*N;

COMMENT GET THE Y MATRIX;
Y←XY←XXY←YY←0; J←0;
FOR I ← INDXL STEP 1 UNTIL INDXH DO BEGIN
	D←RAD[I];
	J←J+1;
	Y←Y+D;
	XY←XY+D*J;
	XXY←XXY+D*J*J;
	YY←YY+D*D;
	END;

COMMENT SOLVE THE EQUATIONS;
D ← N9 - 6*N7 + 9*N5 - 4*N3;
A ← (XXY * 180*(N4 - N2)
    + XY * (E ← 180*(-N5 - N4 + N3 + N2))
    +  Y * (F ← 30*(N6 + 3*N5 + N4 - 3*N3 - 2*N2)))/D;
B ← (XXY * E
    + XY * 12*(16*N6 + 30*N5 - 5*N4 - 30*N3 - 11*N2)
    +  Y * (G ← 18*(-2*N7 - 7*N6 - 5*N5 + 5*N4 + 7*N3 + 2*N2)))/D;
C ← (XXY * F
    + XY * G
    +  Y * 3*(3*N8 + 12*N7 + 14*N6 - 13*N4 - 12*N3 - 4*N2) )/D;

COMMENT COMPUTE THE ERROR;
RETURN (( A↑2 * (6*N5 + 15*N4 + 10*N3 - N) / 30
	+ 2 * A * B * (N4 + 2*N3 + N2) / 4
	+ (B↑2 + 2*A*C) * (2*N3 + 3*N2 + N) / 6
	+ 2 * B * C * (N2 + N) / 2
	+ C↑2 * N
	- 2*A*XXY - 2*B*XY - 2*C*Y + YY) / N);
END "APARABOLA";
INTEGER PROCEDURE AXIS (SHORT INTEGER  GRP;REFERENCE INTEGER BD1,BD2,LOW1F,LOW2F,
        HIGH1F,HIGH2F,BSIGN,PARB,PARE,TERMB,TERME,CIRC;REAL ARRAY AXX,AXY,AXZ,RAD,ANGLE);
BEGIN "AXIS" 
SHORT REAL X1,Y1,Z1,X2,Y2,Z2,XN,YN,ZN,XN1,YN1,ZN1,XN2,YN2,ZN2,COR,XCOR,YCOR,ZCOR,
	   DX,DY,U1,V1,U2,V2,CORTHRESH ,ERR,XNP,YNP,ZNP,A,B,C,RADERR,RADTHRESH;
OWN REAL ARRAY PLANE[1:4],F,L[1:3],AXLINE,INITLIN[1:9] ;
SHORT INTEGER   I,J,K,SIZG,BASG,MIDSEG,N1,N2,N11,N22,NM,SEG1,SEG2,SEGM,ISEG,LINE,
		SIZE1,SIZE2,FITSIZ ,M,I1,PAR,PAR1,PAR2,IPAR,PARBEG,PAREND,
		STEPSIZ,INITNO ,IFIRST,ILAST,I2,I3,BSIGN1,BSIGN2,LEN1,LEN2,OFFSET,
		LOW1,HIGH1,LOW2,HIGH2,MID1,MID2 ,
		LOW11,HIGH11,LOW22,HIGH22,LOW10,LOW20,HIGH10,HIGH20,INDX1,INDX2,
		HSIZE1,HSIZE2,TEMP,SIGN,FINDX,LINDX,LEN ;

BOOLEAN AXDEB,TTEST ;
LABEL L1;
DEFINE TOP= "1" , BOTTOM = "2" ;
DEFINE NEAR_LOW = {0},NEAR_HIGH = {1}, FOUND = {-1} ;
DEFINE EXP_RADIUS(PAR) = {(A*(PAR-OFFSET)↑2 + B*(PAR-OF\X1ε'FSET) + C) };
DEFINE BDDIST(F,L,SIZE)= {(IF (TEMP←ABS(F-L) MOD SIZE)<(SIZE/2)
				THEN (TEMP-1) ELSE (SIZE-TEMP-1))} ;


REAL SIMPLE PROCEDURE DIST_LINE(INTEGER BD,INDX;REAL ARRAY LINE);
BEGIN "DIST_LINE"
SHORT REAL X,Y;
GETXYBD(BD,INDX,X,Y);
RETURN(LINE[1]*X + LINE[2]*Y + LINE[3]);
END "DIST_LINE";


SIMPLE BOOLEAN PROCEDURE TERMTEST(INTEGER BD,INDX1,INDX2,SIZE,SIGN;REAL THRESH;
					REFERENCE INTEGER LEN;REAL ARRAY LINE);
BEGIN "TERMTEST"
	⊃ TEST WHETHER BOUNDARY EXTNDS FAR ENOUGH,BETWEEN INDX1 AND INDX2(MODULO
	  SIZE) IN THE DIRECTION SIGN;
SHORT INTEGER I,SSIZE,NSTEPS,INDX ;
SHORT REAL DIST ;

LEN← ABS(INDX1 - INDX2) MOD SIZE ;
IF LEN> (SIZE/2) THEN LEN← SIZE - LEN ;
IF LEN≤1 THEN RETURN(TRUE); ⊃ OBVIOUSLY NEAR BOUNDARY;
SSIZE← SIGN * (5 MIN (LEN1/5)); ⊃ CHOOSE AN APPROPRIATE SIZE ;
NSTEPS← LEN1/SSIZE ;

INDX←INDX1;
FOR I LOOP NSTEPS DO BEGIN
	DIST ← DIST_LINE(BD,INDX,LINE);
	IF ABS(DIST)>(5.0 MAX THRESH) THEN RETURN(FALSE);
	INDX←INDX + NSTEPS ;
	END;
RETURN(TRUE);
END "TERMTEST";

BOOLEAN PROCEDURE RADCHECK(INTEGER PAR,PAR1,PAR2,STEPSIZ;REAL ARRAY RAD);
BEGIN "RADCHECK"
	⊃ CHECK TO SEE IF THE PRESENT RADIUS FUNCTION IS ACCEPTABLE;
REAL A,B,C,PRAD,LRAD,RADERR,RADTHRESH,DIFF,ERAD ;
SHORT INTEGER F,L,OFFSET;
DIFF← ABS( (PRAD←RAD[PAR])-(LRAD←RAD[PAR-STEPSIZ]));
IF DIFF< (0.15*LRAD MAX 2.0)  THEN RETURN(TRUE) ; ⊃ ITS OK,KEEP IT;
IF DIFF> (LRAD MAX 4.0) THEN RETURN(FALSE);⊃ CHANGE IS TOO LARGE,NEED NOT CHECK FURTHER;

⊃ NOW CHECK FOR MORE GLOBAL EVALUATIONS;
  ⊃ FIRST FIT A PARABOLA TO THE PREVIOUS POINT ;
  IF STEPSIZ>0 THEN BEGIN F←PAR1; L←PAR2-STEPSIZ ;END 
	       ELSE BEGIN F←PAR1-STEPSIZ ;L←PAR2; END;
  APARABOLA(F,L,RAD,A,B,C,OFFSET);
  ERAD← EXP_RADIUS(PAR);
RETURN(IF ABS(RADERR←(ERAD- PRAD))> (RADTHRESH← (0.25*ERAD MAX 5.0)) THEN
	   FALSE ELSE TRUE);

END "RADCHECK";

SIMPLE BOOLEAN PROCEDURE ISTERM(SHORT INTEGER BD,INDX1,INDX2,F,L,SIGN);
BEGIN "ISTERM"
⊃ DETERMINE WHETHER THE SEGMENT FROM F TO L OF THE BD CAN BE CONSIDERED AS
  TERMINATION,INDX1,INDX2 THE INDICES OF END-POINTS;
SHORT INTEGER INDX,SIZB ;
SHORT REAL X1,Y1,X2,Y2,X,Y,XN,YN,DIST,T1,T2,T3,T4,T5,T6,T7,THETA,XT,YT;

DEFINE ATHRESH = {10.0}, XTHRESH = {7.0} ; ⊃ THRESHHOLDS ALONG AND NORMAL TO AXIS;
DEFINE EDDIST(INDX) =
		{ BEGIN
		  GETXYBD(BD,INDX,X,Y);
		 XN← T1* (XT← X- T5) + T2 * (YT← Y - T6);
		 YN← ABS(T3* XT + T4 * YT) ;
		 XN← IF XN< 0 THEN ABS(XN) ELSE IF XN≤ T7 THEN 0.0 ELSE XN - T7 ;
			⊃ SEE WHETHER BETWEEN THE SEGMENT OR DIST FROM THE NEAR END;
		END; } ;

SIZB← SIZBD[BD];
⊃ FIRST FORM A TRANSFORMATION TO ALIGN COORDINATES WITH EDGE JOINING INDX1
 AND INDX2 ;
GETXYBD(BD,INDX1,X1,Y1);
GETXYBD(BD,INDX2,X2,Y2);
THETA← ATAN2(Y2-Y1,X2-X1);
T1← T4← COS(THETA);
T3← -(T2← SIN(THETA)); ⊃ COEFFICIENTS OF THE TRANSFORMATION;
T5← X1; T6← Y1; ⊃ COORDINATES OF THE NEW ORIGIN;
T7← SQRT((X1- X2)↑2 + (Y1-Y2)↑2); ⊃ LENGTH OF THE EDGE ;

INDX← F;
WHILE TRUE DO BEGIN
	EDDIST(INDX,DIST); ⊃ THE DISTANCE FROM EDGE  ;
	IF XN > XTHRESH  OR YN> ATHRESH THEN RETURN(FALSE);
	INDX← INDX + SIGN ;
	IF (ABS(INDX- L) MOD SIZB ) = 0 THEN DONE ;⊃ L=INDX MODULO SIZB ;
	END;
RETURN(TRUE); ⊃ NONE OF THE DIST WAS LARGER THEN THRESH ;
END "ISTERM";


DPYSET(DPYBUF);
SIZG←GRPSIZ[GRP];
BASG←BASGRP[GRP];
N11←1; N22← SIZG;

⊃ Find parameters to relate to the boundaries ;
SEG1← GRPLST[BASG+N11-1];
SEG2← GRPLST[BASG+N22-1];
⊃ First find the boundary corresponding to SEG1 ;
WHILE TRUE DO BEGIN
	IF N22-N11 <3 THEN BEGIN 
		OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
		RETURN(0); END;
	BD1← BDNOBEG(SEG1);BD2← BDNOEND(SEG1);⊃ Pick the opposite boundaries;
	IF BD1=0 OR BD2=0 THEN BEGIN
		N11←N11+1;
		SEG1← GRPLST[BASG+N11-1];END
		ELSE DONE;
	END;

SIZE1← SIZBD[BD1];SIZE2← SIZBD[BD2];
HSIZE1← SIZE1/2; HSIZE2← SIZE2/2; ⊃ Half the size ;
LOW1← BDLOCBEG(SEG1); LOW2← BDLOCEND(SEG1);

⊃ FIND THE BOUNDARY SEGMENT CORRESPONDING TO SEG2.THIS MAY REQUIRE A SEARCH;
WHILE TRUE DO BEGIN
	IF N22-N11 <3 THEN BEGIN 
		OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
		RETURN(0); END;
	HIGH1← BDLOCBEG(SEG2); HIGH2← BDLOCEND(SEG2);
	IF HIGH1= 0 OR HIGH2=0 THEN BEGIN
		N22← N22 -1;
		SEG2← GRPLST[BASG+N22-1];END 
		ELSE DONE;
	END;
⊃ NOW FIND A SEGMENT IN THE MIDDLE ;
NM← N11+1;
SEGM←GRPLST[BASG+NM-1];
WHILE TRUE DO BEGIN
	IF N22-NM <1 THEN BEGIN 
		OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
		RETURN(0); END;
	MID1← BDLOCBEG(SEGM); MID2← BDLOCEND(SEGM);
	IF MID1= 0 OR MID2=0 THEN BEGIN
		NM← NM +1;
		SEGM← GRPLST[BASG+NM-1];END 
		ELSE DONE;
	END;
IF (HIGH1-MID1)*(MID1-LOW1)<0 THEN
	IF LOW1≤HIGH1 THEN LOW1← LOW1 +SIZE1
		      ELSE HIGH1←HIGH1 + SIZE1 ; ⊃ We have just gone around;
IF (HIGH2-MID2)*(MID2-LOW2)<0 THEN 
	IF LOW2≤HIGH2 THEN LOW2← LOW2 +SIZE2
		      ELSE HIGH2←HIGH2 + SIZE2 ; ⊃ We have just gone around;
IF LOW1> HIGH1 THEN BEGIN LOW1↔HIGH1;BSIGN1← -1; END ELSE BSIGN1← 1 ;
IF LOW2> HIGH2 THEN BEGIN LOW2↔HIGH2;BSIGN2← -1; END ELSE BSIGN2← 1 ;
BSIGN← IF BSIGN1= BSIGN2 THEN 1 ELSE -1;


⊃ PICK THE MIDDLE SECTION TO WORK ON ;
IF SIZG> 16 THEN BEGIN
	N1← SIZG /4 ;
	N2← 3*SIZG/4 ;END
	ELSE BEGIN N1←1; N2← SIZG;END;
FITSIZ← N2 - N1 +1 ;
PAR1← 1; PAR2← FITSIZ ;

FOR I← 1 STEP 1 UNTIL FITSIZ DO BEGIN
	ISEG← IF BSIGN1=1 THEN GRPLST[BASG+I+N1-2] ELSE GRPLST[BASG+N2-I];
	⊃ the axis always goes from LOW1 to HIGH1 ;
	U1←XBEG(ISEG);V1←YBEG(ISEG);
	THREED(LINNO[ISEG],U1,V1,X1,Y1,Z1);
	U2←XEND(ISEG);V2←YEND(ISEG); 
	THREED(LINNO[ISEG],U2,V2,X2,Y2,Z2);
	AXX[I]← (X1+X2)/2;AXY[I]← (Y1+Y2)/2;AXZ[I]← (Z1+Z2)/2;
	DX← ((U1+U2)/2)*DPYMULT + XOFFSET ;
  	DY← ((V1+V2)/2)*DPYMULT + XOFFSET ;
	APOINT(DX,DY);
	END;
DPYOUT(1);
⊃ NOW FIT  A LINE TO THIS ;

LINE3D(PAR1,PAR2,AXX,AXY,AXZ,AXLINE);
LINEDPY(AXX,AXY,AXZ,AXLINE,PAR1,FITSIZ);

⊃ NOW ITERATE ON THIS PART ;
⊃ LET US START WORKING WITH THE BOUNDARY RIGHT FROM HERE ;
	

CORTHRESH ← 0.1 ;
FOR J LOOP  5 DO BEGIN "INIT_ITER" 
INTEGER I2,I3,INTS1,INTS2;
	DEFINE ADD1 = {
		⊃ NOW ADD FIVE MORE SEGMENTS ;
		IF INTS1= NEAR_HIGH THEN HIGH1← HIGH1 + 5 ELSE LOW1← LOW1-5;
		IF ( (BD1= BD2) AND ( (ABS((HIGH1-LOW2) MOD SIZE1)<5) OR
			 	      (ABS((HIGH2-LOW1) MOD SIZE1)<5) ))
		   OR (ABS((HIGH1-LOW1) MOD SIZE1)< 5)  THEN DONE; };
	DEFINE ADD2 = {
		⊃ NOW ADD FIVE MORE SEGMENTS ;
		IF INTS2= NEAR_HIGH THEN HIGH2← HIGH2+5 ELSE LOW2← LOW2-5;
		IF ( (BD1= BD2) AND ( (ABS((HIGH2-LOW1) MOD SIZE2)<5) OR
			 	      (ABS((HIGH1-LOW2) MOD SIZE2)<5) ))
		   OR (ABS((HIGH2-LOW2) MOD SIZE2)< 5)  THEN DONE; };
			
	COR← 0 ;M← 0 ;
	DPYSET(DPYBUF);
	FOR I← PAR1 STEP 1 UNTIL PAR2 DO BEGIN
		NORMAL(AXLINE,PLANE,AXX[I],AXY[I],AXZ[I]);
		I2← 0;
	        IF (INT3DBD(BD1,LOW1,HIGH1,PLANE,XN1,YN1,ZN1,INDX1)= FOUND ) AND 
	           (INT3DBD(BD2,LOW2,HIGH2,PLANE,XN2,YN2,ZN2,INDX2)= FOUND ) THEN BEGIN
			⊃ We have found intersections with both boundaries ;
			M← M +1 ;
			IF M= 1 THEN BEGIN PAR1← I ; LOW1F←HIGH1F←INDX1;
					   LOW2F←HIGH2F← INDX2; END 
			    ELSE BEGIN
				HIGH1F←INDX1;
				IF BSIGN>0 THEN HIGH2F←INDX2 ELSE LOW2F←INDX2;END;
			PAR← PAR1 + M-1 ;

			XN← (XN1+XN2)/2; YN ←(YN1+YN2)/2 ; ZN← (ZN1+ ZN2)/2;
		 	XCOR← XN - AXX[I];YCOR← YN - AXY[I];ZCOR← ZN -AXZ[I];
			IF AXDEB THEN BEGIN
 			        ITYPE(I);ITYPE(M);RTYPE(XCOR);RTYPE(YCOR);RTYPE(ZCOR);
				OUTSTR(CRLF); END;
			COR ← COR + SQRT(XCOR↑2 + YCOR↑2 + ZCOR↑2 ) ;
			AXX[PAR]← XN ;AXY[PAR]← YN ; AXZ[PAR]← ZN ;
			RAD[PAR]←SQRT((XN1-XN2)↑2 + (YN1-YN2)↑2 + (ZN1-ZN2)↑2) ;
⊃ *************		ANGLE[PAR]← AXLINE[6];
			TWOD(XN,YN,ZN,DX,DY);
			APOINT(DX*DPYMULT +XOFFSET,DY*DPYMULT + YOFFSET);
			END;
		END;
	DPYOUT(1);
 	PAR2←  PAR ;
	IF M ≤ 1 THEN BEGIN OUTSTR("NOT ENOUGH POINTS LEFT"&CRLF);
			    RETURN(0); END;
	LINE3D(PAR1,PAR2,AXX,AXY,AXZ,AXLINE);
	LINEDPY(AXX,AXY,AXZ,AXLINE,PAR1,PAR2);
	IF (COR/M) < CORTHRESH THEN DONE "INIT_ITER" ;
END "INIT_ITER" ;
⊃ PART OF AXIS ;

ARRTRAN(INITLIN,AXLINE); ⊃ SAVE THE INITIAL LINE FOUND FOR LATER USE ;
INITNO ← M ;
⊃ NOW EXTEND THIS AXIS ;
⊃ VARIABLES PAR1 AND PAR2 STAND FOR THE LOWER AND UPPER INDICES OF SEGMENTS 
  ALREADY INCLUDED.PAR IS THE INDEX OF THE CURRENT SEGMENT.;

APARABOLA(PAR1,PAR2,RAD,A,B,C,OFFSET);
⊃ CONVENTION FOR BOUNDARY VARIBLES .
  LOW10,LOW20 etc. for storing results of initial iterations.
  LOW1,HIGH1 etc. for keeping track of how far the boundaries have been extended,
		  to avoid crossing over .
  LOW11,LOW22 etc. for picking a segment of boundary to pass along for finding
		intersections with(i.e. just the current extensions).
  LOW1F,LOW2F etc. for storing the `finished' results ;

STEPSIZ ← 1 ; PAR← PAR2; ⊃ Initialize for loop EXTND ;
LOW1←LOW1F;HIGH1←HIGH1F;LOW10←LOW1F+1;HIGH10←HIGH1F-1;⊃ START JUST INSIDE
								OF INITIAL GROUP ;
LOW2←LOW2F;HIGH2←HIGH2F;LOW20←LOW2F+1;HIGH20←HIGH2F-1;
LOW11←HIGH11←HIGH10;
IF BSIGN>0 THEN LOW22←HIGH22←HIGH20 ELSE LOW22←HIGH22←LOW20 ;
TERMB←TERME← 0 ; ⊃ More Initialization ;

WHILE TRUE DO BEGIN  "EXTND"
 INTEGER SIGN,SEG,SEGF,SEGL,NSEG,NDX1,NDX2;
 BOOLEAN NEWBD1,NEWBD2;
 DEFINE ADDSEG1 = {
	⊃ NOW ADD FIVE MORE SEGMENTS ;
	IF STEPSIZ> 0 THEN BEGIN
	   IF BD1= BD2 THEN 
		BEGIN IF (LEN1← BDDIST(HIGH1,LOW2,SIZE1))≤0 THEN
				BEGIN ⊃ TERME←-1; DONE "ONEDIR" ;END; END
	        ELSE IF (LEN1← BDDIST(HIGH1,LOW1,SIZE1))≤ 0 THEN
				BEGIN ⊃ CIRC←-1 ;DONE "ONEDIR";END;
	   LOW11← HIGH11 ;
	   HIGH11← HIGH1← HIGH1 + ( 5 MIN LEN1 );
	   END
	 ELSE BEGIN
	   IF BD1= BD2 THEN 
		BEGIN IF (LEN1← BDDIST(HIGH2,LOW1,SIZE1))≤0 THEN
				BEGIN ⊃ TERMB←-1;DONE "ONEDIR" ;END; END
		ELSE IF (LEN1← BDDIST(HIGH1,LOW1,SIZE1))≤ 0  THEN 
				BEGIN ⊃ CIRC←-1 ;DONE "ONEDIR";END;
	   HIGH11← LOW11 ;
	   LOW11← LOW1← LOW1-(5 MIN LEN1);
	   END;  };
 DEFINE ADDSEG2 = {
	⊃ NOW ADD FIVE MORE SEGMENTS ;
	IF STEPSIZ* BSIGN >0 THEN BEGIN
	   IF BD1= BD2 THEN 
		BEGIN IF (LEN2← BDDIST(HIGH2,LOW1,SIZE2))≤0 THEN BEGIN
		 ⊃  IF STEPSIZ>0 THEN TERME←-1 ELSE TERMB←-1 ;DONE "ONEDIR";END; END
		ELSE IF (LEN2← BDDIST(HIGH2,LOW2,SIZE2))≤ 0  THEN BEGIN
				⊃  CIRC ← ¬1 ;DONE "ONEDIR";END;
	   LOW22← HIGH22 ;
	   HIGH22← HIGH2← HIGH2 + ( 5 MIN LEN2); 
	   END
	 ELSE BEGIN
	   IF BD1= BD2 THEN 
		BEGIN IF (LEN2← BDDIST(HIGH1,LOW2,SIZE2))≤0 THEN BEGIN
		  ⊃ IF STEPSIZ>0 THEN TERME←-1 ELSE TERMB←-1 ;DONE "ONEDIR";END; END
		ELSE IF (LEN2← BDDIST(HIGH2,LOW2,SIZE2))≤ 0  THEN BEGIN
		 		⊃ CIRC← -1;DONE "ONEDIR";END;
	   HIGH22← LOW22 ;
	   LOW22← LOW2← LOW2 -( 5 MIN LEN2);
	   END;  };
    WHILE TRUE DO BEGIN "ONEDIR"		
	PAR← PAR + STEPSIZ ;
	IF STEPSIZ> 0 THEN BEGIN PAR2← ILAST← PAR;
				IFIRST←(PAR1 MAX (PAR-4));⊃ FIVE POINTS IF AVAILAIBLE;
				END
		      ELSE BEGIN PAR1←IFIRST← PAR;
				 ILAST← (PAR2 MIN (PAR+4));END;
	M← M+1 ;
	FOR I LOOP 1 DO BEGIN "ONESTEP"
	    ⊃ The loop is set up to iterate only once ;
	    LINEEXP(AXLINE,AXX,AXY,AXZ,RAD,PAR,STEPSIZ,XNP,YNP,ZNP);
            NORMAL(AXLINE,PLANE,XNP,YNP,ZNP);
	    I2← 0;
	    WHILE INT3DBD(BD1,LOW11,HIGH11,PLANE,XN1,YN1,ZN1,INDX1)≠FOUND DO BEGIN
		⊃ PUT MORE SEGMENTS IN THIS GROUP ;
		IF I2=3 THEN DONE "ONEDIR"; ⊃ You have tried to extend three times;
		ADDSEG1;I2← I2+1 ; END;
	    ⊃ RESET THE INDICES ;
	    IF STEPSIZ>0 THEN BEGIN LOW11←HIGH11←INDX1-1;HIGH1←INDX1;END 
			 ELSE BEGIN LOW11←HIGH11←INDX1+1;LOW1←INDX1;END;
	    I3←0;
	    WHILE INT3DBD(BD2,LOW22,HIGH22,PLANE,XN2,YN2,ZN2,INDX2)≠FOUND DO BEGIN
		IF I3=3 THEN DONE "ONEDIR"; ⊃ You have tried to extend three times;
		ADDSEG2;I3← I3+1 ;END;
		XN←(XN1+XN2)/2;YN←(YN1+YN2)/2; ZN← (ZN1+ZN2)/2 ;
		AXX[PAR]← XN ;AXY[PAR]← YN ;AXZ[PAR]← ZN ;
		RAD[PAR]← SQRT((XN1-XN2)↑2 + (YN1-YN2)↑2 +(ZN1-ZN2)↑2 );
		TWOD(XN,YN,ZN,DX,DY);
		APOINT(DX*DPYMULT+XOFFSET,DY*DPYMULT + YOFFSET);
		DPYOUT(1);
		⊃  IF ABS(RADERR← (EXP_RADIUS(PAR)- RAD[PAR]))> (RADTHRESH ←0.5*C)
						THEN DONE "ONEDIR" ;
⊃ ***********	IF NOT RADCHECK(PAR,PAR1,PAR2,STEPSIZ,RAD) THEN DONE "ONEDIR";
⊃ ***********	ANGLE[PAR]← AXLINE[6];
		ERR← SQRT( (XN-XNP)↑2 + (YN-YNP)↑2 +(ZN-ZNP)↑2 );

		IF STEPSIZ>0 THEN BEGIN
		   HIGH1F←INDX1;
		   IF BSIGN>0 THEN BEGIN LOW22←HIGH22←INDX2-1;HIGH2←HIGH2F←INDX2;END 
			      ELSE BEGIN LOW22←HIGH22←INDX2+1;LOW2←LOW2F←INDX2; END;
		   END
		 ELSE BEGIN
		   LOW1F←INDX1;
		   IF BSIGN>0 THEN BEGIN LOW22←HIGH22←INDX2+1;LOW2←LOW2F←INDX2;END
			      ELSE BEGIN LOW22←HIGH22←INDX2-1;HIGH2←HIGH2F←INDX2; END;
		 END;


		⊃ CHECK HERE WHETHER ANY FARTHER TO GO IN THIS DIRECTION ;
		IF STEPSIZ>0 THEN  BEGIN
		     FINDX← HIGH1F;
			IF BD1=BD2 THEN  BEGIN LINDX← LOW1F;SIGN←1 ;END
				ELSE BEGIN LINDX←IF BSIGN>0 THEN HIGH2F ELSE LOW2F;	
						SIGN← BSIGN ; END; END
		ELSE BEGIN
		     FINDX← LOW1F;
			IF BD1=BD2 THEN  BEGIN LINDX← HIGH1F;SIGN←-1 ;END
				ELSE BEGIN LINDX←IF BSIGN>0 THEN LOW2F ELSE HIGH2F;	
						SIGN← - BSIGN ; END; END ;
		IF TTEST THEN 
		  IF TERMTEST(BD1,FINDX,LINDX,SIZE1,SIGN,RAD[PAR],LEN,PLANE)
					THEN  DONE "ONEDIR";
			
			IF ERR < 0.2 THEN DONE "ONESTEP" ELSE BEGIN 
			⊃ FIT A NEW CURVE ;
			LINE3D(IFIRST,ILAST,AXX,AXY,AXZ,AXLINE);
			LINEDPY(AXX,AXY,AXZ,AXLINE,IFIRST,ILAST);
			⊃   APARABOLA(PAR1,PAR2,RAD,A,B,C,OFFSET);
			END  ;
	END "ONESTEP" ;
     END "ONEDIR";
     L1:IF STEPSIZ> 0 THEN BEGIN
		⊃ RESTORE THE INDICES TO HOW FAR WE HAD ACTUALLY GONE;
		HIGH1← HIGH1F;
		IF BSIGN>0 THEN HIGH2← HIGH2F ELSE LOW2← LOW2F;
		⊃  CHECK TO SEE IF YOU WERE NEAR THE TERMINATION;
		IF BD1=BD2 THEN
		      BEGIN IF ISTERM(BD1,HIGH1F,LOW2F,HIGH1F,LOW2F,STEPSIZ)
								THEN TERME←-1;END;
		⊃ CHECK FOR CIRCULARITY WILL BE MADE ONLY WHEN GOING BACKWARDS;
		STEPSIZ← -STEPSIZ;
		⊃ NOW RESTORE;
		ARRTRAN(AXLINE,INITLIN);PAR2← PAR2 -1;M← M-1;
		LOW11←HIGH11←LOW10;
		IF BSIGN>0 THEN LOW22←HIGH22←LOW20 ELSE LOW22←HIGH22←HIGH20;
		PAR← PAR1; END
	   ELSE BEGIN
		IF BD1=BD2 THEN
		      BEGIN IF ISTERM(BD1,LOW1F,HIGH2F,LOW1F,HIGH2F,STEPSIZ)
							THEN TERMB←-1;END
		      ELSE BEGIN 
			IF BSIGN>0 THEN BEGIN NDX1←LOW2F; NDX2← HIGH2F; END
				   ELSE BEGIN NDX1← HIGH2F;NDX2← LOW2F; END;
		        IF ISTERM(BD1,LOW1F,NDX1,LOW1F,HIGH1F,STEPSIZ) AND
			   ISTERM(BD2,LOW1F,NDX1,NDX1,NDX2,BSIGN*STEPSIZ) THEN
							CIRC← -1; END;
		PAR1← PAR1 +1;M← M-1 ;
		DONE "EXTND" ;END;

END  "EXTND" ;

⊃ NOW DISPLAY ALL THE AXIS POINTS ;
  DPYSET(DPYBUF);
  FOR IPAR← PAR1 STEP 1 UNTIL PAR2 DO BEGIN
	TWOD(AXX[IPAR],AXY[IPAR],AXZ[IPAR],DX,DY);
	APOINT(DX*DPYMULT + XOFFSET,DY*DPYMULT + YOFFSET); END;
  DPYOUT(1);	

PARB←PAR1; PARE← PAR2; ⊃ BEGIN AND AND INDICES ;
RETURN(PAR2-PAR1+1); ⊃ NO OF AXIS POINTS ;

END "AXIS" ;

PROCEDURE MK_PIECES;
BEGIN "MK_PIECES"
SHORT INTEGER BD1,BD2,LOW1F,LOW2F,HIGH1F,HIGH2F,BSIGN,PARB,PARE,I,J,
	      BASPC,SIZPC,IGRP,NPIECES,TERMB,TERME,CIRC,SIZ1,SIZ2;
REAL ARRAY AXX,AXY,AXZ,RAD,ANGLE[-100:100];
INTEGER ARRAY FIXUP[1:128];
LABEL L1,L2;
STRING FILNAM,STR ;

SIMPLE PROCEDURE FILEOUT;
BEGIN "FILEOUT"
INTEGER DUMMY ;
LABEL AGAIN;
OPEN(DSKOUT,"DSK",8,2,2,COUNT,BRCHAR,EOF);
AGAIN:OUTSTR("OUTPUT PIECE FILE NAME");
FILNAM← INCHWL;
ENTER(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("ENTER UNSUCCESSFUL !!!!");GO AGAIN ;END;
WORDOUT(DSKOUT,CVSIX("PIECE1"));
WORDOUT(DSKOUT,DUMMY);
WORDOUT(DSKOUT,DUMMY);

END"FILEOUT";

SIMPLE PROCEDURE FRAMEOUT;
BEGIN "FRAMEOUT"


WORDOUT(DSKOUT,BASPC);
WORDOUT(DSKOUT,SIZPC);
ARRYOUT(DSKOUT,AXX[PARB],SIZPC);
ARRYOUT(DSKOUT,AXY[PARB],SIZPC);
ARRYOUT(DSKOUT,RAD[PARB],SIZPC);
ARRYOUT(DSKOUT,ANGLE[PARB],SIZPC);
WORDOUT(DSKOUT,BD1);
WORDOUT(DSKOUT,BD2);
WORDOUT(DSKOUT,LOW1F);
WORDOUT(DSKOUT,HIGH1F);
WORDOUT(DSKOUT,LOW2F);
WORDOUT(DSKOUT,HIGH2F);
WORDOUT(DSKOUT,1);
WORDOUT(DSKOUT,BSIGN);
WORDOUT(DSKOUT,TERMB);
WORDOUT(DSKOUT,TERME);
WORDOUT(DSKOUT,CIRC);

END "FRAMEOUT";


⊃ START OF MK_PIECES ;

FILEOUT;
BASPC← 1 ; NPIECES← 0 ;
WHILE TRUE DO BEGIN
    OUTSTR("TYPE GROUP TO BE TESTED");
    IGRP←INTSCAN(STR←INCHWL,J);
    IF IGRP<1 OR IGRP>20 THEN DONE ;
    IF GRPSIZ[IGRP]≠0 THEN
      IF (SIZPC← AXIS(IGRP,BD1,BD2,LOW1F,LOW2F,HIGH1F,HIGH2F,BSIGN,PARB,PARE,
				TERMB,TERME,CIRC,AXX,AXY,AXZ,RAD,ANGLE))>0 THEN BEGIN
	  OUTSTR("PUT THIS PIECE ON THE OUTPUT FILE?"&CRLF);
	  IF INCHWL = "Y" THEN BEGIN
		IF BSIGN<0 THEN LOW2F↔HIGH2F; ⊃ To confirm to TOB Covention;
		⊃ NOW DO SOME CLEANING UP,NORMALISE INDICES TO LIE BET 1 AND SIZBD;
		SIZ1← SIZBD[BD1] ; SIZ2← SIZBD[BD2] ;
		IF LOW1F≤0 THEN LOW1F← LOW1F + SIZ1
				ELSE IF LOW1F>SIZ1  THEN LOW1F← LOW1F-SIZ1;
		IF HIGH1F<0 THEN HIGH1F← HIGH1F + SIZ1
				ELSE IF HIGH1F>SIZ1 THEN HIGH1F← HIGH1F-SIZ1;
		IF LOW2F≤0 THEN LOW2F← LOW2F + SIZ2
				ELSE IF LOW2F>SIZ2 THEN LOW2F← LOW2F-SIZ2;
		IF HIGH2F<0 THEN HIGH2F← HIGH2F + SIZ2
				ELSE IF HIGH2F>SIZ2 THEN HIGH2F← HIGH2F-SIZ2;
		FRAMEOUT ;
		BASPC← BASPC + SIZPC;
		NPIECES← NPIECES +1 ;
		END;
	  END;
END;

CLOSE(DSKOUT);


⊃ NOW REWRITE THE HEADER ;

L1:LOOKUP(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("LOOKUP OF "&FILNAM&" UNSUCCESSFUL !!!"); INCHWL; END;
L2:ENTER(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("ENTER OF "&FILNAM&" UNSUCCESSFUL !!!"); INCHWL; END;
USETI(DSKOUT,1);
ARRYIN(DSKOUT,FIXUP[1],'200);
FIXUP[2]← NPIECES;
FIXUP[3]← BASPC-1 ; ⊃ NOTE THIS IS  THE NO OF AXPTS UPTO THE POINT ;

USETO(DSKOUT,1);
ARRYOUT(DSKOUT,FIXUP[1],'200);

CLOSE(DSKOUT);

END "MK_PIECES" ;
⊃ MAIN PART OF OUTER ;

READPTS(FILNAM);
 READCAL;
BDIN;
READPTS;
GET_BREAKS;
GROUP;
AGAIN:  DISPLAY;
⊃  MOMENTS ;
MK_PIECES ; GO AGAIN;

END "OUTER" ;
COMMENT MAIN PROGRAM ;

SHORT INTEGER  LINE,I,J,NUM ;
STRING STR;
LABEL AGAIN ;
	
IF ¬RUNBEFORE THEN BEGIN
	POG ← 1 ;
	DISFINAL←DPYTST=0; COMMENT DISPLAY IF THIS IS A III;
	DISPOINTS←0;
	DISCONT←0;
	FMODE←1;
	LMODE←1;
	DPYMULT←3;
	WHICHFRAME←1;
	FIRSTONLY←FALSE;
	PDEN←0.1;
	SEGTHRESH← 12.0;
	RUNBEFORE←-1;END
ELSE DPYCLR;
STEPANGLE← PI/432 ;
FINALPOG←GETPOG; 
IF DISPOINTS THEN PPOG←GETPOG;
⊃ OPEN(OUTFILE,"DSK",0,0,3,0,0,EOF);
⊃ IF EOF ≠ 0 THEN OUTSTR("OPEN FAILS  ON OUTFILE"&CRLF);
⊃ ENTER(OUTFILE,"DEBUG1",FLAG);
⊃ IF FLAG THEN OUTSTR("ENTER FAILS ON OUTFILE"&CRLF);
OUTSTR("DEBUG?");
IF INCHWL = "Y" THEN DEBUG←TRUE ELSE DEBUG ← FALSE ;
OUTSTR("DISPLAY?");
IF INCHWL = "Y" THEN DISP←TRUE ELSE DISP ← FALSE;
FILIN(FILNAM);
OUTER;
END "CURVE" ;