function wvlt_scale,curve,coarse=coarse,_extra=e
;+
;function	findscale
;	returns the lengthscale in pixels at each point on the given curve
;	using the MexHat Wavelet Transforms as a guide
;
;syntax
;	ls=wvlt_scale(curve,/coarse)
;
;parameters
;	curve	[INPUT; required] regularly gridded 1D array of function
;		values to be used to compute the length scales
;		* if scalar, returns 0
;		* if >1D, convert to 1D
;
;keywords
;	coarse	[INPUT; default=2] explore scales at specified coarseness
;		e.g.: if COARSE=1, every scale from 1 to N(CURVE)-1 will
;		be tested
;	_extra	[INPUT] junk; ignore.  here only to prevent crashing program.
;
;history
;	vinay kashyap (Apr97)
;	fixed bug when N(CURVE) was a power of 2 (LiWei Lin/VK; Jan'03)
;-

;	usage
if n_elements(curve) eq 0 then begin
  print,'Usage: ls=wvlt_scale(curve,/coarse)'
  print,'  returns length scales at each point along curve'
  return,0L
endif

;	check inputs
f=curve & szf=size(f)
if szf(0) eq 0 then return,[0L]		;scalar -- return 0
if szf(0) gt 1 then begin		;convert to 1D
  f=[temporary(f(*))] & szf=size(f) & nszf=n_elements(szf)
endif
;
if not keyword_set(coarse) then coarse=2 else coarse=fix(coarse)>1

;	initialize
n=szf(1)				;number of elements in array
fmax=abs(f)				;array of maximum values of f
nn=roofn(n,2)				;change to power of 2 (for FFT)
n2=nn*2L				;double that, for padding (FFT)
x=findgen(n2)				;for defining the MexHat
if nn gt n then ff=[f,f(0:nn-n-1)] else ff=f	;fill out the reminder (periodic bc)...
g=[reverse(ff(0:nn/2-1)),ff,reverse(ff(nn/2:*))]	;...re-position...
f=g					;...rename
scale=lonarr(n)+long(n/sqrt(2))		;the output
	;the sqrt(2) factor because the maximum response of a MexHat to
	;a Gaussian of sigma=d occurs at a scale of d/sqrt(2)

;	catch trivial errors
if n eq 1 then return,[0L]		;scalar masquerading as an array

;	convolve with MexHat wavelet at different scales
for i=1,n-1,coarse do begin		;{shtep through scales
  scal=float(i)						;sigma of the MexHat
  g=(x-n2/2)/scal & g=g^2 & mh=(1.-g) & g=g<138
  norm=3./sqrt(2*scal)/!pi^(0.25)
  mh=mh*exp(-g/2) & mh=shift(mh,n2/2)*norm		;the MexHat
  ff=float(fft(fft(f,1)*fft(mh,1),-1))			;convolve
  ff=ff(n2/4:n2/4+n-1)					;unpad
  oo=where(abs(fmax) lt abs(ff),moo)
  if moo gt 0 then begin				;stronger signal?
    scale(oo)=long(scal/sqrt(2)+0.5) & fmax(oo)=abs(ff(oo))
  endif
endfor					;I=0,N-1,COARSE}

return,scale
end