Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/PINTofALE/pro/setcont.pro Дата изменения: Mon Aug 11 07:02:22 2003 Дата индексирования: Tue Oct 2 00:26:53 2012 Кодировка:

function setcont,x,y,ysig,ycont=ycont,xcont=xcont,const=const,$
xrange=xrange,yrange=yrange,lock=lock,ccol=ccol, _extra=e
;+
;function setcont
; interactively determine a continuum
;
;syntax
; cont=setcont(x,y,ysig,ycont=ycont,xcont=xcont,/const,$
; xrange=xrange,yrange=yrange,ccol=ccol,/lock)
;
;parameters
; x [INPUT; required] abscissa values
; * must be an array with at least 2 elements
; y [INPUT; required] Y(X)
; * length must match X
; ysig [OUTPUT] errors on returned continuum (see description below)
;
;keywords
; ycont [I/O] if set on input, begins work from here
; * overwritten by new values on output
; xcont [I/O] if N(YCONT) .NE. N(X), then look here for the
; actual abscissae values.. YCONT=YCONT(XCONT)
; * if N(XCONT) .NE. N(YCONT), both are ignored
; const [I/O] if set, returns piecewise >constant< continum,
; as opposed to piecewise >linear< curve
; * can be toggled by clicking outside the plot area
; xrange [INPUT] xrange of plot -- passed with no comment to PLOT
; yrange [INPUT] yrange of plot -- passed with no comment to PLOT
; ccol [INPUT] color of continuum line
; * default is 150*!D.N_COLORS/256
; lock [INPUT] if set, locks up the relative positions of YCONT
; * has an effect ONLY if YCONT is set
; * side-effects:
; -- does not increase number of N({XY}CONT)
; -- left or middle click moves the entire curve up and down
; _extra [JUNK] here only to prevent crashing the program
;
;error calculation
; since the continuum is being set interactively, there is no meaningful
; way to define a statistical error on it.
;
;restrictions
; works in X. only?
;
;history
; vinay kashyap (Oct98)
; increased default plotting YRANGE (VK; Nov98)
; added keyword LOCK (VK; Dec98)
; cursor changes shape when mouse button is pressed (VK; FebMM)
; allowed return of zero errors, cleaned up error handling, added CCOL
; cleaned up LOCK (VK; MarMM)
; undid the cursor shape changes (VK; AprMM)
; bug: wasn't holding the poisson v/s zero error toggle (VK; JulMM)
; improved color-scale setting for 24-bit consoles (VK; FebMMI)
; changed default to produce errors=0 rather than Poisson errors
; (VK; Aug03)
;-

; usage
nx=n_elements(x) & ny=n_elements(y)
if nx lt 2 or nx ne ny then begin
print,'Usage: cont=setcont(x,y,ysig,ycont=ycont,xcont=xcont,/const,$'
print,' xrange=xrange,yrange=yrange,ccol=ccol,/lock)'
print,' return interactively determined continuum'
if nx eq 0 then return,[-1L] else return,0*x
endif

; sundry initializations
xmin=min(x,max=xmax) & ymin=min(y,max=ymax)
draglim=1e-4*sqrt((xmax-xmin)^2+(ymax-ymin)^2)
dncolors=256. > !D.N_COLORS ;24-bit color screen temporary fix
if not keyword_set(ccol) then ccol=fix(150.*float(!d.n_colors)/dncolors+1)
nc=0L

; is there an initial guess?
nyc=n_elements(ycont) & nxc=n_elements(xcont)
if nyc eq nx then begin & xcont=x & nxc=nx & endif
if nyc ne 0 and nyc eq nxc then begin
nc=nyc & xc=xcont & yc=ycont & yce=0*(sqrt(abs(yc)+0.75)+1.)
endif

; keyword check
if n_elements(xrange) lt 2 then xrange=[xmin,xmax]
if n_elements(yrange) lt 2 then begin
dy=ymax-ymin
yrange=[ymin-0.1*dy,ymax+0.1*dy]
endif
if not keyword_set(const) then const=0

; are the continuum points in lockstep?
lockstep=0
if keyword_set(lock) then begin
if nc gt 0 then lockstep=1
endif

