pro disp_mdidust, seq_filnam, ilaser, iman, ihr, images, dark_images, $
		posx, posy, dfposx, dfposy, rx, rfac, rms, cpox, hc=hc, one_page=one_page
;
;
;
progver = 'MDIDUST  Ver 1.1 '
;
f = 1.0
if (keyword_set(hc)) then begin
    save_device = !d.name
    color_save = !color
    set_plot, 'ps
    !color = 0
end

nx = 200	;n_elements(cpox(*,0,0))
ny = 200	;n_elements(cpox(0,*,0))
n  = n_elements(cpox(0,0,*))
np = n_elements(posx)
;
nnout = nx*1.2
nxout = nnout*3
nyout = nnout*4
;
case !d.name of
   'X': wdef, 1, nxout, nyout, /already
   'PS': tv2, nxout, nyout, /init
end

man_arr = ['Default Fixed Values', 'Manually Defined', 'Center of Gravity', 'SAT_OBS_LOC']
laser_arr = ['OBS Mode', 'CAL Mode']
hr_arr = ['Full Disk (FD) Pixels', 'High Res (HR) Pixels']
;
step = 15
npage = n / 3
if ((n mod 3) ne 0) then npage = npage + 1
;
npp = 3		;number of images per page
for iimg=0,n-1 do begin
    if ((iimg mod npp) eq 0) then begin
	erase
	ipage = iimg/npp + 1

	xyouts2, 0, nyout-step*1, /dev, seq_filnam, siz=1.5
	xyouts2, 0, nyout-step*3, /dev, progver + ' Program Run: ' + !stime
	xyouts2, nyout, nyout-step*3, /dev, 'Page ' + strtrim(ipage,2) + ' of ' + strtrim(npage,2), align=1.0

	xyouts2, 0, nyout-step*6, /dev, 'Image #''s used: ' + string(images, format='(99i4)')
	xyouts2, 0, nyout-step*7, /dev, 'Dark #''s used:  ' + string(dark_images, format='(99i4)')

	xyouts2, 0, nyout-step*9, /dev, 'OBS X Position: '
	for i=0,np-1 do xyouts2, 150+i*55, nyout-step*9, /dev, string(posx(i), format='(f6.3)'), align=1.0, siz=f*0.8
	for i=0,np-1 do xyouts2, 150+i*55, nyout-step*10, /dev, string(dfposx(i), format='("(",f6.3,")")'), align=1.0, siz=f*0.8
	xyouts2, 0, nyout-step*11, /dev, 'OBS Y Position: '
	for i=0,np-1 do xyouts2, 150+i*55, nyout-step*11, /dev, string(posy(i), format='(f6.3)'), align=1.0, siz=f*0.8
	for i=0,np-1 do xyouts2, 150+i*55, nyout-step*12, /dev, string(dfposy(i), format='("(",f6.3,")")'), align=1.0, siz=f*0.8
        xyouts2, 0, nyout-step*13, /dev, 'Radius Factor: ' + strtrim(rfac, 2)

	xyouts2, 0, nyout-step*15, /dev, 'Image Mode:       ' + laser_arr(ilaser)
	xyouts2, 0, nyout-step*16, /dev, 'Image Resolution:  ' + hr_arr(ihr)
	xyouts2, 0, nyout-step*17, /dev, 'Locations Derived: ' + man_arr(iman)

	xyouts2, 0*nnout, 3*nnout-.2*nx+2, /dev, 'Image Scale +/- 3 sigma', siz=f*.8
	xyouts2, 2*nnout, 3*nnout-.2*nx+2, /dev, 'Image Scale +/- 1 sigma', siz=f*.8
    end
    ;
    ii = iimg mod npp
    y0 = ((npp-ii)-1) * nnout
    x1 = 0
    x2 = 2*nnout
    cpox0 = cpox(*,*,iimg)
    idev = stdev(cpox0, iavg)
    bimg1 = congrid(bytscl(cpox0, -3*rms, 3*rms), nx, ny)
    bimg2 = congrid(bytscl(cpox0, -1*rms, 1*rms), nx, ny)
    tv2, bimg1, x1, y0
    tv2, bimg2, x2, y0
    xyouts2, nx*1.05, y0+ny-2*step, /dev, 'Dust Surface Radius:   ' + strtrim(rx(iimg), 2)
    xyouts2, nx*1.05, y0+ny-3*step, /dev, 'Dust Surface:  ' + rad2location(rx(iimg), ihr)
    xyouts2, nx*1.05, y0+ny-5*step, /dev, 'STDDEV:   ' + strtrim(idev,2)

    if ( ((iimg mod npp) eq npp-1) or (iimg eq n-1)) then begin
	pause
	if (keyword_set(one_page)) then return
	if (!d.name eq 'PS') then pprint
    end
