pro embed9, serial, qm2, array, pos, out, out2, mposout, ccdposout, sub10=sub10, missing=missing, correct=correct, save=save
;
;HISTORY:
;	Written Jan/Feb by M.Morrison
;	 8-Feb-94 (MDM) - Modified to use FITS files instead of IDL save files
;	18-Feb-94 (MDM) - Changed calling sequence to use qm2 instead of
;			  m1=m1 and m2=m2
;	23-Feb-94 (MDM) - Added POSOUT output variable
;			- Added round off for insertion position (it was
;			  truncating before)
;			- Added /CORRECT option to shift the image to
;			  correct for variations due to angle
;	24-Feb-94 (MDM) - added CCDPOSOUT
;	25-Feb-94 (MDM) - Added /SAVE option
;
item_list = ['_m1', '_m2']
item = item_list( keyword_set(qm2) )
filnam = map_filnam(serial, 0, sub10, item=item)
array = rfits(filnam)
nx = n_elements(array(*,0))
ny = n_elements(array(0,*))
;
if (n_elements(missing) eq 0) then missing = 0
array = fltarr(nx, ny, 9)
pos = fltarr(2, 9)
mposout = fltarr(2,9)
ccdposout = fltarr(2,9)
for ipos=0,8 do begin
    filnam = map_filnam(serial, ipos, sub10, item=item)
    print, 'EMBED9: Reading file: ', filnam
    img = rfits(filnam, head=head)
    array(0,0,ipos) = img
    pos(0,ipos) = sxpar(head, 'LASPOS_X')
    pos(1,ipos) = sxpar(head, 'LASPOS_Y')
end
;
xmat = (lindgen(nx, ny) mod nx) - nx/2
ymat = (lindgen(nx, ny) / nx ) - ny/2
rmat = sqrt(xmat^2 + ymat^2)
ss = where(rmat lt nx/2-1)
;
;rad1 = 10*nx/64 	& rad2 = 20*nx/64
;rad2 = 24.06*nx/64 	& rad1 = rad2/2.
rad2 = 20*nx/64 	& rad1 = rad2/2.
if (max(pos(0,*))-min(pos(0,*)) le 25) then begin
    print, 'Small radius circle'
    rad=rad1
end else begin
    print, 'Large radius circle'
    rad=rad2
end
rad0 = (rad * 1024/64.)/512.		;position on CCD in relative units (-1 to 1)
;
if (qm2) then begin
    rad = rad*121./100 
    corr_slope = -0.04415
end else begin
    rad = rad*69.9/100
    corr_slope = 0.06197
end
;
;;if (keyword_set(m1)) then rad = rad*69.9/100
;;if (keyword_set(m2)) then rad = rad*121./100
;
if (pos(0,1) gt pos(0,0)) then begin
    phase = !pi
end else begin
    print, 'Starting position not "normal" - image needs to be rotated?'
    phase = 0
end
;
nout = nx*2		;should be 128 for nominal case
nxy = nx*ny
out = fltarr(nout, nout, 9) + missing
out2 = fltarr(nout, nout)
overmat = intarr(nout, nout)
;
;;x00 = nout/2 - nx/2
;;y00 = nout/2 - ny/2
;;out(x00, x00, 0) = array(*,*,0)
x00 = nout/2
y00 = nout/2
ss2 = xmat(ss)+x00 + (ymat(ss)+y00)*nout
out(ss2) = array(ss)
out2(ss2) = array(ss)
overmat(ss2) = overmat(ss2) + 1
;
for i=0,7 do begin
    x0 = cos(-2*!pi*i/8.+phase)*rad
    y0 = sin(-2*!pi*i/8.+phase)*rad
    mposout(0,i+1) = x0
    mposout(1,i+1) = y0
    ccdposout(0,i+1) = cos(-2*!pi*i/8.+phase)*rad0
    ccdposout(1,i+1) = sin(-2*!pi*i/8.+phase)*rad0
    ;
    corr_val = 0
    if (keyword_set(correct)) then begin
	corr_val = y0*(1024/ny)*corr_slope	;1024/ny is normally 16 since these are rebinned to 64x64
	print, 'Angle offset correction: ', corr_val, ' for ', y0
    end
    x0 = fix(x0 + sgn(x0)*0.5)	;round off
    y0 = fix(y0 + sgn(y0)*0.5)
    print, 'EMBED9: Position shift', i+1, x0, y0
    ;;out(x0+x00, y0+y00, i+1) = array(*,*,i+1)
    ;;out(xmat(ss)+x0+x00,ymat(ss)+y0+y00,i+1) = array(ss+nx*ny*(i+1))
    ss2 = xmat(ss)+x00+x0 + (ymat(ss)+y00+y0)*nout
    out(ss2 + nout*nout*(i+1)) = array(ss + nx*ny*(i+1))
    out2(ss2) = out2(ss2) + array(ss + nx*ny*(i+1)) - corr_val
    overmat(ss2) = overmat(ss2) + 1
end
;
out2 = out2 / (overmat>1)
ss3 = where(overmat eq 0)
if (ss3(0) ne -1) then out2(ss3) = missing
;
if (keyword_set(save)) then begin
    item_list = ['_m1embed', '_m2embed']
    item = item_list( keyword_set(qm2) )
    filnam = map_filnam(serial, 0, sub10, item=item)
    wrt_fits, filnam, dummy_head, out2
end

end