function norh_beam,index,marx=marx,xbeam=xbeam
;+
; NAME:
;       NORH_BEAM
;
; PURPOSE:
;       Create a projected image of the clean beam
;
; CALLING SEQUENCE:
;       beam=norh_beam(index)
;
; INPUT:
;       index: index structure
;
; OUTPUT:
;       returns an float array.
;       The default array size is 21^2 and
;       the center of the beam is (10,10).
;
; OPTIONAL INPUT:
;       marx: array size (default is 21)
;
; MODIFICATION HISTORY:
;       Y. Hanaoka, July, 1993.
;       T. Yokoyama, Sep, 1998.
;       Modified 2000-04-19 TY based on V. Grechnev's suggestion
;-


if n_elements(marx) le 0 then marx=21

solp=get_rb0p(index,/pangle,/quiet)
pmat=norh_gt_pmat(index)

sec_per_pix=norh_gt_pixsz(index)
sec_per_pix_dty=norh_gt_pixsd(index)
ratio_c2d=sec_per_pix/sec_per_pix_dty
rbeam=norh_gt_efl(index)

ns2ew=norh_gt_ns2ew(index)
delta=1./rbeam^2


mx=n_elements(index)

beam=fltarr(marx, marx, mx)

ii=(-marx/2+lindgen(marx))#replicate(1, marx)
jj=replicate(1, marx)#transpose(-marx/2+lindgen(marx))

for m=0, mx-1 do begin

  x=(cos(solp(m))*ii -sin(solp(m))*jj)*ratio_c2d(m)
  y=(sin(solp(m))*ii +cos(solp(m))*jj)*ratio_c2d(m)

  xx= pmat(0,m)*x+pmat(2,m)*y
  yy=(pmat(1,m)*x+pmat(3,m)*y)*ns2ew(m)

  beam(*,*,m)=exp(-delta(m)*(xx^2+yy^2))

endfor

xbeam=sec_per_pix*findgen(marx)

return,beam

end