pro spectra_2d_plot,spec,noblock=noblock,lev=lev,fac=fac,surf=surf,yr=yr,xr=xr,ymarg=ymarg,xmarg=xmarg,title=title,nocol=nocol,noyax=noyax,_extra=e

; Purpose: plot (contour) a 2-d fft spectrum
; Category: graphics statistics
; See-also: spectra_2d (which explains the structure of the input)
; Example:
;	a1 = fltarr(72,12*15) of seaice edge anomalies
;	spectra_2d_plot,spectra_2d(a1,yfac=12),xr=[-6,6],yr=[0,.5],lev=.05+indgen(10)/40.

@comm_error

; Input
if (datatype(spec) ne 8) then message,'Require structure generated by specta_2d as input'

; Keywords, misc
if (n_elements(fac) eq 0) then fac=50
if (not keyword_set(yr)) then yr=makerange(spec.f)
if (not keyword_set(xr)) then xr=makerange(spec.w)
if (not keyword_set(ymarg)) then ymarg=[8,8]
if (not keyword_set(xmarg)) then xmarg=[8,8]

; Compute contour levels
if (n_elements(lev) eq 0) then begin
  nl=10
  mx=max(spec.spec)
  lev=indgen(nl)*mx/nl+mx/float(fac)
endif else begin
  nl=n_elements(lev)
endelse

; Compute colours
  cols=(wmc_guess_col(2*nl+1))(nl:2*nl)

; Find the max value in the plot
mx=max(spec.spec,i)
s=size(spec.spec)
l1=i mod s(1)
l2=i/s(1)

; Plot
; Shaded surface - not terribly useful but might look pretty...
if (keyword_set(surf)) then begin

  case surf of
  1:  shade_surf,spec.spec,spec.w,spec.f,_extra=e 
  2:  surface,spec.spec,spec.w,spec.f,_extra=e
  else: message,'surf=1 (colour) or 2 (wireframe)'
  endcase
 
; More commonly, contour
endif else begin

; Basic contours
  if (not keyword_set(title)) then title='2-d spectral plot!CMax value: '+shtstr(mx)+' at ['+shtstr(spec.w(l1))+','+shtstr(spec.f(l2))+']'
  xtitle='Wavenumber (-ve = E''wards)'
  ytitle='Frequency'
  plot,/nodata,spec.w,spec.f,xs=5,ys=5,ymarg=ymarg,xmarg=xmarg,xtitle=xtitle,ytitle=ytitle,title=title,yr=yr,xr=xr
; Add second y-axis with periods
  for i=0,n_elements(spec.f)-1 do $
    if (spec.f(i) le yr(1) and spec.f(i) ge yr(0)) then begin
      if (spec.f(i) eq 0) then txt='!M$!N' $		; Infinity
      else txt=shtstr(1/spec.f(i),dp=1)
      plots,[xr(1),xr(1)+(xr(1)-xr(0))*1/40.],spec.f(i)*[1,1]
      xyouts,xr(1),spec.f(i),' '+txt,align=0
    endif

; Contour. 
  contour,spec.spec,spec.w,spec.f,c_lev=lev,c_col=cols

; Plot "*" for max value
  plots,spec.w(l1),spec.f(l2),psym=2

; Re-add axes
  axis,xax=2,/xs,xr=xr,xtitle='Wavenumber (-ve = E''wards)',xcharsize=!p.charsize
  if (keyword_set(noyax)) then $
    plots,xr(0)*[1,1],yr $
  else $
    axis,yax=2,/ys,yr=yr,ytitle='Frequency',ycharsize=!p.charsize

endelse

end