pro fft_filter_sin

!p.multi=0
Tx=1
Nx=100
x=findgen(Nx)/(Nx-1)*Tx
Ty=1
Ny=100
y=findgen(Ny)/(Ny-1)*Ty
T=fltarr(Nx,Ny)
for i=0,Nx-1 do begin
   for j=0,Ny-1 do begin
      T(i,j)= sin(5.*!pi*((30.*i/Nx)))*cos(2.*!pi*((30.*i/Nx)+(20.*j/Ny)))
      endfor
   endfor
;
!p.multi=[0,2,2]
shade_surf,T,x,y,xtitle='x',ytitle='y',title='T(x,y)'
;
;
C=fft(T,-1) ; complex fourier coefficients
surface,float(C)
aaa=where(float(c) gt .4)
surface,imaginary(C)
;
a=findgen(Nx/2+1)
b=-reverse(findgen(Nx/2-1)+1)
ns=[a,b] ; this is the array of n's associated with C(n). n goes with x
;print,ns ; n goes from 0,...,Nx/2, -(Nx/2-1),...,-1
subn=sort(ns) ;  n goes with x
n_sort=ns(subn)
;
a=findgen(Ny/2+1)
b=-reverse(findgen(Ny/2-1)+1)
ms=[a,b] ; this is the array of m's associated with C(n,m)
print,ms ; m goes from 0,...,Ny/2, -(Ny/2-1),...,-1
subm=sort(ms) ; m goes with y
m_sort=ms(subm)
;
sub_n_p3=where(ns eq 3)
sub_n_n3=where(ns eq -3)
sub_n_0=where(ns eq 0)
;
sub_m_p2=where(ms eq 2)
sub_m_n2=where(ms eq -2)
sub_m_0=where(ms eq 0)
;
print,'C(3,0),c(-3,0)=',C(sub_n_p3,sub_m_0),C(sub_n_n3,sub_m_0)
;
print,'C(0,2),c(0,-2)=',C(sub_n_0,sub_m_p2),C(sub_n_0,sub_m_n2)
;
; now we need to define CC(n,m) to have normal scaling in n & m.
;
CC=C*0.
;
for n=0,Nx-1 do begin
   for m=0,Ny-1 do begin
      CC(subn(n),subm(m))= C(n,m)
      endfor
   endfor
;
surface,ABS(CC),n_sort,m_sort
;
stop

end