pro splac,x,y,ytol,u,v,w,noband=noband,discon=discon,verbose=verbose, _extra=e
;+
;procedure	splac
;	compute Segmented Piecewise Linear Approximation to Curve.
;	generates piecewise linear approximations with fewest line
;	segments within given tolerances.
;
;about the name
;	have you heard the one about the Stanford Linear Accelerator (SLAC)
;	and what would happen after the Big One hits?  Obviously, it would
;	have to be renamed SPLAC, the Stanford Piecewise Linear Accelerator..
;
;syntax
;	splac,x,y,ytol,u,v,w,/noband,/discon,verbose=verbose
;
;parameters
;	x	[INPUT; required] points at which data are defined
;		* must be in ascending order
;	y	[INPUT; required] values of data array, Y(X)
;		* must match size of X
;	ytol	[INPUT; required] tolerance of approximation: allowed
;		>deviation< from the true curve
;		* if scalar, single absolute-value tolerance for entire curve.
;		* if vector, tolerance at each X -- size must match X
;		* if -ve, ABS(YTOL) is taken to be the fractional factor
;		  relative to the true value.
;	u	[OUTPUT] a partition of the interval [X(0),X(N[X]-1)]
;		with U(0)=X(0) and U(N[U]-1)=X(N[X]-1)
;	v	[OUTPUT] ordinates, V(U)
;		* the approximating segment is (u[i],v[i]) -> (u[i+1],v[i+1])
;	w	[OUTPUT] ordinates if the approximation may be discontinuous
;		* the approximating segment is (u[i],w[i]) -> (u[i+1],v[i+1])
;
;keywords
;	noband	[INPUT] indicates whether the approximation is to be
;		restricted to the 'tolerance band' about the data
;		(the 'tolerance band' is a piecewise linear band
;		centered at the data whose width is determined by the
;		tolerances at the data points.)
;		* if set, indicates no band restriction
;		* if not set, or set to 0, indicates apply this restriction
;	discon	[INPUT] indicates whether or not the approximation must be
;		continuous
;		* if set, indicates continuity not required
;		* if not set, or set to 0, indicates continuity is required
;	verbose	[INPUT] controls chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;examples
;	x=findgen(100) & y=sqrt(x)
;	splac,x,y,0.1,u,v,w,/noband,verbose=10
;	splac,x,y,0.3,u,v,w,/discon,verbose=10
;	splac,x,y,-0.1,u,v,w & plot,x,y & oplot,u,v,psym=-1,col=100
;
;subroutines
;	NONE
;
;restrictions
;	Will not work with IDL versions < 5
;
;these are all the smoothing tools in PINTofALE
;	ALALOESS() makes a loess curve
;	CLTSMOOTH() removes -ves from background subtracted spectra or light curves
;	CONV_RMF convolves with specified RMF
;	HAARTRAN() smoothes by filtering on Haar wavelet coefficients
;	LINEREM() removes lines from spectrum by iterative S/N filtering
;	NOISMOOTH() does boxcar accumulation a la Ebeling's asmooth
;	REGROUP() accumulates from one end
;	SCRMF does fast convolution using optimized RMF
;	SMOOTHIE does peak-descent accumulation to build up S/N
;	SPLAC computes a segmented piecewise linear approximations
;	UNKINK() removes sharp kinks from a curve
;	VARSMOOTH() does local smoothing at specified scales
;	VOORSMOOTH() does fast local smoothing at specified scales
;
;history
;	From D.G. Wilson, 1976, ACM Transactions on Mathematical Software,2,388
;	http://www.netlib.org/toms/510
;	  converted to IDL by Vinay Kashyap, who is very unhappy about all
;	  those GOTOs, but doesn't have the patience to rewrite from scratch
;	  (FebMM)
;	bug correction for when YTOL was -ve; made into IDL5 (VK; JunMM)
;	bug correction in plotting with /DISCON (VK; AprMMVII)
;-

;	usage
ok='ok'
np=n_params() & nx=n_elements(x) & ny=n_elements(y) & ntol=n_elements(ytol)
if np lt 4 then ok='Insufficient parameters' else $
 if nx eq 0 then ok='Input data points undefined' else $
  if ny eq 0 then ok='Input data values undefined' else $
   if ntol eq 0 then ok='Tolerance undefined' else $
    if nx ne ny then ok='Input data mismatch' else $
     if ntol gt 1 and ntol ne nx then ok=$
	'Tolerance is a virtue sadly lacking here'
if ok ne 'ok' then begin
  print,'Usage: splac,x,y,ytol,u,v,w,/noband,/discon,verbose=verbose'
  print,'  compute a Segmented Piecewise Linear Approximation to Curve'
  if np ne 0 then message,ok,/info
  return
endif

