; file: plotpowermap.pro = plot Fourier power map ; init: Aug 24 2010 ; last: Aug 30 2010 ; note: needs ssw ; ;+ ; NAME: plotpowermap ; ; PURPOSE: plot a power map ; Format described in Krijger et al. 2001A&A...379.1052K ; ; CALLING SEQUENCE: ; plotpowermap,cube,cadence,fmin,fmax,$ ; olddata=olddata,mode=mode,normalization=normalization,$ ; apod=apod, pxdetrend=pxdetrend, meandetrend=meandetrend,$ ; _extra=frimagekeywords ; ; INPUTS: ; cube: (x,y,t) data cube, any type ; cadence: [regular] sampling cadence in sec ; fmin,fmax: frequency band in Hz ; optional keywords: ; olddata=1: re-use the previous cube (its FT is a file on /tmp) ; mode=0: plot log(power) (default) ; mode=1: plot linear power ; mode=2: plot sqrt(power) = amplitude ; normalization=0: no power normalization (default) ; normalization=1: divided by power at zero frequency (modulation) ; normalization=2: divided by amplitude at zero frequency (Leahy) ; apod: fractional extent of apodization edges; default 0.1 ; pxdetrend=1: subtract linear trend with time per pixel (a bit slow) ; pxdetrend=2: subtract linear trend with time per pixel (better but slower) ; meandetrend: subtract linear trend with time for the image means ; _extra=frimagekeywords: type IDL> frimage ; ; MODIFICATION HISTORY ; ; Aug 2010: Rob Rutten (RR) assembly of Alfred de Wijn's routines. ;- ;----------------------------------------------------------------------- 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 by ssw's mpfitfun in apodtempcube ymod=p[0] + x * p[1] ;RR p[i] are mpfitfun's return parameters return, ymod end ;----------------------------------------------------------------------- function apodtempcube,cube,apod,meandetrend,pxdetrend ; apodizes the temporal (x,y,*) columns of an (x,y,t) data cube ; optional detrending, either mean image sequence or per (x,y) pixel sizecube=size(cube) nx=sizecube[1] ny=sizecube[2] nt=sizecube[3] apocube=cube tf=findgen(nt) + 1 col=fltarr(nt) apodt = fltarr(nt)+1 if (apod ne 0) then begin apodrim = apod*nt apodt[0] = (sin(!pi/2.*findgen(apodrim)/apodrim))^2 apodt = apodt*shift(rotate(apodt,2),1) ;RR had ik nooit verzonnen endif ; meandetrend: get spatially-averaged trend fit=0 if (meandetrend) then begin avgseq=fltarr(nt) for it=0,nt-1 do avgseq[it]=total(cube[*,*,it]) avgseq=avgseq/(double(nx)*double(ny)) meanfitp = mpfitfun('linear',tf,avgseq,fltarr(nt)+1,[1000.,0.],/quiet) ;RR AdW; ssw; Levenberg-Marquardt least-squares fit to IDL function ;RR fltarr(nt)+1 = 1 sigma error = set to 1; what units? meanfit=meanfitp[0]+tf*meanfitp[1]-total(avgseq)/double(nt) endif ; apodize per [x,y] temporal column = time sequence per px for ix=long(0),long(nx)-1 do begin for iy=long(0),long(ny)-1 do begin col=cube[ix,iy,*] meancol=avgstd(col) if (meandetrend) then col=col-meanfit if (pxdetrend eq 1) then begin pxfitp=poly_fit(tf,col,1) col=col-pxfitp[0]-tf*pxfitp[1]+meancol endif if (pxdetrend eq 2) then begin ;RR slow pxfitp = mpfitfun('linear',tf,col,fltarr(nt)+1,[meancol, 0.],/quiet) col=col-pxfitp[0]-tf*pxfitp[1]+meancol endif apocube[ix,iy,*]=(col-meancol)*apodt+meancol endfor if (pxdetrend ne 0) then $ writeu,-1,string(format='(%"\r == detrend next row... ",i5,"/",i5)',$ ix,nx) endfor return,apocube end ;----------------------------------------------------------------------- function getftcube,cube,olddata ;RR Aug 24 2010 make and write, or read, ftcube ;RR after detrending set olddata to use the /tmp/ftcube again sizecube=size(cube) nx=sizecube[1] ny=sizecube[2] nt=sizecube[3] ftcube=complexarr(nx,ny,nt/2+1) if (olddata eq 0) then begin for ix=0, nx-1 do begin for iy=0, ny-1 do begin col=cube[ix,iy,*] ftcube[ix,iy,*]=(fft(col,-1))[0:nt/2] endfor endfor openw, unit, '/tmp/ftcube', /get_lun writeu, unit, ftcube close, unit free_lun, unit endif else begin openr, unit, '/tmp/ftcube', /get_lun readu, unit, ftcube close, unit free_lun, unit endelse return,ftcube end ;------------------------------------------------------------------------- function getpowermap,ftcube,cadence,fmin,fmax,normalization ;RR cadence in sec, min and max in Hz ;RR normalization: none if 0, other values see below sizecube=size(ftcube) nx=sizecube[1] ny=sizecube[2] nf=sizecube[3] nt=2*nf-2 frequencies=1./(cadence*2)*findgen(nt/2+1)/(nt/2) sf=where(frequencies ge fmin and frequencies le fmax) powermap=fltarr(nx,ny) slice=complexarr(nx,ny) for i=0, n_elements(sf)-1 do begin slice[*,*]=ftcube[*,*,sf[i]] powermap=powermap+abs(slice)^2 endfor slice[*,*]=ftcube[*,*,0] if (normalization eq 0) then print,' normalization = 0: none' if (normalization eq 1) then begin powermap=powermap/abs(slice)^2 print,' normalization=1: divide by power at f=0' endif if (normalization eq 2) then begin powermap=powermap/abs(slice) ; Leahy print,' normalization=2: divide by amplitude at f=0 (Leahy)' endif if (normalization eq 3) then begin powermap=powermap/abs(slice)^(8./9.) ; Alfred?? print,' normalization=3: funny' endif return, powermap end ; ====================== MAIN ROUTINE ================================ pro plotpowermap,cube,cadence,fmin,fmax,$ outmap=outmap,$ olddata=olddata,mode=mode,normalization=normalization,$ apod=apod,pxdetrend=pxdetrend,meandetrend=meandetrend,$ _extra=frimagekeywords ;RR Aug 24 2010 make and write, or read, ftcube ;RR detrending is slow, set olddata to re-use /tmp/ftcube if not keyword_set(outmap) then mapout=0 else mapout=1 if not keyword_set(olddata) then olddata=0 if not keyword_set(normalization) then normalization=0 if not keyword_set(mode) then mode=0 if (n_elements(apod) ne 0) then apod=apod else apod=0.1 if (n_elements(pxdetrend) ne 0) then pxdetrend=pxdetrend else pxdetrend=0 if not keyword_set(meandetrend) then meandetrend=0 if n_params() lt 4 then begin print,' plotpowermap, cube, cadence, fmin, fmax, outmap=outmap,$' print,' olddata=olddata, mode=mode, normalization=normalization,$' print,' apod=apod, pxdetrend=pxdetrend, meandetrend=meandetrend,$' print,' _extra=frimagekeywords' return endif ; apodize the cube apocube=apodtempcube(cube,apod,meandetrend,pxdetrend) ; make or read Fourier transform cube ftcube=getftcube(apocube,olddata) ; get powermap powermap=getpowermap(ftcube,cadence,fmin,fmax,normalization) ; plot power map, with optional frimage keywords if (mode eq 0) then begin powermap=alog10(powermap) print,' mode = 0: log(power)' endif if (mode eq 1) then print,' mode=1: linear power' if (mode eq 2) then begin powermap=sqrt(powermap) print,' mode = 2: sqrt(power)' endif frimage,bytscl(powermap),_extra=frimagekeywords ; output powermap if (mapout eq 1) then outmap=powermap end