; ****************** ; * * ; * SPDYNPS.PRO * ; * * ; ****************** ;-------------------------------------------------------------------- pro FOND, tab, fon,sigma ;-------------------------------------------------------------------- ; Background and sigma of a 1-D distribution of intensities. ; tab = (INPUT) 1-D array of intensities ; fon,sigma = (OUTPUT) background and fluctuations (1 sigma level) if n_elements(tab) gt 1 then sigma=10.*stdev(tab,fon) else begin $ fon=tab(0) & sigma=0 endelse encore: test=where(abs(tab-fon) lt 2.5*sigma) if n_elements(test) eq 1 then return sigma=stdev(tab(test),moy) if moy eq fon then return fon=moy goto,encore end ;-------------------------------------------------------------------- pro FAIT_FOND, tab,fichier, f,s ;-------------------------------------------------------------------- ; Calculation of Background and Sigma arrays of a 2-D distribution of ; intensities, frequency by frequency. ONLY values >min(tab) are used. ; A save set can be put in 'fichier'. ; tab = (INPUT) 2-D array of intensities, function of (time,freq) ; fichier = (INPUT) name of save_set containing f & s ; f,s = (OUTPUT) arrays of background and fluctuations (1 sigma level) nf=n_elements(tab(0,*)) f=fltarr(nf) & s=fltarr(nf) mintab=min(tab) for i=0,nf-1 do begin ; if (i mod 50) eq 0 then print,i test=where(tab(*,i) gt mintab) if test(0) ne -1 then begin FOND,tab(test,i),ff,ss f(i)=ff & s(i)=ss endif endfor if fichier ne '' then save,f,s,filename=fichier return end ;-------------------------------------------------------------------- pro RETIRE_FOND, f,s,seuil, tab ;-------------------------------------------------------------------- ; Substraction of (background + seuil * sigmas) from a 2-D array of raw ; data, and setting of resulting negative values to 0. ; f,s = (INPUT) arrays of background and fluctuations (1 sigma level) ; seuil = (INPUT) threshold value for the number of sigmas above background ; tab = (INPUT/OUTPUT) 2-D array of raw data (time,freq) -> tab-(f+seuil*s) for i=0,n_elements(tab(*,0))-1 do tab(i,*)=(tab(i,*)-f-seuil*s) >0 return end ;-------------------------------------------------------------------- pro DYN, tab,fracmin,fracmax, tabmin,tabmax ;-------------------------------------------------------------------- ; Adjustment of dynamic range tabmax=-1 & tabmin=-1 test=where(tab gt min(tab)) if test(0) eq -1 then begin print,'Null dynamic range' return endif test=float(tab(test)) dh=max(test)-min(test) if dh eq 0 then begin print,'Null dynamic range > min(tab)' tabmin=min(tab) & tabmax=max(tab) return endif dh=dh/1000. h=histogram(float(test),min=min(test),max=max(test),binsize=dh) xh=findgen(n_elements(h))*dh+min(test) th=total(h) & nh=0 for i=0,n_elements(h)-1 do begin nh=nh+h(i) if nh le fracmin*th then tabmin=i if nh le fracmax*th then tabmax=i endfor if tabmin eq -1 then tabmin=0 tabmin=xh(tabmin) if tabmax eq -1 then tabmax=0 tabmax=xh(tabmax) return end ;-------------------------------------------------------------------- pro SPDYNPS, image,xmin,xmax,ymin,ymax,legende_x,legende_y,titre, $ grid,fond,seuil,fracmin,fracmax,color,bar_titre, log=log, screen=screen ;-------------------------------------------------------------------- ; PostScript / Screen plot of dynamic spectrum. ; image = (INPUT) dynamic spectrum ; xmin,xmax,ymin,ymax = (INPUT) ranges of x & y coordinates ; legende_x,legende_y,titre = (INPUT) x & y captions + title of plot ; grid = (INPUT) =1 to superimpose a grid to the plot, =0 otherwise ; fond = (INPUT) =1 to substract a background to the data before ; plotting, =0 otherwise ; seuil = (INPUT) =number of sigmas -in addition to fond- to substract ; from data ; fracmin,fracmax = (INPUT) min and max threshold for optimization of the ; dynamic range of the plot ; color = (INPUT) grey scale of plot: =0 -> black=intense, ; =1 -> white=intense ; bar_titre = (INPUT) title of dynamic bar (unit = dB, except if '.') ; Last parameters can be omitted (beginning anywhere after 'image') ; Default values = 0,n_spectra,0,n_frequencies,' ',' ',' ',0,0,0,.05,.95,0,' ' if n_elements(xmin) eq 0 or n_elements(xmax) eq 0 then begin xmin=0 & xmax=n_elements(image(*,0)) endif if n_elements(ymin) eq 0 or n_elements(ymax) eq 0 then begin ymin=0 & ymax=n_elements(image(0,*)) endif if n_elements(legende_x) eq 0 then legende_x=' ' if n_elements(legende_y) eq 0 then legende_y=' ' if n_elements(titre) eq 0 then titre=' ' if n_elements(grid) eq 0 then grid=0 if n_elements(fond) eq 0 then fond=0 if n_elements(seuil) eq 0 then seuil=0 if n_elements(fracmin) eq 0 or n_elements(fracmax) eq 0 then begin fracmin=0.05 & fracmax=0.95 endif if n_elements(color) eq 0 then color=0 if n_elements(bar_titre) eq 0 then bar_titre=' ' tab=image ; Fond if fond eq 1 then begin FAIT_FOND, tab,'', f,s RETIRE_FOND, f,s,seuil, tab endif ; Dynamic range if (fracmin ne 0.) or (fracmax ne 1.) then begin DYN, tab,fracmin,fracmax, tabmin,tabmax if (tabmax eq -1) or (tabmin eq tabmax) then RETURN tab=bytscl(tab,min=tabmin,max=tabmax) endif else begin test=where(tab gt min(tab)) tabmin=min(tab(test)) & tabmax=max(tab(test)) tab=bytscl(tab,min=tabmin,max=tabmax) endelse if color eq 0 then tab=255-tab ; 0->noir=intense, 1->blanc=intense if not keyword_set(screen) then print,'Writing in Postscript file ...' tl=-0.015 if grid eq 1 then tl=1.0 nx=!p.multi(1) & ny=!p.multi(2) if nx eq 0 then nx=1 if ny eq 0 then ny=1 x=!p.multi(0) if x ne 0 then x=nx*ny-x i=x mod nx j=ny-(x-i)/nx-1 dx=(0.9-0.1*nx)/nx & dy=(0.9-0.1*ny)/ny posxmin=0.08+i*(0.12+dx) & posymin=0.14+j*(0.1+dy) posxmax=posxmin+dx & posymax=posymin+dy cm=1. & if nx ge 3 or ny ge 3 then cm=1.5 if keyword_set(screen) then cm=1.4 if keyword_set(log) and not keyword_set(screen) then $ plot_io, [xmin,xmax], [ymin,ymax], /nodata, xra=[xmin,xmax], xstyle=13, $ yra=[ymin,ymax], ystyle=13, title=titre, ticklen=tl, charsize=1.1*cm, $ pos = [posxmin,posymin,posxmax,posymax] if keyword_set(log) and keyword_set(screen) then $ plot_io, [xmin,xmax], [ymin,ymax], /nodata, xra=[xmin,xmax], xstyle=13, $ yra=[ymin,ymax], ystyle=13, title=titre, ticklen=tl, charsize=1.4*cm, $ /device, position=[90,50,90+n_elements(image(*,0)),50+n_elements(image(0,*))] if not keyword_set(log) and not keyword_set(screen) then $ plot, [xmin,xmax], [ymin,ymax], /nodata, xra=[xmin,xmax], xstyle=13, $ yra=[ymin,ymax], ystyle=13, title=titre, ticklen=tl, charsize=1.1*cm, $ pos = [posxmin,posymin,posxmax,posymax] if not keyword_set(log) and keyword_set(screen) then $ plot, [xmin,xmax], [ymin,ymax], /nodata, xra=[xmin,xmax], xstyle=13, $ yra=[ymin,ymax], ystyle=13, title=titre, ticklen=tl, charsize=1.4*cm, $ /device, position=[90,50,90+n_elements(image(*,0)),50+n_elements(image(0,*))] tv, tab, !x.window(0), !y.window(0), xsize=!x.window(1)-!x.window(0), $ ysize=!y.window(1)-!y.window(0),/normal axis, xaxis=0, xra=[xmin,xmax], xstyle=1, xticklen=tl, $ xtitle=legende_x, charsize = 1.*cm axis, xaxis=1, xra=[xmin,xmax], xstyle=1, xticklen=tl, $ xtickname=replicate(' ',9) axis, yaxis=0, yra=[ymin,ymax], ystyle=1, yticklen=tl, $ ytitle=legende_y, charsize = 1.*cm axis, yaxis=1, yra=[ymin,ymax], ystyle=1, yticklen=tl, $ ytickname=replicate(' ',9) if not keyword_set(screen) then begin posxmin=posxmax+0.03*dx/0.8 & posymin=posymin+0.05*dy/0.8 posxmax=posxmin+0.03*dx/0.8 & posymax=posymax-0.05*dy/0.8 bar_tick = [' ','dB',' '] if bar_titre eq '.' then begin bar_titre=' ' & bar_tick = [' ',' ',' '] endif plot, [0,1], [0,1], /nodata, /noerase, xst=1, yst=12, $ pos = [posxmin,posymin,posxmax,posymax], $ /norm, xticklen=0, xticks=2, xtickv=[0,0.5,1], $ xtickname=bar_tick, xtitle=bar_titre, charsize=1.1*cm axis, /yaxis, yra=[tabmin,tabmax], ticklen=-0.02, /ysty, charsize=1.*cm ; dynamique=alog10(tabmax-tabmin) ; if dynamique ge 0. then dynamique=fix(dynamique) ; if dynamique lt 0. then dynamique=fix(dynamique)-1 ; dynamique=10.^dynamique ; ndyn = fix((tabmax-tabmin)/dynamique) ; axis, /yaxis, yra=[tabmin,tabmax], ticklen=-0.02, /ysty, yticks=ndyn-1, $ ; ytickv=tabmin+dynamique-(tabmin mod dynamique)+indgen(ndyn)*dynamique, $ ; yminor=-1 bcb = bytscl(replicate(1,10)#bindgen(256)) if color eq 0 then bcb=255-bcb tv, bcb, !x.window(0), !y.window(0), xsize = !x.window(1) - !x.window(0), $ ysize = !y.window(1) - !y.window(0), /normal oplot, [1,0,0,1], [1,1,0,0] endif RETURN END ;-------------------------------------------------------------------- pro HELP_SPDYNPS ;-------------------------------------------------------------------- print,'SPDYNPS, image, xmin,xmax,ymin,ymax, legende_x,legende_y,titre, $' print,' grid,fond,seuil,fracmin,fracmax,color,bar_titre, [/log], [/screen]' return end