SUBROUTINE GRAF1A ( XC, YC, NC, IC, ILINE, THICK, IPENC,
	1  XS, YS, NS, IS, DIA, ISYMB, IPENS )
C
C	 This routine is capable of producing a curve that will not pass
C	through any of the specified symbols.  The symbols need not be
C	located on the curve.
C	 There is no restriction on the order of the points in either
C	array, except of course that the curve be continuous.
C	For example, the curve may be mulivalued, and the symbols
C	may be randomly distributed.
C	 It assumes that the axes have already been drawn, for example
C	by GRAF1 or GRAF2.  The symbols are optionally drawn, but may
C       already be drawn, again by GRAF1 or GRAF2;  they are still
C       jumped over by the curve.
C	 Note that this routine is quite general with regard to being
C	able to draw lines which jump over regions.  To use an arbitrary
C	shape for the regions only requires that GRCHKX do that appropriate
C	arithmetic and return the entry and exit coordinates nearest X0.
C	This will even work for arbitrary polygons.  However this routine
C	only considers circular regions around the symbol markers.
C
C	XC, YC		Arrays which contain the coordinates for the
C			curve in units of the axes.
C	NC    		Number of points in the above arrays.
C	IC    		Interleave factor, ie every IC'th point is used.
C	ILINE		Type of line, if 0 then only symbols will be
C	     		generated, see GRDSHL for line types.
C	THICK		Thickness of the line in inches.
C	IPENC  		Pen colour to use when drawing the curve.
C	XS, YS		Arrays which contain the coordinates for the
C	      		symbols in units of the axes.
C	NS    		Number of points in the above arrays.
C			NS=0 means there are not symbols to plot.
C	IS    		Interleave factor, ie every IS'th point is used.
C	DIA   		Diameter of the symbols in inches, if 0. then
C	      		only the curve will be generated. If < 0 then
C	      		the symbols have already been drawn, but the
C	      		line will still avoid them.
C	ISYMB		Type of symbol, See SYMBOL.
C	IPENS  		Pen colour to use when drawing the symbols.
C
C
	DIMENSION XC(1), YC(1), XS(1), YS(1)
C
C GENERATE THE SYMBOLS IF DIA IS POSITIVE, OTHERWISE ASSUME THEY'RE ALREADY THERE
C
	IFF (DIA .GT. 0.0 .AND. NS .GT. 0)
	THEN	{
		CALL NEWPEN (IPENS)
		FOR K=1,NS,IS
			{
			CALL GRLOC (XP, YP, XS(K), YS(K))
			CALL SYMBOL (XP, YP, DIA, ISYMB, 0.0, -1)
			}
		}
C
C MOVE TO FIRST POINT
C
	DIAM = 1.4*DIA			!ALLOW SOME WHITE SPACE AROUND SYMBOLS
	CALL NEWPEN (IPENC)
	CALL GRLOC  (X1, Y1, XC(1), YC(1))
	CALL GRDSHL (X1, Y1, 3, ILINE, THICK)
	ICOUNT = IC+1			!INDEX INTO XC, YC
C
C LOOP THROUGH EACH POINT ON CURVE, GENERATE EACH LINE SEGMENT FROM (X0,Y0) TO (X1,Y1)
C
	WHILE (ICOUNT .LE. NC)
		{
		X0 = X1			!NEW LEFT POINT
		Y0 = Y1
		CALL GRLOC (X1, Y1, XC(ICOUNT), YC(ICOUNT))
C
C CHECK IF LINE SEGMENT INTERSECTS WITH ANY SYMBOL.
C IF SO, ENTRY AND EXIT LOCATIONS ARE RETURNED.
C
100		SX = SIGN( 1.0, X1-X0 )		!DIRECTION OF LINE SEGMENT
		CALL GRCHKX (X0, Y0, X1, Y1, SX, XIN, YIN, XOUT, YOUT,
	1			XS, YS, NS, IS, DIAM)
C
C DETERMINE HOW TO DRAW LINE SEGMENT
C
		IFF (ABS(XIN) .GT. 0.9E30)
		THEN	{		!NO INTERSECTION, JUST DRAW
			CALL GRDSHL (X1, Y1, 2, ILINE, THICK)
			}
		ELSE	{		!THERE WAS AN INTERSECTION
			IFF (XIN*SX .GT. X0*SX)
			THEN	{	!PART OF CURVE IS TO LEFT OF SYMBOL
				CALL GRDSHL (XIN, YIN, 2, ILINE, THICK)
				}
			IFF (XOUT*SX .LT. X1*SX)
			THEN	{	!MORE OF CURVE AFTER SYMBOL
				X0 = XOUT + 1.0E-7*SX
				Y0 = YOUT
				CALL GRDSHL (X0, Y0, 3, ILINE, THICK)
				GO TO 100	!SEE IF REMAINING SEGMENT OK
				}
			}
		ICOUNT = ICOUNT + IC
		}
	CALL PLOTF
	RETURN
	END
	SUBROUTINE GRCHKX (X0, Y0, X1, Y1, SX, XIN, YIN, XOUT, YOUT,
	1			XS, YS, NS, IS, DIAM)