;	check inputs
; there must be at least 3 points in curve
if nx lt 3 then begin
  message,'Nothing to approximate here, is there?',/info
  u=x & v=y & w=y & return
endif
; X must be in ascending order
dx=x[1:*]-x & oo=where(dx lt 0,moo)
if moo gt 0 then begin
  message,'Input data points not in ascending order; doing nothing.',/info
  u=x & v=y & w=y & return
endif
; YTOL must be +ve, and if not is taken to be fractional deviations
toly=0.*y
if ntol eq 1 then begin
  if ytol[0] lt 0 then toly=abs(ytol[0]*y) else toly[*]=ytol[0]
endif else begin
  oo=where(ytol lt 0,moo)
  if moo gt 0 then toly[oo]=abs(ytol[oo]*y[oo])
  oo=where(ytol ge 0,moo)
  if moo gt 0 then toly[oo]=ytol[oo]
endelse

;	keywords
if not keyword_set(noband) then noband=0	;noband=(0,1)::i2=(2,1)
if not keyword_set(discon) then discon=0	;discon=(0,1)::i3=(3,1)
i2=2 & if keyword_set(noband) then i2=1
i3=3 & if keyword_set(discon) then i3=1
itch=i2*i3
vv=0 & if keyword_set(verbose) then vv=fix(verbose) > 1

;	outputs
u=x & v=y & w=y

