SUBROUTINE LOWESS (X,Y,N, F, NSTEPS, DELTA, YS, RW, RES)

C+-------------------------------------------------------------------------
C
C LOWESS.FOR          LOWESS Curve fitting
C
C Program family:     Numerical and Statistical
C For use by:         Public
C Keywords:           get,fit,lowess,curve,scatterplot,regression,distribution
C Author:             William Cleveland
C
C Complete Description:
C    LOWESS -- Locally weighted scatterplot smoother.
C
C Logicals used:      None.
C Files used:         None.
C
C Special Notes:      
C    Cleveland, W.S.  1979.  Robust locally weighted regression and smoothing
C         scatterplots.  J.A.S.A. 74:829-836
C    Cleveland, W.S.  1985.  Graphical perception and graphical methods for
C         analyzing scientific data.  Science.  229:828-833
C
C Modification History:
C  Date       Person    Modification
C  ---------  --------  ---------------------------------------------------
C  22 Dec 93  Larkin    Start describing the arguments in comments
C			Remove type declarations for IFIX,FLOAT, MIN0, MAX0.
C--------------------------------------------------------------------------
C+
C       Arguments not described in original source; description below by R.L.
C
C Input.  Length of all arrays (all of which are arguments)
	INTEGER N	    
C Input.  Independent variable (abcissa)
	REAL	X(N)	   
C Input.  Dependent variable (ordinate)
	REAL	Y(N)	   
C Input.  "the neighborhood of influential points"
C	  Increasing F --> more smoothing
C	  Used (only) in setting NS.
C         F*N is upper bound of NS.
C         Cleveland recommends F be in the range
C         0.2 to 0.8 & to use F=0.5 as a starting
C	  value with new data.
	REAL	F	   
C Input.  Number of iterations ("t" in Cleveland 
C         1979).  He recommends NSTEPS=2 (p. 834).
	INTEGER NSTEPS	   
C Input.  In units of X.  I'm going to assume
C         DELTA is the "d" of Cleveland 1979 p.833
C	  and that therefore DELTA=1 is good.
	REAL	DELTA	   
