PRO k_w1,vel_xyt,xoy,dt,kw_dia,nux,nut ; Calculo del 2-D diagrama k_w a partir de una serie temporal. ; Cambiar mediante editor el parametro "apo" segun se quiera (=1) ; o no (=0) apodizacion. Also change optionally the % of tapered ; pixels by the apodization y la escala espatial. ; ; From the 3-D data volume VEL_XYT (x,y,t) is produced a 2-D (!) ; k-w diagram corresponding to a chosen coordinate (x or y), which ; is an average of individual k-w diagrams computed for different ; values of the other coordinate. ; ;Inputs: ; ; vel_xyt = 3-D (x,y,t) imagen a estudiar. ; ; xoy = X or Y coordinate option. ; X coord ... xoy le 0 (e.g. 0) ; Y coord ... xoy gt 0 (e.g. 1) ; ; dt = intervalo de tiempo entre dos observaciones (intervalo de ; muestreo temporal) en secundos (v sekundach). ; ;Outputs: ; ; kw_dia = diagrama k-w calculado ; ; nux,nut = arrays of Fourier coordinates ; on_error,1 s=size(vel_xyt) ; 1-x,2-y,3-t if s(0) ne 3 then message,'Input must be a 3-D array' if xoy gt 0 then begin ; Y-axis chosen sx=s(2) nn=s(1) endif else begin ; X-axis chosen sx=s(1) nn=s(2) endelse st=s(3) nyqt=fix(st/2) ; Nyquist frequencies nyqx=fix(sx/2) dx=0.062 ; !!!!!! SCALE (arcsec) !!!!!!!!! ; Determinacion de los valores de muestreo en el dominio de Fourier nux=findgen(nyqx+1)/(sx*dx) nut=findgen(nyqt+1)/(st*dt) print,'Intervalos de muestreo: ',nux(1),' arcsec^-1',$ nut(1),' Hz' print,'Frecuencias de Nyquist: ',nux(nyqx),' arcsec^-1',$ nut(nyqt),' Hz' print,'Averaged over:',fix(nn),' time-slices' ; Apodization mask (apma) perc=10. ; !!!!!!!!!% of the smoothed edges !!!!!!!!! npx=fix(round(sx*perc/100.)) npt=fix(round(st*perc/100.)) u=fltarr(sx)+1. v=fltarr(st)+1. u(0)=0. v(0)=0. for i=1,npx do begin u(i)=(1.-cos(!pi*i/npx))/2. endfor for j=1,npt do begin v(j)=(1.-cos(!pi*j/npt))/2. endfor u(sx-npx:sx-1)=reverse(u(1:npx)) v(st-npt:st-1)=reverse(v(1:npt)) apma=u#v ; Loop over individual time-slices kw_dia=0 ; summation variable for i=0,nn-1 do begin print,i if xoy gt 0 then vel=reform(vel_xyt(i,*,*)) $ ; Y-axis else vel=reform(vel_xyt(*,i,*)) ; X-axis ; Correccion de trend temporal. meancol=total(vel,1)/float(sx) ;avg. over space coef=poly_fit(findgen(st),meancol,2,yfit) mask=(fltarr(sx)+1)#yfit ;mask to correct for trend in t dimension vel=vel-mask ; Apodization apo=1 ;set apo to zero if no apodization is desired if (apo ne 0) then begin mea=mean(vel) vel=vel-mea ;subtracting the mean value vel=vel*apma endif ; Calcula la transformada de Fourier (FFT) kw=abs(fft(vel,-1)) ; amplitude spectrum kw_dia=kw_dia+kw(0:nyqx,0:nyqt)^2 ; power summation endfor kw_dia=kw_dia/nn nux=nux(0:nyqx) nut=nut(0:nyqt) END