function powerspec,fov,Nbin,tmap

;author:  Chris Vale, UC Berkeley, June 02 2003
;program to return a power spectrum for a 2-d input map
;returns L,L^2C_L/2pi for an input 2-d map.  return is of the form [2,*]
;where [0,*] is L, and [1,*] is L^2C_L/2pi.  values are returned for 
;as many as Nbin locations, but zero values are not included
print,'computing ell, ell squared CL on 2 pi'

d_el = 360.0/fov
N    = round(sqrt(n_elements(tmap)))
CL2d = (abs(fft(tmap,/double)))^2
kx   = (dist(N,1)#replicate(1,N))
ky   = (replicate(1,N)#dist(N,1))
k2   = (dist(N,1)#replicate(1,N))^2 + (replicate(1,N)#dist(N,1))^2
k2   = k2*2.0*!pi*2.0*!pi

lkmn = 2.0*alog(2.0*!pi)
lkmx = 2.0*alog(0.5*!pi*N)
k2mn = exp(lkmn)
k2mx = exp(lkmx)
lk2  = k2
lk2[where(k2 lt k2mn/2.0)] = 0.1
lk2  = alog(lk2)
ibin = round((lk2-lkmn)*Nbin/(lkmx-lkmn) + 0.499999)
if min(ibin) lt 0 then begin
 ibin[where(ibin le 0)] = 0
endif
ibin[where(ibin ge Nbin)] = Nbin-1
ibin[where(k2 le k2mn)] = -1
ibin[where(k2 ge k2mx)] = Nbin

;print,ibin(999:1023,1)

;print,k2(0:23,99)
;print
;print,CL2d(0:23,99)
;print
;print,k2(999:1023,1)
;print
;print,CL2d(999:1023,1)
;print
;print,size(where(ibin eq 1))
;print,size(where(ibin eq 2))
;print,size(where(ibin eq 3))
;print,size(where(ibin eq 5))
;print,size(where(ibin eq Nbin))

;print,mean(sqrt(k2(where(ibin eq 5))))

CL   = fltarr(Nbin+1)
ell  = fltarr(Nbin+1)
CL2d = CL2d*k2/(2.0*!pi)
for indx = 0,Nbin do begin
 wher = where(ibin eq indx)
 print,size(wher)
 if n_elements(wher) gt 1 then begin
  ell[indx] = mean(sqrt(k2[wher]))
  CL[indx]  = mean(CL2d[wher])
 endif
 if n_elements(wher) eq 1 then begin
  if wher ge 0 then begin
   ell[indx] = sqrt(k2[wher])
   CL[indx]  = CL2d[wher]
  endif
 endif
endfor


print,CL

wher = where(ell gt 0.1)
CL   = CL[wher]
ell  = ell[wher]*d_el/(2.0*!pi)

nel  = n_elements(ell) - 2
el_CL= fltarr(2,nel+1)
el_CL[0,*] = ell[0:nel]
el_CL[1,*] = CL[0:nel]

print,nel+1,' non zero bins'

return,el_CL

end