end


if (!d.name eq 'PS') then begin
    !color = color_save
    set_plot, save_device
end
;
if (keyword_set(qstop)) then stop
end
;------------------------------------------------------------------------------
pro mdidust, seq_filnam, pox, cpox, obsmode=obsmode, fd=fd, $
	hole_manual=hole_manual, hole_centroid=hole_centroid, $
	rad_arr=rad_arr, hc=hc, one_page=one_page, $
	corr_image=corr_image, force_nseries=force_nseries
;+
;NAME:
;	mdidust
;PURPOSE:
;	Program to make dust images at various positions through
;	the instrument.
;SAMPLE CALLING SEQEUENCE:
;	mdidust, '/md66/log/sci160k/940402_033500'
;	mdidust, '950630_214927', corr_image=[0,0,0,0,0,0,-1,-1,-1], force_nseries=245
;INPUT:
;	seq_filnam
;OPTIONAL KEYWORD INPUT:
;	obsmode	- Set if the observations are in obs mode (default is
;		  cal mode)
;	fd	- Set if you wish to analyze the Full Disk images
;		  (default is to analyze the high res (HR) images)
;	hole_man- If set, allow the user to manually select the hole
;		  positions
;	hole_centroid - If set, derive the hole positions by finding
;		  the center of gravity of the images
;	rad_arr	- The radius array of the positions through the 
;		  instrument which should be checked in "funny" units
;		  (see Jesper Schou)
;	corr_image - An array of the correction offsets to get to the
;		  correct image (for cases when an image is missing)
;OUTPUT:
;	pox	- input images used
;	cpox	- dust image for each radius location 
;METHOD:
;	Corrects for number of images used at each point
;	reads raw images and removes gradients and limb-darkening
;HISTORY:
;	Originally Jesper Schou's DUST2.PRO
;	Written 4-Apr-94 by M.Morrison
;V1.0	 4-Apr-94 (MDM) - Made it a procedure, not a program
;			- Removed all of the questions and
;			  added keywords for the options
;			- Changed format of output quite a bit
;V1.1	 5-Apr-94 (MDM) - Changed format/info printed
;V1.2	25-Apr-94 (MDM) - Added /ONE_PAGE option
;V1.3	26-Jul-95 (MDM) - Added CORR_IMAGE optional input and
;			  FORCE_NSERIES
;-
;
;---------------------------------------- Initialization
;
ilaser = 1	;default is calmode
if (keyword_set(obsmode)) then ilaser = 0
;
iman = 3	;use SAT_OBS_LOC
if (keyword_set(hole_fixed)) then iman = 0	;set to fixed location
if (keyword_set(hole_manual)) then iman = 1	;set hole positions manually
if (keyword_set(hole_centroid)) then iman = 2	;to use center of gravity of hole images
;
ihr = 1
if (keyword_set(fd)) then ihr = 0
;
;---------------------------------------- Check the data
;
read_mdi, seq_filnam, 0, nseries, /nimages
if (keyword_set(force_nseries)) then nseries = force_nseries
;
case nseries of
    65: begin
	  nskip=7
	; set offset in series to find appropriate images
	  icaldark=1
	  caldarkoffset=0
	  caloffset=1
	  darkoffset=7
	  if (ihr eq 0) then begin
	    imoffset=2
	  end else begin
	    imoffset=6
	  end
	  nim=(nseries-2)/nskip
	end
    245: begin
	  nskip=27
	; set offset in series to find appropriate images
	  icaldark=1
	  caldarkoffset=0
	  caloffset=1
	  darkoffset=7
	  if (ihr eq 0) then begin
	    imoffset=2
	  end else begin
	    imoffset=6
	  end
	  nim=9
	end
    else: begin
	  print,'Incorrect number of images:',nseries
	  print,'Type .con if you want to try to see what you can get'
	  stop

	  nskip=7
	; set offset in series to find appropriate images
	  icaldark=1
	  caldarkoffset=0
	  caloffset=1
	  darkoffset=7
	  if (ihr eq 0) then begin
	    imoffset=2
	  end else begin
	    imoffset=6
	  end
	  nim=(nseries-2)/nskip
	end
