function findscale,curve,dim,crunch=crunch,half=half,pick=pick,choice=choice,\$ eps=eps,_extra=e ;+ ;function findscale ; returns the lengthscale in pixels at each point on the given curve ; ;syntax ; ls=findscale(curve,dim,/crunch,/half,pick=pick,choice=choice,eps=eps) ; ;parameters ; curve [INPUT; required] regularly gridded array of function ; values to be used to compute the length scales ; * if scalar, returns 0 ; * if 2D, ; -- use DIM to specify primary dimension ; -- compute lengthscales separately along each projection ; * if >2D, convert to 1D ; dim [INPUT; default=1] primary dimension in case of 2D array ; (e.g., if CURVE=CURVE(NX,NY), DIM=2 returns SCALE=SCALE(NY)) ; ;keywords ; crunch [INPUT] if set, and CURVE is 2D, collapses the array along ; the secondary dimension to generate 1D curve ; half [INPUT] if set, returns the half-scale ; pick [INPUT; default=0] if 2D, specifies how to combine the ; scales computed at the different cuts ; 0: pick the smallest scale ; 1: pick the largest scale ; 2: get the average ; choice [INPUT; default=0] what algorithm to use to find the scale? ; 0: MexicanHat wavelet ; 1: use inverse of 1st derivative ; 2: radius of curvature ; 3: stepped toggle ; eps [INPUT; default=1e-7] small number ; _extra [JUNK] ignore. here only to prevent crashing program. ; ;subroutines ; WVLT_SCALE [ROOFN] ; ;history ; vinay kashyap (Apr97) ; added CHOICE option 3 (VK; Feb03) ;- ; usage if n_elements(curve) eq 0 then begin print,'Usage: ls=findscale(curve,dim,/crunch,/half,pick=pick,choice=choice,eps=eps)' print,' returns length scales at each point along curve' return,0L endif ; save inputs f=curve & if keyword_set(dim) then d=fix(dim) else d=1 ; check dimensions szf=size(f) & nszf=n_elements(szf) if szf(0) eq 0 then return,[0L] ;scalar -- return 0 if szf(0) gt 2 then begin ;convert to 1D f=[temporary(f(*))] & szf=size(f) & nszf=n_elements(szf) endif if szf(0) ne 2 then d=1 ;only 1 D, see? nx=szf(1) & if szf(0) eq 1 then ny=1L else ny=szf(2) ; if primary dimension is the 2nd D, then transpose the matrix if d eq 2 then begin nx=szf(2) & ny=szf(1) & f=transpose(temporary(f)) endif ; catch trivial errors if nx lt 3 then return,lonarr(NX) ;scalar masquerading as array ; collapse to 1D if szf(0) eq 2 and keyword_set(crunch) then begin g=reform(f(*,0)) for ix=0,nx-1 do g(ix)=total(f(ix,*)) f=g & ny=1 endif ; initialize if not keyword_set(choice) then choice=0 ;how to make the scales if not keyword_set(pick) then pick=0 ;how to combine scales across dimensions scale=lonarr(nx)+nx & if pick eq 2 then scale(*)=0 ;the output if not keyword_set(eps) then eps=1e-7 ;"epsilon" norm=lonarr(nx) ;for PICK=2 ; get length scale scl=lonarr(nx) for iy=0,ny-1 do begin ;{shtep through secondary dimensions g=reform(f(*,iy)) ;la function dg=deriv(g) ;derivative d2g=deriv(dg) ;2nd derivative gb=intarr(nx)+1 ;coverage function ok=where(g lt 0.01*max(g),mok) & if mok gt 0 then gb(ok)=0 if not keyword_set(choice) then begin scl=wvlt_scale(g,_extra=e) endif else begin if choice(0) eq 1 then scl=ceil(abs(g)/(abs(dg)>eps)) else \$ if choice(0) eq 2 then scl=ceil((1.+dg^2)^(1.5)/(abs(d2g)>eps)) else \$ if choice(0) eq 3 then scl=scl+gb else \$ scl=wvlt_scale(g,_extra=e) endelse if choice(0) ne 3 then begin for ix=0,nx-1 do begin ;if 2D, we gotta pick s=scale(ix) if g(ix) gt eps*max(f) then begin if pick eq 0 then scale(ix)=scl(ix) < s if pick eq 1 then scale(ix)=scl(ix) > s if pick eq 2 then begin scale(ix)=scl(ix) + s norm(ix)=norm(ix)+1L endif endif endfor endif else scale=max(scl)-scl+1L endfor ;IY=0,NY-1} norm=norm>1 & if pick eq 2 then scale=ceil(float(scale)/float(norm)) ; return the half-scale if asked if keyword_set(half) then scale=scale/2 return,scale end