; until otherwise stated, go and mark points for the continuum
go_on=1
if not keyword_set(lockstep) then begin
print,'LEFT CLICK to mark, LEFT DRAG to move'
print,'MIDDLE CLICK/DRAG to delete, RIGHT CLICK/DRAG to exit'
print,'CLICK ABOVE plot area to toggle piecewise constant curve'
print,'CLICK BELOW plot area to toggle Poisson errors v/s zero errors'
endif else begin
print,'LEFT to move, MIDDLE to undo, RIGHT to exit'
endelse

while go_on eq 1 do begin ;{endless loop

; plot
plot,x,y,/xs,/ys,psym=10,xrange=xrange,yrange=yrange
; reset xmin,xmax,ymin,ymax from the plot
xmin=min(!X.CRANGE) & xmax=max(!X.CRANGE)
ymin=min(!Y.CRANGE) & ymax=max(!Y.CRANGE)

; is there a continuum to overplot?
if lockstep then oplot,xc,yc,col=ccol else begin
nxc=n_elements(xc) & nyc=n_elements(yc) & nyce=n_elements(yce)
if nyce eq 0 and nyc ne 0 then yce=0*yc
if nxc eq 0 then begin & ycc=0*y & ycce=0*y & endif
if nxc eq 1 then begin & ycc=0*y+yc(0) & ycce=0*y+yce(0) & endif
if nxc gt 1 then begin
if keyword_set(const) then begin
xcc=fltarr(2L*nc) & ycc=xcc & ycce=xcc
dxc=xc(1:*)-xc
for i=0,nc-1 do begin
if i eq 0 then xcc(2*i)=xmin else xcc(2*i)=xc(i)-0.5*dxc(i-1)
if i lt nc-1 then xcc(2*i+1)=xc(i)+0.5*dxc(i) else xcc(2*i+1)=xmax
ycc([2*i,2*i+1])=yc(i)
ycce([2*i,2*i+1])=yce(i)
endfor
ycc=interpol(ycc,xcc,x)
ycce=interpol(ycce,xcc,x)
endif else begin
ycc=interpol(yc,xc,x) & ycce=interpol(yce,xc,x)
endelse
endif
oplot,x,ycc,col=ccol
oplot,x,ycc+ycce,col=ccol,line=1 & oplot,x,ycc-ycce,col=ccol,line=1
if nc gt 0 then oplot,xc,yc,psym=1
endelse

; mark
cursor,x0,y0,/down,/data
mbutton=!mouse.button
if !D.NAME eq 'X' then begin
;if mbutton eq 1 then device,cursor_standard=74
;if mbutton eq 2 then device,cursor_standard=82
;if mbutton eq 4 then device,cursor_standard=100
endif
cursor,x1,y1,/up,/data
;if !D.NAME eq 'X' then device,/cursor_original
drag=sqrt((x1-x0)^2+(y1-y0)^2)
if drag lt draglim then drag=0.
y0e=0*(sqrt(abs(y0)+0.75)+1.)
y1e=0*(sqrt(abs(y1)+0.75)+1.)

outclik=0 ;click was inside plot area?
if x0 gt xmax then outclik=outclik+2^0 ;right of plot area
if y0 gt ymax then outclik=outclik+2^1 ;above plot area
if x0 lt xmin then outclik=outclik+2^2 ;left of plot area
if y0 lt ymin then outclik=outclik+2^3 ;below plot area
;
if keyword_set(lockstep) then outclik=0
if outclik ne 0 then begin
mbutton=0 ;mark/delete no points
; toggle piecewise constant if click is above plot area
if outclik ge 2 and outclik le 6 then begin
const=1-const
c1='changing from piecewise '
if const eq 0 then c1=c1+'constant to piecewise linear' else $
c1=c1+'linear to piecewise constant' & message,c1,/info
endif
; toggle poisson errors v/s zero errors if click is below plot area
if outclik ge 8 then begin
if n_elements(yc) gt 0 then begin
c1='changing errors from '
if n_elements(yce) eq 0 then yce=0*yc
if yce(0) eq 0 then begin
message,c1+'zero to Poisson',/info
yce=sqrt(abs(yc)+0.75)+1.
endif else begin
message,c1+'Poisson to zero',/info
yce(*)=0.
endelse
endif else message,'Nothing to make an error out of yet!',/info
endif
endif

if mbutton eq 1 then begin ;(LEFT
if drag eq 0 then begin ;(click
; mark this point
if nc eq 0 then begin
xc=[x0] & yc=[y0] & yce=[y0e]
endif else begin
if lockstep then begin ;(just move the whole curve
;find nearest point to X0 and peg it to Y0
tmp=min(abs(xc-x0),ic) & dely=y0-yc(ic)
yc=yc+dely
endif else begin ;)(else add a new point
xc=[xc,x0] & yc=[yc,y0]
if yce(0) ne 0 then yce=[yce,y0e] else yce=[yce,0]
ox=sort(xc) & xc=xc(ox) & yc=yc(ox) & yce=yce(ox)
endelse ;lockstep?)
endelse
endif else begin ;)(click+drag
; find nearest point and move it to new position
if nc gt 0 then begin
if lockstep then begin ;(move curve as a whole
tmp=min(abs(xc-x0),ic) & dely=y1-yc(ic)
yc=yc+dely
endif else begin ;)(move just this point
dd=sqrt((xc-x0)^2+(yc-y0)^2) & tmp=min(dd,ic)
xc(ic)=x1 & yc(ic)=y1 & if yce(0) ne 0 then yce(ic)=y1e
ox=sort(xc) & xc=xc(ox) & yc=yc(ox) & yce=yce(ox)
endelse ;lockstep?)
endif
endelse ;drag)
nc=n_elements(xc)
endif ;LEFT)

if mbutton eq 2 then begin ;(MIDDLE
if lockstep then begin ;(undo previous shift
if not keyword_set(dely) then dely=0.
yc=yc-dely
dely=-dely
endif else begin ;)(delete nearest point
; delete nearest point
if nc gt 0 then begin
dd=sqrt((xc-x0)^2+(yc-y0)^2) & tmp=min(dd,ic)
ii=lindgen(nc) & oi=where(ii ne ic,moi)
if moi gt 0 then begin
xc=xc(oi) & yc=yc(oi) & yce=yce(oi) & nc=moi
endif else nc=0L
endif
endelse ;lockstep?)
endif ;MIDDLE)

if mbutton eq 4 then begin ;(RIGHT
go_on=0 ;exit
endif ;RIGHT)

endwhile ;end of endless loop}