C
C This is an internal routine used by GRAF1A which determines if the
C line segment from (X0,Y0) to (X1,Y1) intersects with any of the
C plotting symbols located at (Xs,Ys) of diameter DIAM.
C If there is no intersection, XIN is returned as a big number,
C otherwise the points of entry and exit are returned in
C (XIN,YIN) and (XOUT,YOUT).
C  This assumes that the symbols are circular in shape, but this
C routine can use any general shape as long as it returns the
C first entry and exit points along the line.
C  Note that X1 may be changed slightly by this routine if it is
C very close to X0.  This works correctly for any arbitrary line
C orientation.
C
	REAL XS(1), YS(1)
C
C AVOID CASE OF VERTICAL LINE WITH INFINITE SLOPE
C
	IF (ABS(X1-X0) .LT. 1.0E-5) X1 = X0 + 1.0E-5*SX

	XIN = 1.0E30*SX
	SLOPE = (Y1-Y0) / (X1-X0)
	R2 = (DIAM/2.0)**2
	I = 1
C
C DETERMINE QUADRADRIC EQUATION OF INTERSECTION OF LINE AND CIRCLE
C
	WHILE (I .LE. NS)
		{
		CALL GRLOC (XC, YC, XS(I), YS(I))
		X = X0-XC
		Y = Y0-YC
		A = SLOPE*SLOPE + 1.0
		B = SLOPE*Y + X
		C = X*X + Y*Y - R2
		D = B*B - A*C
C
C IF DETERMINANT IS POSITIVE, SOLUTION EXISTS
C
		IFF (D .GT. 1.0E-4)
		THEN	{
			XL0 = (-B-SQRT(D)*SX) / A
			XR0 = (-B+SQRT(D)*SX) / A
			XL = XL0 + X0
			XR = XR0 + X0
			IFF (XL*SX .LT. XIN*SX
	1		 .AND. (XL*SX .LT. X1*SX .AND. XR*SX .GT. X0*SX) )
			THEN	{
				XIN  = XL
				YIN  = SLOPE*XL0 + Y0
				XOUT = XR
				YOUT = SLOPE*XR0 + Y0
				}
			}
		I = I+IS
		}
	RETURN
	END
	SUBROUTINE GRDSHL (X, Y, IPEN, ITYPE, THICK)
C
C	This routine will draw various styles of lines
C	depending upon the value of ITYPE.  If ITYPE equals:
C
C	0 - No Line
C	1 - Solid
C	2 - Short dashes
C	3 - Long and Short dashes
C	4 - Long dashes
C	5 - Long and 2 short dashes
C	6 - Dotted
C	7 - Long and 3 short dashes
C	8 - 2 Long and 1 short dashes
C	9 - 2 Long and 2 Short dashes
C	10 - 2 Long and 3 Short dashes
C
C	IPEN equal to +/-3 will reset the arc length gauge to zero.
C	THICK is the thickness of the line in inches.
C
C This routine is a modified version of [370,30]DSHPLT obtained 
C from the FALL 1980 RSX DECUS Symposium tape, in a submission
C by Bob Grandle from NASA Langley .
C
C
	INTEGER UPDOWN(29), START(10), LENGTH(10)
	DATA UPDOWN/ 8*2, 3, 8*2, 3, 2,2, 3, 2,2, 3, 2,2, 3, 2, 3/
	DATA START/ 0, 25, 10, 10, 10, 28, 10, 1, 1, 1/
	DATA LENGTH/ 0, 3, 12, 9, 15, 2, 18, 21, 24, 27/
	DATA ARC, XLAST, YLAST/3*0./, I/0/, JLAST/0/
	DATA DARC/ 0.06 /	!LENGTH OF EACH LINE SEGMENT IN INCHES