;	INITIALIZATION FOR ENTIRE PROGRAM
form="($,a)"
epsln=toly[1-1l] & sgn=1.0 & keep=1L & i=1L & u[1-1l]=x[1-1l] & j=2L
init=1L & indc=0L
goto, L30					;(30->

;	INITIALIZATION FOR EACH SEGMENT
L20: if vv ge 5 then print,form=form,'-20'		;(20<-

j=j+1L & init=i & indc=0L
if itch lt 3 then keep=i	;DISCON=1
if abs(itch-4) ne 2 then $	;NOBAND=1
	goto, L30				;(30->

;	RESTRICTED TO TOLERANCE BAND
xeye=u[j-1L-1l] & yeye=v[j-1L-1l]
temp1=epsln+(sgn*toly[i-1L-1l]-epsln)*(x[i-1l]-u[j-1L-1l])/(x[i-1l]-x[i-1L-1l])
yinit=yeye-temp1-temp1
goto, L40					;(40->

L30: if vv ge 5 then print,form=form,'-30'		;->30))
if vv ge 1 and vv lt 5 then print,form=form,'.'

;	NOT RESTRICTED TO TOLERANCE BAND
xeye=x[i-1l] & yeye=y[i-1l]+epsln & yinit=y[i-1l]-epsln
if itch eq 1 or i eq 1 then $	;NOBAND=DISCON=1
	goto, L40				;(40->
temp1 = sgn*toly[i+1L-1l]
smin = (y[i+1L-1l]-yeye-temp1)/(x[i+1L-1l]-xeye)
temp1 = sgn*toly[i-1L-1l]
smax = (yeye-y[i-1L-1l]+temp1)/(xeye-x[i-1L-1l])
if keep eq i-1L then goto, L50			;(50->
it=i-2L & xinit=xeye & ipiv=i & igraze=i & i=i+1L
goto, L150					;(150->

L40: if vv ge 5 then print,form=form,'-40'		;->40))

if xeye ge x[i-1l] then i=i+1L
epsln = sgn*toly[i-1l]
dx = x[i-1l]-xeye
smax = (y[i-1l]+epsln-yeye)/dx
smin = (y[i-1l]-epsln-yeye)/dx

L50: if vv ge 5 then print,form=form,'-50'		;->50)

xinit=xeye & ipiv=i & igraze=i

;	DETERMINATION OF INDIVIDUAL SEGMENT
L60: if vv ge 5 then print,form=form,'-60'		;((60<-
if vv eq 4 then print,form=form,'.'

if i eq nx then goto, L260			;(260->
i = i+1L

L70: if vv ge 5 then print,form=form,'-70'		;((70<-

;	TEST FOR NEW *MAX* SLOPE
dx=x[i-1l]-xeye & epsln=sgn*toly[i-1l]
temp1=(y[i-1l]+epsln-yeye)/dx
test = temp1-smax
if sgn le 0 then test=-test
if test eq 0 then goto, L90			;(90->
if test gt 0 then goto, L100			;(100->

;	TEST FOR END OF CANDIDATE SEGMENT
test = temp1-smin
if sgn le 0 then test=-test
if test lt 0 then goto, L210			;(210->
smax = temp1

L90: if vv ge 5 then print,form=form,'-90'		;->90)

;	TEST FOR NEW *MIN* SLOPE
ipiv = i

L100: if vv ge 5 then print,form=form,'-100'		;->100)

temp2 = (y[i-1l]-epsln-yeye)/dx
test = temp2 - smax
if sgn le 0 then test=-test
if test lt 0 then goto, L110			;(110->
if test eq 0 then goto, L120			;(120->
if test gt 0 then goto, L140			;(140->
L110: if vv ge 5 then print,form=form,'-110'		;->110)
test=smin-temp2
if sgn le 0 then test=-test
if test lt 0 then goto, L120			;(120->
if test eq 0 then goto, L130			;(130->
if test gt 0 then goto, L60			;->60)
L120: if vv ge 5 then print,form=form,'-120'		;->120))
smin=temp2
L130: if vv ge 5 then print,form=form,'-130'		;->130))
igraze=i
goto, L60					;->60)

;	CHECK FOR PIVOT AT NEW EYE POINT
L140: if vv ge 5 then print,form=form,'-140'		;->140)

if xeye eq x[ipiv-1l] then goto, L220		;(220->
epsln = sgn*toly[ipiv-1l]
indc=1
svx=xeye & svy=yeye & svmn=smin & svmx=smax
xeye=x[ipiv-1l] & yeye=y[ipiv-1l]+epsln
smin=smax & smax=(yinit-yeye)/(xinit-xeye)
if keep ge ipiv then goto, L170			;(170->
it = ipiv-1L

L150: if vv ge 5 then print,form=form,'-150'		;->150)

temp2=yeye+epsln
for l=keep,it do begin				;{160->
  temp2=yeye+sgn*toly[l-1l]
  temp1=(y[l-1l]-temp2)/(x[l-1l]-xeye)
  test=temp1-smax
  if sgn le 0 then test=-test
  if test lt 0 then smax=temp1
endfor						;L=KEEP,IT ->160}

L170: if vv ge 5 then print,form=form,'-170'		;->170)

if ipiv ge i-1L then goto, L70			;70->)
it=i-2L
temp2=yeye-epsln
idiot=ipiv
for l=idiot,it do begin				;{200->
  dx=x[l+1L-1l]-xeye
  temp2=yeye-sgn*toly[l+1L-1l]
  temp1=(y[l+1L-1l]-temp2)/dx
  test=temp1-smax
  if sgn le 0 then test=-test
  if test lt 0 then smax=temp1			;(->180<-)
  if test le 0 then ipiv=l+1L			;(->190<-)
endfor						;L=IDIOT,IT ->200}

goto, L70					;70->)

;	END OF CURRENT SEGMENT
L210: if vv ge 5 then print,form=form,'-210'		;->210)
temp2=smin
if i eq nx then goto, L240			;(240->
keep=igraze
goto, L250					;(250->

L220: if vv ge 5 then print,form=form,'-220'		;->220)

temp2=smax
if i eq nx then goto, L230			;(230->
sgn=-sgn & epsln=-epsln & keep=ipiv
goto, L250					;(250->

L230: if vv ge 5 then print,form=form,'-230'		;->230)

if indc eq 0 or xeye ne x[nx-2L-1l] then goto, L240;(240->
xeye=svx & yeye=svy & smin=svmn & smax=svmx

L240: if vv ge 5 then print,form=form,'-240'		;->240))

u[j-1l]=x[nx-1L-1l] & yinit=y[nx-1L-1L]
goto, L270					;(270->

L250: if vv ge 5 then print,form=form,'-250'		;->250))

if abs(itch-4) ne 2 then $	;NOBAND=1
	goto, L300				;(300->

;	DETERMINE KNOT ON EDGE OF TOLERANCE BAND
temp1=0.0
temp1=epsln-sgn*toly[i-1L-1l]
temp1=(y[i-1l]-y[i-1L-1l]+temp1)/(x[i-1l]-x[i-1L-1l])
u[j-1l] = (y[i-1l]+epsln-yeye-temp1*x[i-1l]+temp2*xeye)/(temp2-temp1)
goto, L310					;(310->

L260: if vv ge 5 then print,form=form,'-260'		;->260)

u[j-1l]=x[nx-1l] & yinit=y[nx-1l]

L270: if vv ge 5 then print,form=form,'-270'		;->270)

;	CONTINUITY CHECK FOR LAST SEGMENT
if itch ge 3 or init eq 0 then $	;DISCON=0
	goto, L290				;(290->
it=init-1L & svmx=smax+sgn & temp2=yeye+epsln
for l=kp,it do begin				;{280->
  temp2=yeye+sgn*toly[l-1l]
  temp1=(y[l-1l]-temp2)/(x[l-1l]-xeye)
  test=temp1-svmx
  if sgn le 0 then test=-test
  if test lt 0 then svmx=temp1
endfor						;L=KP,IT ->280}

if abs(svmx-smax+svmx-smin) le abs(smax-smin) then smax=svmx

L290: if vv ge 5 then print,form=form,'-290'		;->290)

;	NEARNESS CHECK FOR LAST SEGMENT
temp2=smax & temp1=yeye+smax*(u[j-1l]-xeye) & test=yinit-temp1
if sgn lt 0 then test=-test
if test gt 0 then goto, L310			;(310->
temp2=smin & temp1=yeye+smin*(u[j-1l]-xeye) & test=yinit-temp1
if sgn lt 0 then test=-test
if test gt 0 then goto, L310			;(310->
temp2=(yinit-yeye)/(u[j-1l]-xeye) & v[j-1l]=yinit
goto, L320					;(320->

L300: if vv ge 5 then print,form=form,'-300'		;->300)

if itch ge 3 then $	;DISCON=0
	goto, L330				;(330->
u[j-1l]=0.5*(x[i-1l]+x[i-1L-1l])
L310: if vv ge 5 then print,form=form,'-310'		;->310)))
v[j-1l]=yeye+temp2*(u[j-1l]-xeye)
L320: if vv ge 5 then print,form=form,'-320'		;->320)

if xeye ne xinit then goto, L330		;(330->
if itch eq 2 then $	;NOBAND=0 && DISCON=1
	goto, L360				;(360->
if itch ne 6 then $	;NOBAND=1 || DISCON=1
	goto, L330				;(330->
if j le 2 then goto, L380			;(380->
goto, L390					;(390->

L330: if vv ge 5 then print,form=form,'-330'		;->330)))

;	RECOMPUTATION OF KNOT FOR CONTINUITY
if j le 2 then goto, L370			;(370->
if slope eq temp2 then goto, L360		;(360->
yinit = v[j-2L-1l]
if itch lt 3 then yinit=w[j-2L-1l]	;DISCON=1
temp1=(xeye*temp2-u[j-2L-1l]*slope+yinit-yeye)/(temp2-slope)
if itch ge 3 then $	;DISCON=0
	goto, L350				;(350->
if temp1 gt xinit then goto, L360		;(360->
test=abs(epsln) & idiot=init-kp
for l=1,idiot do begin				;{340->
  it=init-l
  if temp1 ge x[it-1l] then goto, L350		;(350->
  dx=y[it-1l]-yeye-temp2*(x[it-1l]-xeye)
  test=toly[it-1l]
  if abs(dx) gt test then goto, L360		;(360->

endfor						;L=1,IDIOT ->340}

L350: if vv ge 5 then print,form=form,'-350'		;->350))

u[j-1L-1l]=temp1 & v[j-1L-1l]=yeye+temp2*(u[j-1L-1l]-xeye)
if itch lt 3 then w[j-1L-1l]=v[j-1L-1l]	;DISCON=1
goto, L390					;(390->

L360: if vv ge 5 then print,form=form,'-360'		;(->360))))

w[j-1L-1l]=yeye+temp2*(u[j-1L-1l]-xeye)
goto, L390					;(390->

L370: if vv ge 5 then print,form=form,'-370'		;->370)

if itch lt 3 then $	;DISCON=1
	goto, L360				;<-360)

L380: if vv ge 5 then print,form=form,'-380'		;->380)

v[1-1l]=yeye+temp2*(u[1-1l]-xeye)

L390: if vv ge 5 then print,form=form,'-390'		;->390)))