; outputs
nxc=n_elements(xc) & nyc=n_elements(yc)
if nxc gt 0 and nxc eq nyc then begin
xcont=xc & ycont=yc & ysig=yce
endif
;
if nxc eq 0 then begin & cont=0*y & ysig=0*y & endif
if nxc eq 1 then begin & cont=0*y+yc(0) & ysig=0*y+yce(0) & endif
if nxc gt 1 then begin ;(interpolate YC onto X
if keyword_set(const) then begin
xcc=fltarr(2L*nc) & cont=xcc & ysig=xcc
dxc=xc(1:*)-xc
for i=0,nc-1 do begin
if i eq 0 then xcc(2*i)=xmin else xcc(2*i)=xc(i)-0.5*dxc(i-1)
if i lt nc-1 then xcc(2*i+1)=xc(i)+0.5*dxc(i) else xcc(2*i+1)=xmax
cont([2*i,2*i+1])=yc(i)
ysig([2*i,2*i+1])=yce(i)
endfor
cont=interpol(cont,xcc,x)
ysig=interpol(ysig,xcc,x)
endif else begin
cont=interpol(yc,xc,x) & ysig=interpol(yce,xc,x)
endelse
endif ;YC(XC))

;; error-bars
;if ysig(0) ne 0 then ysig=sqrt(abs(cont)+0.75)+1.

; plot
plot,x,y,/xs,/ys,psym=10,xrange=xrange,yrange=yrange
oplot,x,cont,col=ccol
oplot,x,cont+ysig,col=ccol,line=1 & oplot,x,cont-ysig,col=ccol,line=1
if nc gt 0 then oplot,xc,yc,psym=1

return,cont
end