; file: plotkopower.pro = plot "k-omega" = k-f power diagram ; init: Sep 3 2010 ; last: Aug 8 2018 Rob Rutten Deil ; note: test per Hyper-C at bottom ; add mean power spectra along sides as Fig 3 2003A&A...407..735R ;+ ; NAME: plotkopower ; ; PURPOSE: plot "k-omega" diagram = Fourier power at temporal frequency ; f against horizontal spatial wavenumber k_h ; ; CALLING SEQUENCE: ; plotkopower, cube, arcsecpx, cadence, plotfilename, $ ; apod=apod, $ ; kmax=kmax, fmax=fmax, minpower=minpower, maxpower=maxpower, $ ; contours=contours, lamb=lamb, fundamental=fundamental ; ; INPUTS: ; cube: (x,y,t) data cube, any type ; arcsecpx = angular image scale in arcsec/px ; cadence: [regular] sampling cadence in sec ; plotfilename: postscript output file name ; optional keywords: ; apod: fractional extent of apodization edges; default 0.1 ; kmax: maximum k_h axis as fraction of Nyquist value, default 0.2 ; fmax: maximum f axis as fraction of Nyquist value, default 0.5 ; minpower: minimum of power plot range, default maxpower-5 ; maxpower: maximum of power plot range, default alog10(max(power)) ; contours: set true to plot contours ; lamb: value > 0 overplots Lamb line omega = c_s kn at c_s = lamb km/s ; fundamental: set true to overplot fundamental mode omega=sqrt(g k_h) ; ; MODIFICATION HISTORY: ; Sep 2010: Rob Rutten (RR) assembly of Alfred de Wijn's routines ; Oct 2010: RR optional Lamb line and fundamental mode ; Mar 2011: RR -1 > +1 line 164 from Alfred de Wijn ;- ;----------------------------------------------------------------------- function avgstd, array, stdev=stdev ; get array average + standard deviation ;RR Aug 23 2010 found in Mabula Haverkamp's IDL, later also AdW ;RR not the same as avg.pro in ssw ;RR not the same as avg.pro in Pit Suetterlin DOT software ;RR so renamed to avgstd avrg = total(array)/n_elements(array) stdev = sqrt(total((array-avrg)^2)/n_elements(array)) return, avrg end ;---------------------------------------------------------------------------- function linear, x, p ;RR used in temporal detrending ymod = p[0] + x * p[1] return, ymod end ;---------------------------------------------------------------------------- function gradient, x, y, p ;RR used in spatial detrending zmod = p[0] + x * p[1] + y * p[2] return, zmod end ;---------------------------------------------------------------------------- function apod3dcube,cube,apod ; apodizes cube in all three coordinates, with detrending ; get cube dimensions sizecube=size(cube) nx=sizecube[1] ny=sizecube[2] nt=sizecube[3] apocube = fltarr(nx,ny,nt) ; define temporal apodization apodt = fltarr(nt)+1 if (apod ne 0) then begin apodrimt = nt*apod apodt[0] = (sin(!pi/2.*findgen(apodrimt)/apodrimt))^2 apodt = apodt*shift(rotate(apodt,2),1) ;RR had ik nooit verzonnen endif ; temporal detrending, not per px, only mean-image trend ttrend = fltarr(nt) tf = findgen(nt) + 1 for it=0, nt-1 do begin img = cube[*,*,it] ttrend[it] = avgstd(img) endfor fitp = mpfitfun('linear', tf, ttrend, fltarr(nt)+1, [1000.,0.],/quiet) fit = fitp[0] + tf * fitp[1] ; temporal apodization per (x,y) column ;RR do not reinsert trend to keep [0,0] Fourier pixel from dominating for it=0, nt-1 do begin img = cube[*,*,it] apocube[*,*,it] = (img-fit[it])*apodt[it] ;RR + ttrend[it] endfor ; define spatial apodization apodx = fltarr(nx)+1 apody = fltarr(ny)+1 if (apod ne 0) then begin apodrimx=apod*nx apodrimy=apod*ny apodx[0] = (sin(!pi/2.*findgen(apodrimx)/apodrimx))^2 apody[0] = (sin(!pi/2.*findgen(apodrimy)/apodrimy))^2 apodx = apodx*shift(rotate(apodx,2),1) apody = apody*shift(rotate(apody,2),1) apodxy = apodx # apody endif ; spatial gradient removal + apodizing per image xf = fltarr(nx,ny)+1. yf = xf for it=0, nt-1 do begin img = apocube[*,*,it] avg = avgstd(img) ;RR mpfit2dfun = ssw/gen/idl/fitting/mpfit/mpfit2dfun.pro fitp = mpfit2dfun('gradient',xf,yf,img,fltarr(nx,ny)+1,[1000.,0.,0.],$ /quiet) fit = fitp[0]+xf*fitp[1]+yf*fitp[2] apocube[*,*,it] = (img-fit)*apodxy + avg endfor ; done return,apocube end ;--------------------------------------------------------------------------- function ko_dist, sx, sy, double=double ; set up Pythagoras distance array from origin ;RR from Alfred de Wijn email Aug 30 2010 dx = rebin(dindgen(sx/2+1)/(sx/2),sx/2+1,sy/2+1) dy = rotate(rebin(dindgen(sy/2+1)/(sy/2),sy/2+1,sx/2+1),1) dxy = sqrt(dx^2+dy^2)*(min([sx,sy])/2+1.) distances = dblarr(sx,sy) distances[0,0] = dxy ; get other quadrants distances[sx/2,0] = rotate(dxy[1:*,*],5) ;RR 5 = 90 deg distances[0,sy/2] = rotate(dxy[*,1:*],7) ;RR 7 = 270 deg distances[sx/2,sy/2] = rotate(dxy[1:*,1:*],2) ;RR 2 = 180 deg if not keyword_set(double) then distances = fix(round(distances)) return, distances end ;--------------------------------------------------------------------------- function averpower,cube ; compute 2D (k_h,f) power array by circular averaging over k_x, k_y ; get cube dimensions sizecube=size(cube) nx=sizecube[1] ny=sizecube[2] nt=sizecube[3] ; forward fft and throw away half of it ; perform fft in time direction first fftcube = (fft(fft((fft(cube,-1, dimension=3))[*,*,0:nt/2],-1,$ dimension=1),dimension=2)) ; set up distances fftfmt = size(fftcube) distances = ko_dist(nx,ny) ;RR integer-rounded Pythagoras array ;; maxdist = min([nx,ny])/2-1 ;RR largest quarter circle ;RR corrected Alfred de Wijn email answering Julius Koza Jan 19 2012 maxdist = min([nx,ny])/2+1 ;RR largest quarter circle ; get average power over all k_h distances, building power(k_h,f) avpow = fltarr(maxdist+1,nt/2+1) for i=0, maxdist do begin waar = where(distances eq i) for j=0, nt/2 do begin w1 = (fftcube[*,*,j])[waar] avpow[i,j] = total(w1*conj(w1))/n_elements(waar) endfor ;; writeu,-1,string(format='(%"\rcomputing avpow... ",i6,"/",i6)',i,maxdist) endfor ; done return, avpow end ;--------------------------------------------------------------------------- pro koplotpow,avpow,arcsecpx,cadence,kmax,fmax,minpower,maxpower,contours,$ lamb,fundamental,plotfilename ; plotting program sizepow = size(avpow) ; select extent = fractions of Nyquist values plotrange = [fix(sizepow[1]*kmax),fix(sizepow[2]*fmax)] plotpow = avpow[0:plotrange[0]-1,0:plotrange[1]-1] ;RR 5x5 resizing, I guess for better tick positioning and contours xas = 2.*!pi/(arcsecpx*2)*findgen(plotrange[0])/(sizepow[1]-1) rexas = 2.*!pi/(arcsecpx*2)*findgen(plotrange[0]*5)/(sizepow[1]*5-1) yas = 1./(cadence*2)*findgen(plotrange[1])/(sizepow[2]-1)*1e3 reyas = 1./(cadence*2)*findgen(plotrange[1]*5)/(sizepow[2]*5-1)*1e3 plotpowrebin = convol(rebin(plotpow,plotrange[0]*5,plotrange[1]*5),$ fltarr(6,6)+1/(6.*6.),/edge_truncate) xrange = [min(xas),max(xas)] yrange = [min(yas),max(yas)] ; plot startup: what windows? windowsdevice='X' if (!d.name eq 'WIN') then windowsdevice='WIN' ;RR inserted Alfred's startplot parameters instead of calling startplot set_plot, 'ps' !p.font=1 !p.thick=3 !x.thick=3 !y.thick=3 !z.thick=3 aspect=3./2 device,filename=plotfilename,xsize=7,ysize=7,$ set_font='Times', /tt_font,font_size=10,bits_per_pixel=8 ; plot power image logarithmically tv, bytscl(alog10(plotpowrebin) > minpower < maxpower),$ 0.15, 0.15, /normal, xsize=0.7, ysize=0.48 ; define plotframe for overplots plot, xrange, yrange, /nodata, xstyle=13, ystyle=13,$ position=[0.15,0.15,0.85,0.63], /noerase ; overplot contours if requested ; RR contours are white; full range OK since only max power is white if (contours ne 0) then begin contour,alog10(plotpowrebin)-minpower, rexas, reyas, $ xstyle=9, ystyle=9, /overplot, levels=[0,1,2,3,4,5,6,7], $ c_labels=[1,1,1,1,1,1,1,1], c_thick=[2],charsize=0.65, color=255 contour, alog10(plotpowrebin)- minpower, rexas, reyas, $ xstyle=9, ystyle=9, /overplot, levels=[1,3,5,7,9,11,13]/2., $ c_labels=[0], c_thick=[1], color=255 endif ; overplot Lamb line (horizontal propagation at Lamb value = sound speed) ; 725 is km/arcsec on the sun at mean distance from Earth ; 1/2pi from omega to f; 1E3 from Hz to mHz if (lamb ne 0) then oplot, rexas, 1./(2.*!pi)*1E3*lamb/725.*rexas, color=255 ; overplot fundamental mode ; g = 2.74E4 cm/s^2 = 2.74E-1 km/s^2 = 2.741E-1/725. arcsec/s^2 if (fundamental ne 0) then $ oplot, rexas, 1./(2.*!pi)*1E3*sqrt(2.741E-1/725.*rexas),$ color=255, linestyle= 2 ; plot x,y axes plot, xrange, yrange, /nodata, xstyle=9, ystyle=9,$ position=[0.15,0.15,0.85,0.63], $ yticklen=-0.015/0.7, xticklen=-0.015/0.53, /noerase, $ xminor=1, xtitle='horizontal wavenumber [arcsec!U-1!N]', $ yminor=1, ytitle='frequency [mHz]' ; plot wavelength axis along top if (xrange[1] lt 10) then begin ;RR I wonder how to automate this wavtickspos = [10, 5, 3, 2, 1.5, 1] wavticksn = ['10','5','3','2','1.5','1'] endif else if (xrange[1] lt 20) then begin wavtickspos = [10, 5, 3, 2, 1.5, 1, 0.5] wavticksn = ['10','5','3','2','1.5','1','0.5'] endif else if (xrange[1] lt 50) then begin wavtickspos = [5.0, 2.0, 1.0, 0.5, 0.2] wavticksn = ['5','2','1','0.5','0.2'] endif else begin wavtickspos = [5.0, 2.0, 1.0, 0.5, 0.2, 0.1, 0.05] wavticksn = ['5','2','1','0.5','0.2','0.1','0.05'] endelse wavticks=n_elements(wavtickspos)-1 wavticksv = 2.*!pi/wavtickspos ;RR wavelength wav num = circle frequency ;; tickpos = [25,5,3,2,1.5,1,0.8,0.7,0.6] ;; ticknum = n_elements(tickpos)-1 ;; ticknames = ['25','5','3','2','1.5','1','0.8','0.7','0.6'] ;; tickvals = 2.*!pi/tickpos axis, /xaxis, xticks=wavticks, xtickv=wavticksv, xtickname=wavticksn, $ ticklen=-0.015/0.53, xminor=1, xtitle='wavelength [arcsec]' ; plot period axis along righthand side if (yrange[1] lt 10) then begin ;RR I wonder how to automate this pertickspos = [20, 10, 5, 3, 2, 1] perticksn = ['20','10','5','3','2','1'] endif else if (yrange[1] lt 20) then begin pertickspos = [10, 5, 3, 2, 1.5, 1, 0.5] perticksn = ['10', '5','3','2','1.5','1','0.5'] endif else if (yrange[1] lt 50) then begin pertickspos = [10, 5, 2, 1, 0.5, 0.2, 0.1] perticksn = ['10','5','2','1','0.5','0.2','0.1'] endif else if (yrange[1] lt 100) then begin pertickspos = [2, 1, 0.5, 0.2, 0.1] perticksn = ['2','1','0.5','0.2','0.1'] endif else begin pertickspos = [0.5, 0.2, 0.1, 0.05, 0.02] perticksn = ['0.5','0.2','0.1','0.05','0.02'] endelse perticks=n_elements(pertickspos)-1 perticksv = 1e3/pertickspos/60. ;RR period, from mHz to min axis,/yaxis,yticks=perticks, ytickv=perticksv, ytickname=perticksn, $ ticklen=-0.015/0.7, yminor=1, ytitle='period [min]' ; add grayscale bar on top grayscale = bytscl(findgen(256,1)) tv, grayscale, 0.15, 0.91, /normal, xsize=0.7, ysize=0.02 xtickv=findgen(maxpower-minpower+1) plot, [0,maxpower-minpower], [0,1], /noerase, /nodata, xstyle=1, ystyle=1, $ position=[0.15,0.91,0.85,0.93],yticklen=0.0,xticklen=-0.015/0.02, $ xminor=1, xtitle='log power [arbitrary]',$ yminor=1,yticks=1,ytickname=[' ',' '] ; done device, /close set_plot, windowsdevice ; convert the IDL ad so that gv shows the filename as window label if (windowsdevice eq 'X') then begin spawn,'cat '+plotfilename+$ '| sed "s|Graphics produced by IDL|'+plotfilename+$ '|" > idltemp.ps; mv idltemp.ps '+plotfilename endif end ; --------------------------- main part ------------------------------ pro plotkopower,cube,arcsecpx,cadence,plotfilename,$ apod=apod,kmax=kmax,fmax=fmax,minpower=minpower,maxpower=maxpower,$ contours=contours,lamb=lamb,fundamental=fundamental ; wrapper calling the above subroutines if (n_elements(apod) ne 0) then apod=apod else apod=0.1 if (n_elements(kmax) ne 0) then kmax=kmax else kmax=1.0 if (n_elements(fmax) ne 0) then fmax=fmax else fmax=1.0 if not keyword_set(contours) then contours=0 else contours=1 if (n_elements(lamb) ne 0) then lamb=lamb else lamb=7. if not keyword_set(fundamental) then fundamental=0 else fundamental=1 if (kmax gt 1) then kmax=1 if (fmax gt 1) then fmax=1 if n_params() lt 4 then begin print, ' plotkopower,cube,arcsecpx,cadence,plotfilename,' print, ' apod=apod (default 0.1),' print, ' kmax=kmax (default 1.0), print, ' fmax=fmax (default 1.0), print, ' minpower=minpower (default maxpower-5),' print, ' maxpower=maxpower (default alog10(max(cube))' print, ' contours=contours (default 0)' print, ' lamb=lamb (sound speed km/s, default 7)' print, ' fundamental=fundamental (default 0)' return endif ; apodize the cube apocube=apod3dcube(cube,apod) ; compute radially-averaged power avpow=averpower(apocube) ; set min, max cutoffs maxavpow=alog10(max(avpow)) minavpow=alog10(min(avpow)) print, ' max log(power) = ', maxavpow, ' min log(power) = ',minavpow if (n_elements(maxpower) ne 0) then maxpower=maxpower else maxpower=maxavpow if (n_elements(minpower) ne 0) then minpower=minpower else minpower=maxpower-5 ; plot the diagram koplotpow,avpow,arcsecpx,cadence,kmax,fmax,minpower,maxpower,contours,$ lamb,fundamental,plotfilename ; done print,' wrote file ', plotfilename end ; -------------------------- test per Hyper-C -------------------------- cd,'/home/rutten/rr/edu/manuals/idl-cubes/idl-cube' gb=readfits('DOT-2005-10-14-gb.fits') ca=readfits('DOT-2005-10-14-ca.fits') plotkopower,ca,0.071,30,'/tmp/ko-power-ca.ps',maxpower=1,kmax=0.3,/contours,$ lamb=7,/fundamental ; inspect with: $gv /tmp/ko-power-ca.ps end