C
C - IF PEN UP, RESET LINE PARAMETERS AND MOVE.
	IF( IABS(IPEN) .NE. 3 ) GO TO 100
	ARC = 0.0
	I = 0
	CALL GRTHKL (X, Y, IPEN, THICK)
	GO TO 1200
C
C - NO LINE?
100	IF( ITYPE .NE. 0 ) GO TO 200
	CALL GRTHKL (X, Y, 3, THICK)
	GO TO 1200
C
C - WHICH TYPE OF LINE?
200	JTYPE = MOD( ITYPE-1, 10 ) + 1
	IF( JTYPE .NE. 1 ) GO TO 300
	CALL GRTHKL (X, Y, IPEN, THICK)
	GO TO 1200
C
C - DASHED LINE.
300	DELTAX = X - XLAST
	DELTAY = Y - YLAST
	DELTA = SQRT( DELTAX*DELTAX + DELTAY*DELTAY + 1.E-30 )
	JSTYLE = START( JTYPE )
	NUMBER = LENGTH( JTYPE )
	IF( JTYPE .NE. JLAST ) I = 0
	JLAST = JTYPE
	IF( DELTA .GE. ARC ) GO TO 400
	ARC = ARC-DELTA
	CALL GRTHKL (X, Y, UPDOWN(JSTYLE+I), THICK)
	GO TO 1200
C
400	COSTH = DELTAX/DELTA
	SINTH = DELTAY/DELTA
	DX = ARC*COSTH+XLAST
	DY = ARC*SINTH+YLAST
	CALL GRTHKL (DX, DY, UPDOWN(JSTYLE+I), THICK)
	XLAST = DX
	YLAST = DY
	DELTA = DELTA-ARC
	ARC = 0.0
	I = MOD( I+1, NUMBER )
C
	INT = DELTA/DARC
	IF( INT .LT. 1 ) GO TO 1100
C
	DX = DARC*COSTH
	DY = DARC*SINTH
	DO 1000 J=1,INT
	XLAST = XLAST+DX
	YLAST = YLAST+DY
	CALL GRTHKL (XLAST, YLAST, UPDOWN(JSTYLE+I), THICK)
	DELTA = DELTA-DARC
1000	I = MOD( I+1, NUMBER )
C
1100	ARC = DARC - DELTA
	CALL GRTHKL (X, Y, UPDOWN(JSTYLE+I), THICK)
C
1200	XLAST = X
	YLAST = Y
	IF( IPEN .GT. 0 ) RETURN
	CALL GRTHKL (X, Y, -3, THICK)
	XLAST = 0.0
	YLAST = 0.0
	RETURN
	END
	SUBROUTINE GRTHKL (X1, Y1, IPEN, THICK)
C
C Draws a line from the present pen position to point (X1, Y1) in inches
c of thickness THICK inches.  the parameter IPEN is like its counterpart
c in the PLOT subroutine, IPEN=3 to lift pen before the move, IPEN=2 to
c lower the pen before the move.  Unlike PLOT, this will use multiple
c pen strokes to make a thicker line.  The width of the pen is entered
c here as a constant and should be the thinnest pen likely to be used.
C
	DATA PENWID / 0.015 /	!PEN WIDTH IN INCHES
C
	IF (THICK .LE. 0.0) GO TO 1000
	IF (IABS(IPEN) .NE. 2) GO TO 1000
C
C CALCULATE WIDTH OF PEN STROKES AND NUMBER OF STROKES.
C DX AND DY FORM A VECTOR ORTHOGONAL TO THE DIRECTION OF THE
C LINE WHICH IS THE DIRECTION OF THE SHIFT FOR EACH STROKE.
C
	R = SQRT( (X1-X0)**2 + (Y1-Y0)**2 )
	IF (R .EQ. 0.0) GO TO 1000			!ZERO LENGTH LINE
	NSTROK = ( THICK/PENWID - 1.0 ) * 0.5 + 0.9	!# OF EXTRA STROKE PAIRS
	PENWI2 = THICK / ( 2.0*FLOAT(NSTROK) + 1.0 )	!EFFECTIVE PEN WIDTH
	DX = -PENWI2 * (Y1-Y0) / R
	DY =  PENWI2 * (X1-X0) / R
C
C DRAW CENTRAL LINE
C
	CALL PLOT (X0, Y0, 3)
	CALL PLOT (X1, Y1, 2)