endcase
;
;---------------------------------------- Read the data
;
images = indgen(nim)*nskip + imoffset
if (n_elements(corr_image) ne 0) then images = images - corr_image	;MDM 26-Jul-95
read_mdi, seq_filnam, images, pox
;
dark_images = indgen(nim)*nskip + darkoffset
if (n_elements(corr_image) ne 0) then dark_images = dark_images - corr_image	;MDM 26-Jul-95
read_mdi, seq_filnam, dark_images, dark
;
pox = pox - dark
pox = float(pox)
;
q = size(pox)
size1=q(1)
size11=size1-1
size12=size11/2.
; set variables depending on mode and output format
; Bill's notation for the target radii (rx) is used
if (ilaser eq 0) then begin
  rfac=-1.0
  rcal=50./61.*136.
  if (ihr eq 0) then begin
    rx=[-34,-10,0,15,23,40]
  end else begin
    rx=[-36,-35,0,6,7,8,9,10,44]
  end
end else begin ; laser
  rfac=128./136.*61./50.
  rcal=128.
  if (ihr eq 0) then begin
    rx=[-97,-60,-49,-39,-36,-25,-10,0,15,97]
;   rx=[-11,-10,0,3,4,15]
  end else begin
    rx=[-126,-91,-35,0,7,9,46]
  end
end
if (keyword_set(rad_arr)) then rx = rad_arr
;
q=size(rx)
nrx=q(1)
if (ihr eq 0) then begin
  r0=114.*size1/256
end else begin
  r0=300.