C Output. Smoothed Y (I'm pretty sure)
	REAL	YS(N)	   
C Output. Robust weights?
	REAL	RW(N)	   
C Output. Residuals,  Y(i) - Yhat(i)
	REAL	RES(N)	   
C-
	INTEGER NRIGHT,I,J
	INTEGER ITER,LAST,M1,M2,NS,NLEFT
	REAL ABS,CUT,CMAD,R,D1,D2
	REAL C1,C9,ALPHA,DENOM
	LOGICAL OK

	IF (N .GE. 2) GOTO 1
	  YS(1) = Y(1)
	  RETURN
C We have AT LEAST TWO, AT MOST N POINTS
   1	NS = MAX0(MIN0(IFIX(F*FLOAT(N)),N),2)
	ITER = 1
 	  GOTO 3
   2	  ITER = ITER+1
   3	  IF (ITER .GT. NSTEPS+1) GOTO 22   
C ROBUSTNESS ITERATIONS
	  NLEFT = 1
	  NRIGHT = NS
C INDEX OF PREV ESTIMATED POINT
	  LAST = 0
C INDEX OF CURRENT POINT
	  I = 1
   4	    IF (NRIGHT .GE. N) GOTO 5
C MOVE NLEFT,NRIGHT TO RIGHT IF RADIUS DECREASES
    	      D1 = X(I)-X(NLEFT)
C IF D1<=D2 WITH X(NRIGHT+1) == X(NRIGHT),LOWEST FIXES
	      D2 = X(NRIGHT+1)-X(I)
	      IF (D1 .LE. D2) GOTO 5
C RADIUS WILL NOT DECREASE BY MOVE RIGHT
	      NLEFT = NLEFT+1
	      NRIGHT = NRIGHT+1
	      GOTO 4
C FITTED VALUE AT X(I)
   5	    CALL LOWEST (X,Y, N, X(I), YS(I), NLEFT, NRIGHT, RES,
     1         (ITER.GT.1), RW, OK)
	    IF (.NOT. OK) YS(I) = Y(I)
C ALL WEIGHTS ZERO - COPY OVER VALUE (ALL RW==0)
	    IF (LAST .GE. I-1) GOTO 9
	      DENOM = X(I)-X(LAST)
C SKIPPED POINTS -- INTERPOLATE
C NON-ZERO - PROOF?
      	      J = LAST+1
	        GOTO 7
   6		J = J+1
   7		IF (J .GE. I) GOTO 8
		ALPHA = (X(J)-X(LAST))/DENOM
		YS(J) = ALPHA*YS(I)+(1.0-ALPHA)*YS(LAST)
		GOTO 6
   8	      CONTINUE
C LAST POINT ACTUALLY ESTIMATED
   9        LAST = I
C X COORD OF CLOSE POINTS
	    CUT = X(LAST)+DELTA
	    I = LAST+1
	      GOTO 11
  10	      I = I+1
  11	      IF (I .GT. N) GOTO 13
C FIND CLOSE POINTS
	      IF (X(I) .GT. CUT) GOTO 13
C I ONE BEYOND LAST PT WITHIN CUT
	      IF (X(I) .NE. X(LAST)) GOTO 12
		YS(I) = YS(LAST)
C EXACT MATCH IN X
		LAST = I
 12	      CONTINUE
	      GOTO 10
C BACK 1 POINT SO INTERPOLATION WITHIN DELTA, BUT ALWAYS GO FORWARD
 13	    I = MAX0(LAST+1,I-1)
 14	    IF (LAST .LT. N) GOTO 4
C RESIDUALS
	  DO 15 I = 1,N
	    RES(I) = Y(I) - YS(I)
 15	    CONTINUE
	  IF (ITER .GT. NSTEPS) GOTO 22
C COMPUTE ROBUSTNESS WEIGHTS EXCEPT LAST TIME
	  DO 16 I = 1,N
	    RW(I) = ABS(RES(I))
 16	    CONTINUE
	  CALL SSORT(RW,N)
	  M1 = N/2+1
	  M2 = N-M1+1
C 6 MEDIAN ABS RESID
	  CMAD = 3.0*(RW(M1)+RW(M2))
	  C9 = .999*CMAD
	  C1 = .001*CMAD
	  DO 21 I = 1,N
	    R = ABS(RES(I))
	    IF (R .GT. C1) GOTO 17
	      RW(I) = 1.
C NEAR 0, AVOID UNDERFLOW
	      GOTO 20
 17	      IF (R .LE. C9) GOTO 18
		RW(I) = 0.
C NEAR 1, AVOID UNDERFLOW
		GOTO 19
 18		RW(I) = (1.0-(R/CMAD)**2)**2
 19	    CONTINUE
 20	    CONTINUE
 21	    CONTINUE
	  GOTO 2

 22	RETURN
	END

C LOWEST -- LOWESS ESTIMATE AT SINGLE X VALUE
	SUBROUTINE LOWEST (X,Y, N, XS, YS, NLEFT, NRIGHT, W, USERW, RW, OK)

	INTEGER N
	INTEGER NLEFT,NRIGHT
	REAL X(N),Y(N),XS,YS,W(N),RW(N)
	LOGICAL USERW,OK
	INTEGER NRT,J
	REAL ABS,A,B,C,H,R
	REAL H1,SQRT,H9,AMAX1,RANGE

	RANGE = X(N)-X(1)
	H = AMAX1(XS-X(NLEFT),X(NRIGHT)-XS)
	H9 = .999*H
	H1 = .001*H
C SUM OF WEIGHTS
	A = 0.0
	J = NLEFT
	  GOTO 2
  1	  J = J+1
  2	  IF (J .GT. N) GOTO 7
C COMPUTE WEIGHTS (PICK UP ALL TIES ON RIGHT)
	  W(J) = 0.
	  R = ABS(X(J)-XS)
	  IF (R .GT. H9) GOTO 5
	    IF (R .LE. H1) GOTO 3
	      W(J) = (1.0-(R/H)**3)**3
C SMALL ENOUGH FOR NON-ZERO WEIGHT
	      GOTO 4
  3	      W(J) = 1.
  4	    IF (USERW) W(J) = RW(J)*W(J)
	    A = A+W(J)
	    GOTO 6
  5	    IF (X(J) .GT. XS) GOTO 7
C GET OUT AT FIRST ZERO WT ON RIGHT
  6	  CONTINUE
	  GOTO 1
C RIGHTMOST PT (MAY BE GREATER THAN NRIGHT BECAUSE OF TIES)
  7	NRT = J-1
	IF (A .GT. 0.0) GOTO 8
	  OK = .FALSE.
	  GOTO 16
  8	  OK = .TRUE.
C WEIGHTED LEAST SQUARES
	  DO 9 J = NLEFT,NRT
C MAKE SUM OF W(J) ==1
	    W(J) = W(J)/A
  9	    CONTINUE
	  IF (H .LE. 0.) GOTO 14
	    A = 0.0
C USE LINEAR FIT
	    DO 10 J = NLEFT,NRT
C WEIGHTED CENTER OF X VALUES
	      A = A+W(J)*X(J)
 10	      CONTINUE
	    B = XS-A
	    C = 0.0
	  DO 11 J = NLEFT,NRT
	    C = C+W(J)*(X(J)-A)**2
 11	    CONTINUE
	  IF (SQRT(C) .LE. .001*RANGE) GOTO 13
	    B = B/C
C POINTS ARE SPREAD OUT ENOUGH TO COMPUTE SLOPE
	    DO 12 J = NLEFT,NRT
	      W(J) = W(J)*(B*(X(J)-A)+1.0)
 12	      CONTINUE
 13	    CONTINUE
 14	  YS = 0.0
	  DO 15 J = NLEFT,NRT
	    YS = YS+W(J)*Y(J)
 15	    CONTINUE
 16	RETURN
	END

C SORTP

	SUBROUTINE SORTP (A, N, B)

	REAL	A(N)
	INTEGER B(N)
	INTEGER	IU(16),IL(16)
	INTEGER P

  1	I = 1
	J = N
	M = 1
  5	IF(I .GE. J) GO TO 70
C FIRST ORDER A(I),A(J),A((I+J)/2), AND USE MEDIAN TO SPLIT THE DATA
 10	K = I
	IJ = (I+J)/2
	T=A(IJ)
	IT=B(IJ)
	IF(A(I) .LE. T) GOTO 20
	A(IJ)=A(I)
	B(IJ)=B(I)
	A(I)=T
	B(I)=IT
	T=A(IJ)
	IT=B(IJ)
 20	L=J
	IF(A(J) .GE. T) GOTO 40
	A(IJ) = A(J)
	B(IJ) = B(J)
	A(J)=T
	B(J)=IT
	T=A(IJ)
	IT=B(IJ)
	IF(A(I) .LE.T) GOTO 40
	A(IJ) =A(I)
	B(IJ)=B(I)
	A(I)=T
	B(I)=IT
	T=A(IJ)
	IT=B(IJ)
	GOTO 40
 30	A(L) = A(K)
	B(L) = B(K)
	A(K)=TT
	B(K)=ITT
 40	L=L-1
	IF(A(L) .GT. T) GOTO 40
	TT=A(L)
	ITT=B(L)
C SPLIT THE DATA INTO A(I TO L) .LT. T, A(K TO J) .GT. T
 50	K =K+1
	IF (A(K) .LT. T) GOTO 50
	IF (K .LE. L) GOTO 30
	P = M
	M = M+1
C SPLIT THE LARGER OF THE SEGMENTS
	IF (L-I .LE. J-K) GOTO 60
	IL(P)=I
	IU(P)=L
	I=K
	GOTO 80
 60	IL(P)=K
	IU(P)=J
	J=L
	GOTO 80
 70	M=M-1
	IF(M .EQ. 0) RETURN
	I = IL(M)
	J = IU(M)
C SHORT SECTIONS ARE SORTED BY BUBBLE SORT
 80	IF (J-I .GT. 10) GOTO 10
	IF (I .EQ. 1) GOTO 5
	I =I-1
 90	I =I+1
	IF(I .EQ. J) GOTO 70
	T = A(I+1)
	IT = B(I+1)
	IF(A(I) .LE. T) GOTO 90
	K=I
 100	A(K+1)=A(K)
	B(K+1)=B(K)
	K=K-1
	IF (T .LT. A(K)) GOTO 100
	A(K+1)=I
	B(K+1)=IT
	GOTO 90

	END
C
C SSORT

	SUBROUTINE SSORT(A,N)

C SORTING BY HOARE METHOD, C.A.C.M. (1961) 321, MODIFIED BY SINGLETON,
C C.A.C.M. (1969) 185.
	REAL	A(N)
	INTEGER	IU(16),IL(16)
	INTEGER P

	I =1
	J = N
	M = 1
  5	IF (I.GE.J) GOTO 70
C FIRST ORDER A(I),A(J),A((I+J)/2), AND USE MEDIAN TO SPLIT THE DATA
 10	K=I
	IJ=(I+J)/2
	T=A(IJ)
	IF(A(I) .LE. T) GOTO 20
	A(IJ)=A(I)
	A(I)=T
	T=A(IJ)
 20	L=J
	IF(A(J).GE.T) GOTO 40
	A(IJ)=A(J)
	A(J)=T
	T=A(IJ)
	IF(A(I).LE.T) GOTO 40
	A(IJ)=A(I)
	A(I)=T
	T=A(IJ)
	GOTO 40
 30	A(L)=A(K)
	A(K)=TT
 40	L=L-1
	IF(A(L) .GT. T) GOTO 40
	TT=A(L)
C SPLIT THE DATA INTO A(I TO L) .LT. T, A(K TO J) .GT. T
 50	K=K+1
	IF(A(K) .LT. T) GOTO 50
	IF(K .LE. L) GOTO 30
	P=M
	M=M+1
C SPLIT THE LARGER OF THE SEGMENTS
	IF (L-I .LE. J-K) GOTO 60
	IL(P)=I
	IU(P)=L
	I=K
	GOTO 80
 60 	IL(P)=K
	IU(P)=J
	J=L
	GOTO 80
 70	M=M-1
	IF(M .EQ. 0) RETURN
	I =IL(M)
	J=IU(M)
C SHORT SECTIONS ARE SORTED BY BUBBLE SORT
 80	IF (J-I .GT. 10) GOTO 10
	IF (I .EQ. 1) GOTO 5
	I=I-1
 90	I=I+1
	IF(I .EQ. J) GOTO 70
	T=A(I+1)
	IF(A(I) .LE. T) GOTO 90
	K=I
 100	A(K+1)=A(K)
	K=K-1
	IF(T .LT. A(K)) GOTO 100
	A(K+1)=T
	GOTO 90

	END