slope=temp2 & kp=keep
if i lt nx then goto, L20			;<-20)
if x[nx-1l] eq u[j-1l] then goto, L400		;(400->
if itch lt 3 then w[j-1l]=v[j-1l]	;DISCON=1
j=j+1L & u[j-1l]=x[nx-1l] & v[j-1l]=y[nx-1l]

L400: if vv ge 5 then print,form=form,'-400'		;->400)
if vv ge 1 and vv lt 5 then print,form=form,'.'

if j ge 2 and itch lt 3 then v[1-1l]=w[1-1l]	;DISCON=1
k=j-1L

u=u[0L:k] & v=v[0L:k] & w=w[0L:k-1l]

tt='Compression: '+strtrim(nx,2)+' -> '+strtrim(n_elements(u),2)
if vv ge 1 then message,tt,/info
if vv gt 0 then print,''
if vv ge 3 then begin
  if vv ge 10 then begin
    print,"plot,x,y,/xs,/ys,xtitle='X',ytitle='Y',title='"+tt+"',line=1,thick=3"
    if not keyword_set(discon) then print,'oplot,u,v,psym=-1' else $
    	print,'for i=0L,k do oplot,[u[i],u[i+1L]],[w[i],v[i+1L]],psym=-1'
  endif
  plot,x,y,/xs,/ys,xtitle='X',ytitle='Y',title=tt,line=1,thick=3
  if not keyword_set(discon) then oplot,u,v,psym=-1 else $
  	for i=0L,k-1L do oplot,[u[i],u[i+1L]],[w[i],v[i+1L]],psym=-1
endif
if vv ge 100 then stop,'HALTING.  type .CON to continue'

return
end