end
;
;set up various variables
x0=size12
y0=size12
dx=rebin(reform(findgen(size1)-x0,size1,1),size1,size1)
dy=rebin(reform(findgen(size1)-y0,1,size1),size1,size1)
r=sqrt(dx^2+dy^2)
;
;---------------------------------------- set hole positions
posx=fltarr(nim)
posy=fltarr(nim)
caldark=0
if (icaldark eq 0) then read_mdi, seq_filnam, caldarkoffset, caldark
;
fposx = posx
fposy = posy
angle=findgen(8)*3.1416/4.
fposx(0:7)=0.62*cos(-angle)		;fixed positions
fposy(0:7)=0.62*sin(-angle)
fposx = [0.0, fposx]
fposy = [0.0, fposy]
fposx=fposx*size1/256
fposy=fposy*size1/256
case iman of
    0: begin		;use fixed positions
	  posx = fposx
	  posy = fposy
	end
    1: begin		;manually select (click)
	    mag=2
	    print,'Please click on the position of the hole on each image'
	    for i=0,nim-1 do begin
	      read_mdi, seq_filnam, nskip*i+caloffset, cal
	      cal = cal - caldark
	      tvscl,rebin(cal,mag*256,mag*256)
	      print,i
	      cursor,xc,yc,/wait,/device
	      posx(i)=(float(xc)/mag-size12)/rcal
	      posy(i)=(float(yc)/mag-size12)/rcal
	      print,posx(i),posy(i)
	    end
	end
    2: begin		;center of gravity
	    print,'The center position of each hole/laser spot will be marked'
	    print,'by a cross on each image. If the crosses are off use manual'
	    print,'mode. If some images do not show a hole, give up!'
	    for i=0,nim-1 do begin
	      read_mdi, seq_filnam, nskip*i+caloffset, cal
	      cal = cal - caldark
	      acal=2*total(cal)/size1^2
	      cal1=0*cal+0.0
	      cal1(where(cal ge acal))=1
	      xc=total(dx*cal1)/total(cal1)+x0
	      yc=total(dy*cal1)/total(cal1)+y0
	      cal(xc-5:xc+5,yc)=0
	      cal(xc,yc-5:yc+5)=0
	      tvscl,cal,i
	;     print,i
	      posx(i)=(xc-size12)/rcal
	      posy(i)=(yc-size12)/rcal
	;     print,posx(i),posy(i)
	    end
	end
    3: begin		;use SAT_OBS_LOC
	    for i=0,nim-1 do begin
	      read_mdi, seq_filnam, nskip*i+caloffset, cal
	      cal = cal - caldark
	      sat_obs_loc, cal, x00, y00
	      posx(i) = (x00-size1/2.) / (size1/2.)
	      posy(i) = (y00-size1/2.) / (size1/2.)
	    end
	end
endcase
;
;---------------------------------------- Do the work
; set up mask
mask=fix((r0-r)>0<1)
; remove the edges
mask(0,*)=0
mask(size11,*)=0
mask(*,0)=0
mask(*,size11)=0
;set up functions to fit for
nfmax=10
nf=6
f=fltarr(size1,size1,nfmax)
f(*,*,0)=1
f(*,*,1)=dx
f(*,*,2)=dy
f(*,*,3)=dx^2
f(*,*,4)=dy^2
f(*,*,5)=r^4
f(*,*,6)=r^6
f(*,*,7)=dx*dy
;normalize the functions
for i=0,nf-1 do f(*,*,i)=f(*,*,i)/sqrt(total(f(*,*,i)^2))
for i=0,nim-1 do begin
  q=reform(pox(*,*,i))
  removegrad1,q,f,nf,mask
  pox(*,*,i)=q
end
for i=0,nim-1 do pox(*,*,i)=pox(*,*,i)*mask
rms=sqrt(total(pox^2)/nim/total(mask))
pox=pox>(-6*rms)<6*rms
cpox=fltarr(size1,size1,nrx)

; make the images at different positions
print,'target rad #, nominal radius, index in movie'
for j=0,nrx-1 do begin
  print,j,rx(j),rfac*rx(j)
  r=rfac*rx(j)
  spox=fltarr(size1,size1)
  smask=fltarr(size1,size1)
  for i=0,nim-1 do begin
    sx=fix(r*posx(i))
    sy=fix(r*posy(i))
    x1=sx>0
    x2=(size11+sx)<size11
    x3=(-sx)>0
    x4=(size11-sx)<size11
    y1=sy>0
    y2=(size11+sy)<size11
    y3=(-sy)>0
    y4=(size11-sy)<size11
    spox(x1:x2,y1:y2)=spox(x1:x2,y1:y2)+pox(x3:x4,y3:y4,i)
    smask(x1:x2,y1:y2)=smask(x1:x2,y1:y2)+mask(x3:x4,y3:y4)
  end
  cpox(*,*,j)=spox/(smask>1)
end
;
disp_mdidust, seq_filnam, ilaser, iman, ihr, images, dark_images, $
			posx, posy, posx-fposx, posy-fposy, rx, rfac, rms, cpox, hc=hc, one_page=one_page
;
end