C
C DRAW THICK LINE USING PAIRS OF LINES, ONE ON EACH SIDE OF CENTER
C
	DXS = 0.0
	DYS = 0.0
	IF (NSTROK .LE. 0) GO TO 1000	!ONE STROKE SUFFICIENT
	DO 100 J=1,NSTROK
		DXS = DXS+DX
		DYS = DYS+DY
		CALL PLOT (X1+DXS, Y1+DYS, 2)
		CALL PLOT (X0+DXS, Y0+DYS, 2)
		CALL PLOT (X0-DXS, Y0-DYS, 2)
		CALL PLOT (X1-DXS, Y1-DYS, 2)
100	CONTINUE
C
C MOVE TO END OF VECTOR, REMEMBER NEW STARTING POINT
C
1000	CALL PLOT (X1, Y1, IPEN)
	X0 = X1
	Y0 = Y1
	RETURN
	END
	SUBROUTINE GRSPLI (XOUT, YOUT, NOUT, XOUT0, XOUT1, XIN, YIN, NIN, WORK)
C
C This takes a set of points XIN,YIN, fits a piecewise cubic spline to
C the points, and puts the curve into XOUT,YOUT.  A cubic spline is a
C cubic equation which passes through the 4 nearest points.  Thus the
C spline curve will pass through every point and will be fairly smooth.
C This is not the same as a least squares fit, which does not have to
C pass through the points.  This is best suited to fairly smooth data
C points which are fairly sparse; if you get too many points (say >20),
C then the curve will look too jagged.
C  This routine can also extrapolate outside the range of XIN,YIN, but
C since it fits a cubic, it quickly gets bent, so don't extrapolate too
C far.  In fact, if you want a nice smooth curve at both ends, add an
C extra point beyond each end, which will define the slope of the curves
C as they pass through the inner points.
C
C XOUT, YOUT  Arrays which will contain the points defining the cubic
C             spline which is fitted to the data.  These must be
C             dimensioned to at least NOUT elements.
C NOUT        Number of points to generate in the fitted curve.
C XOUT0,XOUT1 Minimum and maximum x values for which to fit the curve.
C             Thus XOUT(1) = XOUT0 and XOUT(NOUT) = XOUT1.
C XIN, XIN    The data points to which the curve will be fitted.
C             These must be provided.  These must be dimensioned to
C             at least NIN elements.
C NIN         Number of data points provided, must be at least 4.
C WORK        A work array which must be dimensioned to at least 3*NIN.
C
	REAL XOUT(1), YOUT(1), XIN(1), YIN(1), WORK(1)
C
	DX = (XOUT1-XOUT0) / (FLOAT(NOUT)-1.0)
	X = XOUT0
	INTRVL = -1
C
C FIT POINTS TO CURVE
C
	DO 100 J=1,NOUT
		CALL GRPCS (XIN, YIN, NIN, X, Y, SLOPE, INTRVL, WORK)
		XOUT(J) = X
		YOUT(J) = Y
		X = X+DX
100	CONTINUE
	RETURN
	END
	SUBROUTINE GRPCS ( X, Y, N, X0, FX0, FPX0, INTRVL, WORK )
C
C	This routine will interpolate a set of points using a
C	piecewise-cubic spline.
C	The minimum N should be 4.
C
C	X	Independent array.
C	Y	Dependent array.
C	N	Number of points in the above arrays.
C	X0	Value at which to calculate the interpolate.
C	FX0	Interpolated value.
C	FPX0	Derivative of the interpolate at X0.
C	INTRVL	Equals interval number of where X0 is.
C		Set to 0 or N if extrapolation, must equal -1
C		on first call to initialize interpolation.
C	WORK	Working array, dimensioned at least 3*N.
C
C This routine is a modified version of [370,30]PCS obtained 
C from the FALL 1980 RSX DECUS Symposium tape, in a submission
C by Bob Grandle from NASA Langley .
C
	DIMENSION X(1), Y(1), WORK(1)
C
C - INTERPOLATION ALREADY CALCULATED.
	IF( INTRVL .GE. 0 ) GO TO 5100
C
C - CONSTABLES.
	N2 = 0
	N3 = N
	N4 = N + N3
	ND = N3
	NS = N4
C
C - CALCULATE SLOPES FOR END POINTS.
	CALL GRLSQF ( X, Y, 4, WORK(1), 3, WORK(5), WORK(11) )
	S1 = WORK(2) + X(1)*( 2.*WORK(3) + 
	1  3.*X(1)*WORK(4) )
