; file: scatter_aw.pro from Alfred de Wijn Oct 29 2012 ; last: Sep 19 2013 Rob Rutten, Deil ; note: contours are factor 2 apart ; level = set point density for outer contour ; _extra for oplot and contour parameters function scfunc, p ; determines intersection of moment curves common scattercommon, xmomval, xmomind, ymomind, ymomval xmvint = interpol(xmomval, xmomind, p[0]) ymvint = interpol(ymomval, ymomind, p[1]) return, sqrt((xmvint-p[1])^2+(ymvint-p[0])^2) end ;------------------------------------------------------------------------- ; OPLOT_DL.PRO v1.50 (26 August 2007) --- (c) 1999-2007 H.W. de Wijn ; Macro to draw dashed curve ;------------------------------------------------------------------------- ; Parameters: ; x,y = supporting points ; length_1 = length of dash in x-axis data units ; length_2 = length of blank in x-axis data units ; Keywords: ; /y_units -> use y-axis instead of x-axis data units ; /quiet -> quiet operation ; Other keywords are passed on to OPLOT, which draws the dashes. Keywords ; supported include: thick, clip, and noclip. Note that psym is ignored. ;------------------------------------------------------------------------- PRO OPLOT_DL,x,y,length_1,length_2,y_units=YS,quiet=QU,psym=p,_extra=e ; occurence of keyword psym ensures its absence in _extra=e ; psym is set to zero in calls of oplot below IF n_elements(QU) eq 0 THEN QU=0 ars_y=FLOAT((!y.crange[1]-!y.crange[0])/(!x.crange[1]-!x.crange[0])$ *(!d.x_size*(!x.window[1]-!x.window[0]))$ /(!d.y_size*(!y.window[1]-!y.window[0]))) ;aspect ratio in x,y data units ars_x=1.0 IF n_elements(YS) eq 0 THEN YS=0 IF YS EQ 1 THEN BEGIN ars_x=1.0/ars_y ars_y=1.0 ENDIF size_x=size(x) & i_max=size_x[1]-1 IF QU NE 1 THEN print,'OPLOT_DL: # of supporting points =',i_max+1 xa=x[0] & ya=y[0] ;starting values i=1l BEGIN_LOEP: L=0 & i=i-1 iaa=i & xaa=xa & yaa=ya ;starting values of a dash plus blank ;----- dash ---- REPEAT BEGIN ;loop until length_1 is overshot xx=xa & yy=ya ;starting values of step i=i+1 ;see for which i we surpass length_1 IF i gt i_max THEN goto,END_DASH xxx=x[i] & yyy=y[i] ;end values of step xa=xxx & ya=yyy ;starting values of next step DL=sqrt(((xxx-xx)/ars_x)^2+((yyy-yy)/ars_y)^2) L=L+DL ENDREP UNTIL (L gt length_1) xa=xxx-(L-length_1)/DL*(xxx-xx) ;shift back by excess part of last step ya=yyy-(L-length_1)/DL*(yyy-yy) IF iaa+1 le i-1 THEN $ oplot,[xaa,x[iaa+1:i-1],xa],[yaa,y[iaa+1:i-1],ya],psym=0,_extra=e $ ELSE $ oplot,[xaa,xa],[yaa,ya],psym=0,_extra=e ;----- blank ----- L=0 & i=i-1 REPEAT BEGIN ;loop until length_2 is overshot xx=xa & yy=ya ;starting values of step i=i+1 ;see for which i we surpass length_2 IF i gt i_max THEN goto,END_LOEP xxx=x[i] & yyy=y[i] ;end values of step xa=xxx & ya=yyy ;starting values of next step DL=sqrt(((xxx-xx)/ars_x)^2+((yyy-yy)/ars_y)^2) L=L+DL ENDREP UNTIL (L gt length_2) xa=xxx-(L-length_2)/DL*(xxx-xx) ;shift back by excess part of last step ya=yyy-(L-length_2)/DL*(yyy-yy) GOTO,BEGIN_LOEP ;----- last dash ----- END_DASH: oplot,[xaa,x[iaa+1:i_max]],[yaa,y[iaa+1:i_max]],psym=0,_extra=e END_LOEP: END ;=========================================================================== pro scatter_aw, xdata, ydata, xr, yr, pos, $ moments=moments, hists=hists, gray=gray, level=level, _extra=extra moments = keyword_set(moments) hists = keyword_set(hists) if strcmp(!d.name, 'ps', /fold_case) then begin realxsize = (!x.window[1]-!x.window[0]) * !d.x_size/!d.x_px_cm realysize = (!y.window[1]-!y.window[0]) * !d.y_size/!d.y_px_cm dx = !x.crange[1]-!x.crange[0] dy = !y.crange[1]-!y.crange[0] ; empirical facts. if not keyword_set(level) then level = 30 binwidth = 0.075 ; number of bins nx = fix(realxsize/binwidth + 0.5) ny = fix(realysize/binwidth + 0.5) ; binsize in x and y bx = float(dx)/(nx-1.) by = float(dy)/(ny-1.) ; scale data if !x.crange[0] ne 0 then xd = xdata - !x.crange[0] else xd = xdata if bx ne 1 then xd /= bx if !y.crange[0] ne 0 then yd = ydata - !y.crange[0] else yd = ydata if by ne 1 then yd /= by ; create histogram array h = long(nx) * fix(yd) + fix(xd) ; set points outside of !x.crange and !y.crange to -1 range = xd ge 0 and xd lt nx and yd ge 0 and yd lt ny h = (temporary(h) + 1) * temporary(range) - 1 ; compute histogram h = histogram(h, min=0, max=long(nx)*long(ny)-1l) h = reform(h, nx, ny, /overwrite) xas = findgen(nx) * bx + !x.crange[0] yas = findgen(ny) * by + !y.crange[0] nl = fix(alog10(max(h)/float(level))/alog10(2)) if keyword_set(gray) then begin tv, 255b-bytscl(h,min=0), pos[0], pos[1], xsize=pos[2]-pos[0], ysize=pos[3]-pos[1], /normal endif else begin densatpos = interpolate(h, temporary(xd), temporary(yd)) wltlevel = where(densatpos lt level) oplot, xdata[wltlevel], ydata[wltlevel], psym=3, _extra=extra contour, h, xas, yas, xstyle=13, ystyle=13, /overplot, $ levels=2^(findgen(nl+1))*level, _extra=extra endelse if moments or hists then begin rowtot = total(h,1) coltot = total(h,2) endif if moments then begin common scattercommon, xmomval, xmomind, ymomind, ymomval rowxtot = fltarr(ny) colytot = fltarr(nx) for i=0,ny-1 do rowxtot[i] = total(h[*,i]*xas) for i=0,nx-1 do colytot[i] = total(h[i,*]*yas) w = where(rowtot gt 0) xmomval = rowxtot[w]/rowtot[w] xmomind = yas[w] xmiavg = avg(xmomind) w = where(coltot gt 0) ymomval = colytot[w]/coltot[w] ymomind = xas[w] ymiavg = avg(ymomind) xi = transpose([[1.0,0.0],[0.0,1.0]]) r = [xmiavg,ymiavg] powell, r, xi, 1e-4, fmin, 'scfunc' xmvint = interpol(xmomval, xmomind, r[0], /spline) ymvint = interpol(ymomval, ymomind, r[1], /spline) wx = reverse(where(xmomind lt r[0])) wy = reverse(where(ymomind lt r[1])) oplot_dl, [xmvint,xmomval[wx]], [r[0],xmomind[wx]], $ (!x.crange[1]-!x.crange[0])/40., $ (!x.crange[1]-!x.crange[0])/30., /quiet, _extra=extra oplot_dl, [r[1],ymomind[wy]], [ymvint,ymomval[wy]], $ (!y.crange[1]-!y.crange[0])/40., $ (!y.crange[1]-!y.crange[0])/30., /y_units, /quiet, _extra=extra wx = where(xmomind gt r[0]) wy = where(ymomind gt r[1]) oplot_dl, [xmvint,xmomval[wx]], [r[0],xmomind[wx]], $ (xr[1]-xr[0])/40., (xr[1]-xr[0])/30., /quiet, _extra=extra oplot_dl, [r[1],ymomind[wy]], [ymvint,ymomval[wy]], $ (yr[1]-yr[0])/40., (yr[1]-yr[0])/30., /y_units, /quiet, $ _extra=extra endif if hists then begin oplot, xr[1]-rowtot/max(rowtot)*dx/5, yas, _extra=extra oplot, xas, yr[1]-coltot/max(coltot)*dy/5, _extra=extra endif endif else begin oplot, xdata, ydata, psym=3, _extra=extra endelse end