@gaussfit3
pro plot_lyot, code, ser_arr, lyot_map, trans_map, unc_ok=unc_ok, dev_ok=dev_ok, hc=hc, qstop=qstop, outfil=outfil
;
;
;Ver 1.0
;Ver 2.0  8-Feb-94 (MDM) - Modified to use FITS files
;			 - Cleaned up and label nicely
;Ver 2.1  8-Feb-94 (MDM) - Use GAUSSFIT instead of center
;			 - Normalized for Laser power changes
;			 - Plot line traces too
;Ver 2.2  9-Feb-94 (MDM) - Corrected formating errors
;			 - Sorted by wavelength (fix in GET_SCAN_SER)
;			 - Added /QSTOP
;			 - On the multi plot, plot every 8 instead of every one
;			 - Added OUTFIL option for hardcopy (to be able to run twice)
;Ver 2.3 10-Feb-94 (MDM) - Stopped inverting the color table (now black is low signal)
;			 - List the series in increasing order (not by wavelength)
;
progver = 'PLOT_LYOT.PRO Ver 2.3'
;
f = 0.75	;NCD
;
if (n_elements(outfil) eq 0) then outfil = 'idl.ps'
;
if (keyword_set(hc)) then begin
    set_plot, 'ps'
    device, file=outfil
end
clearplot
erase
;
ser_arr = get_scan_ser(code, ver=tit1, tit=tit, dev_ok=dev_ok, unc_ok=unc_ok)
laser_info, ser_arr, wave, lpow
;
n = n_elements(ser_arr)
trans_map = fltarr(64,64,n)
for i=0,n-1 do begin
    infil = map_filnam(ser_arr(i), 0)
    print, 'Reading: ', infil
    img = rfits(infil)
    trans_map(0,0,i) = img
end
;
tit = 'LYOT Map using Series: ' + arr2str(ser_arr(sort(ser_arr)))
;
nx = 64L
ny = nx
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)
;
plot, wave, trans_map(0,0,*), /nodata, yrange=[0,max(trans_map)/min(lpow)], tit=tit
xarr = indgen(n)*nx*nx
;for i=0,63 do for j=0,63 do oplot, wave, trans_map(i,j,*)
for i=0,n_elements(ss)-1,8 do oplot, wave, trans_map(xarr+ss(i))/lpow
plottime, 0, 0, progver & plottime & pause
;
if (keyword_set(qstop)) then stop
;
;; !p.multi=[0,4,4] & for i=0,n_elements(ss)-1 do plot, wave, trans_map(xarr+ss(i))
;
wave = double(wave)
lyot_map = dblarr(64,64)
for i=0,63 do begin
    print, 'i = ', i
    for j=0,63 do begin
	;lyot_map(i,j) = total(wave*trans_map(i,j,*))/total(trans_map(i,j,*))
	yfit=gaussfit3(wave, reform(trans_map(i,j,*))/lpow,coeff)
	lyot_map(i,j) = coeff(1)
    end
end
ss = where(rmat lt nx/2-1)
avg = total(lyot_map(ss)-6767)/n_elements(ss) + 6767

ss = where(rmat ge nx/2-1)
lyot_map(ss) = avg
;
!p.multi = [0,1,2]
plot, lyot_map(*,32)-6767, tit=tit, xtit='X Pixels (row #32 plot)', ytit='Lyot Wavelength (A) (from 6767)', /ynoz
plot, lyot_map(32,*)-6767, tit=tit, xtit='Y Pixels (col #32 plot)', ytit='Lyot Wavelength (A) (from 6767)', /ynoz
plottime, 0, 0, progver & plottime & pause
!p.multi = 0
if (keyword_set(hc)) then pprint
;
if (!d.name eq 'PS') or (getenv('ys_batch') ne "") then begin
    set_plot, 'Z'
end
nxout = 640
nyout = 512
case !d.name of
   'X': wdef, 1, nxout, nyout
   'Z': device, set_res=[nxout, nyout]
end
;
imin = 6767.86D
imax = 6767.93D
bimg = bytscl(lyot_map, min=imin, max=imax)
;
erase
!color = 0
tv, bytarr(nxout, nyout)+255
tvscl, rebin(bimg, 64*7, 64*7, /sample)
mk_refbar, 64*7.5, 0, imin, imax, /delta, size=f
xyouts, 0, 64*7.5, tit, /device, siz=1.5*f
xyouts, 0, 64*7, progver + '   Program Run: ' + !stime, /device, siz=f
xyouts, 64*7.5, 300, 'Lyot Average: ' + strtrim(avg,2), /dev, siz=f
!color = 255

if (keyword_set(hc)) then begin
    img = tvrd()
    olddev = !d.name
    set_plot, 'ps'
    XPAGE=7.0 & YPAGE=9.5 & XOFFS=0.75 & YOFFS=1  ;INCHES
    DEVICE, /INCHES,XSIZE=XPAGE,YSIZE=YPAGE,XOFFSET=XOFFS, YOFFSET=YOFFS, file=outfil
    tv, img
    pprint
    set_plot, olddev
end

end