C
	CALL GRLSQF ( X( N-3 ), Y( N-3 ), 4, WORK(1), 3, WORK(5),
	1  WORK(11) )
	WORK( N2+N ) = WORK(2) + X(N)*( 2.*WORK(3) +
	1  3.*X(N)*WORK(4) )
	WORK( N2+1 ) = S1
	WORK( ND+1 ) = 0.
	WORK( NS+1 ) = 1.
C
C - CALCULATE DIFFERENCES AND SLOPES.
	DO 1000 I=2,N
	WORK( ND+I ) = X( I ) - X( I-1 )
1000	WORK( NS+I ) = ( Y( I ) - Y( I-1 ) )/WORK( ND+I )
C
	DO 2000 I=2,N-1
	WORK( N2+I ) = 3.*( WORK( ND+I)*WORK( NS+I+1 ) +
	1  WORK( ND+I+1 )*WORK( NS+I ) )
2000	WORK( NS+I ) = 2.*( WORK( ND+I ) + WORK( ND+I+1 ) )
C
	DO 3000 I=2,N-1
	G = -WORK( ND+I+1)/WORK( NS+I-1 )
	WORK( NS+I ) = WORK( NS+I ) + G*WORK( ND+I-1 )
3000	WORK( N2+I ) = WORK( N2+I ) + G*WORK( N2+I-1 )
C
	DO 4000 I=N-1,1,-1
4000	WORK( N2+I ) = ( WORK( N2+I ) - WORK( ND+I )*WORK( N2+I+1 ) )/
	1  WORK( NS+I )
C
	DO 5000 I=1,N-1
	DX = X( I+1 ) - X( I )
	F1 = ( Y( I+1 ) - Y( I ) )/DX
	F3 = WORK( N2+I ) + WORK( N2+I+1 ) - 2.*F1
	WORK( N3+I ) = ( F1 - WORK( N2+I ) - F3 )/DX
5000	WORK( N4+I ) = F3/( DX*DX )
	INTRVL = 1
C
C - FIND WHICH INTERVAL POINT IS WITHIN.
5100	IF( X0 .GT. X(1) ) GO TO 5200
	INTRVL = 0
	ISEG = 1
	GO TO 5500
5200	IF( X0 .LT. X(N) ) GO TO 5300
	INTRVL = N
	ISEG = N-1
	GO TO 5500
C
C - MOVE TO PROPER SEGMENT.
5300	INTRVL = MIN0( N-1, MAX0( 1, INTRVL ) )
	IF( ( X(INTRVL) .LE. X0 ) .AND. ( X0 .LE. X(INTRVL+1) ) )
	1  GO TO 5400
	INTRVL = INTRVL + SIGN( 1., X0 - X(INTRVL) )
	GO TO 5300
C
C - CALCULATE INTEROPLATE AND DERIVATIVE.
5400	ISEG = INTRVL
5500	DX = X0 - X( ISEG )
	FX0 = Y( ISEG ) + DX*( WORK( N2+ISEG ) +
	1  DX*( WORK( N3+ISEG ) + DX*WORK( N4+ISEG ) ) )
	FPX0 = WORK( N2+ISEG ) + DX*( 2.*WORK( N3+ISEG ) +
	1  DX*3.*WORK( N4+ISEG ) )
C
	RETURN
	END
	SUBROUTINE GRLSQF (X,Y,NPTS,COEFF,N,ERROR,WORK)
C
C	This subroutine calculates the coefficients of the Nth order
C	polynomial which provides the best least-squares fit to the
C	set of points, (X,Y).
C
C	X	Array of independent variables.
C	Y	Array of dependent variables.
C	NPTS	In absolute value, the number of data points.  If negative
C		the error will be expressed as a percentage.
C	COEFF	Coefficients of the fitted polynomial, stored such that
C		Y(X) = COEFF(1) + ... + COEFF(N+1)*X**N.  COEFF must
C		be dimensioned at least N+1.
C	N	Order of the polynomial, N < 6.
C	ERROR	Array containing the differences between the actual
C		values and those calculated by the polynomial.  
C		Percentage error are stored if NPTS < 0.  Also
C		stored are the maximum error in ERROR( IABS(NPTS) + 1 )
C		and the sigma error in ERROR( IABS(NPTS) + 2 ).  ERROR
C		must be dimensioned at least IABS(NPTS)+2.
C	WORK	Working array, dimensioned at least 2*IABS(NPTS).
C
C This routine is a modified version of [370,30]ORTHO obtained 
C from the FALL 1980 RSX DECUS Symposium tape, in a submission
C by Bob Grandle from NASA Langley .
C
	REAL A(6,6),B(6),C(6),D(6),S(6)
	REAL X(1),Y(1),COEFF(1),ERROR(1),WORK(1)
	MPTS = IABS(NPTS)
	NORDER = N+1
C
C - ZERO WORKING ARRAYS.
	DO 1000 J=1,NORDER
	B(J) = 0.0
	COEFF(J) = 0.0
	D(J) = 0.0
	S(J) = 0.0
	DO 1000 I=1,NORDER
1000	A(I,J) = 0.0
C
C - CALCULATE SUM OF X AND Y VALUES FOR AVERAGE.
	DO 2000 I=1,MPTS
	B(1) = B(1)+X(I)
2000	D(1) = D(1)+Y(I)
	C(1) = 0.0
	S(1) = MPTS
	D(1) = D(1)/S(1)
C
C - CALCULATE DIFFERENCE (ERROR) BETWEEN Y VALUES AND Y AVERAGE.
	DO 3000 I=1,MPTS
3000	ERROR(I) = Y(I)-D(1)
	IF( NORDER .EQ. 1 ) GO TO 7100
	B(1) = B(1)/S(1)
C
C - CALCALUTE FIRST ORDER MOMENT ARMS.
	DO 4000 I=1,MPTS
	WORK( MPTS+I ) = 1.
4000	WORK(I) = X(I)-B(1)
	J=1
4100	J = J+1
C
C - CALCULATE JTH ORDER MOMENTS.
	DO 5000 I=1,MPTS
	D(J) = D(J)+ERROR(I)*WORK(I)
	B(J) = B(J)+X(I)*WORK(I)**2
5000	S(J) = S(J)+WORK(I)**2
	D(J) = D(J)/S(J)
C
C - CALCULATE ERROR FROM JTH ORDER CURVE.
	DO 6000 I=1,MPTS
6000	ERROR(I) = ERROR(I)-D(J)*WORK(I)
	IF( J .EQ. NORDER ) GO TO 7100
	B(J) = B(J)/S(J)
	C(J) = S(J)/S(J-1)
C
C - CALCULATE NEXT ORDER MOMENT ARMS.
	DO 7000 I=1,MPTS
	WW = WORK(I)
	WORK(I) = (X(I)-B(J))*WORK(I)-C(J)*WORK( MPTS+I )
7000	WORK( MPTS+I ) = WW
	GO TO 4100
C
C - SET UP COEFFICIENT ARRAY.
7100	A(1,1) = 1.0
	IF( NORDER .EQ. 1 ) GO TO 8100
	A(2,2) = A(1,1)
	A(2,1) = -B(1)*A(1,1)
	IF( NORDER .EQ. 2 ) GO TO 8100
	DO 8000 J=3,NORDER
	A(J,1) = -B(J-1)*A(J-1,1)-C(J-1)*A(J-2,1)
	DO 8000 M=2,NORDER
8000	A(J,M) = A(J-1,M-1)-B(J-1)*A(J-1,M)-C(J-1)*A(J-2,M)
C
C - REDUCE COEFFICIENT MATRIX TO OBTAIN COEFFICIENTS.
8100	DO 9000 I=1,NORDER
	DO 9000 J=I,NORDER
9000	COEFF(I) = COEFF(I)+D(J)*A(J,I)
C
C - CALCULATE ERROR.
	ERROR(MPTS+1) = -1.
	ERROR(MPTS+2) = 0.0
	DO 10000 I=1,MPTS
	ERROR(MPTS+2) = ERROR(MPTS+2)+ERROR(I)**2
	IF( NPTS .GT. 0 ) GO TO 9200
C
C - CALCULATE PERCENTAGE ERROR.
	IF( ABS(Y(I)) .GT. 1.E-6 ) GO TO 9100
	ERROR(I) = 1.E10
	GO TO 10000
9100	ERROR(I) = 1.0-ERROR(I)/Y(I)
9200	IF( ABS(ERROR(I)) .LE. ERROR(MPTS+1) ) GO TO 10000
C
C - FIND MAXIMUM ERROR.
	ERROR(MPTS+1) = ABS(ERROR(I))
10000	CONTINUE
C
C - CACULATE SIGMA.
	ERROR(MPTS+2) = SQRT(ERROR(MPTS+2)/MPTS)
C
	RETURN
	END