function smooc,ar,nh,nv,xsize,ysize,minx,maxx,miny,maxy,difmax,i2
;
; This function smooths an image file within the boundaries of the windowed
; area so that point glitches, cosmic rays, especially are removed from
; the images.
;
	i2 = ar
;
    iavx = xsize
	iavy = ysize
;
    ivv  = iavy-1
	if miny gt ivv then ivv = miny
;
	ivp  = nv - iavy - 1
	if maxy lt ivp then ivp = maxy
;
    ihh  = iavx-1
	if minx gt ihh then ihh = minx
;
	ihp  = nh - iavx - 1
	if maxx lt ihp then ihp = maxx
;
	print, ' Enter Smoothing'
;
  nt = 1
  over:
	  nx1 = 0L
        FOR iyy=ivv,ivp do BEGIN
          FOR ixx=ihh,ihp do BEGIN
		    arsur = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $
			+ ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1) $
			+ ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0
		    ardiff = abs(ar(ixx,iyy) - arsur)
		    if ardiff gt difmax then begin
			  i2(ixx,iyy) = arsur
			  nx1 = nx1 + 1L
		    endif
       	  ENDFOR
        ENDFOR
	  if nx1 gt 0 then print, 'Pass', nt, '  found', nx1, ' points to smooth'
	  ar = i2
	  nt = nt + 1
	if nx1 lt nh*nv*0.001 then goto, endcos
  goto, over
 endcos:
;
end
;
  function cosrmv,ar,a2,a3,nh,nv,difmax,i2
;
; This function takes 3 image files within the boundaries of the windowed
; area and removes point glitches, cosmic rays, especially are removed from
; the images.
;
	i2 = ar
;
	print, ' Enter Cosmic Ray Removal'
;
  nt = 1
  over:
	  nx1 = 0L
        FOR iyy=1,nv-2 do BEGIN
          FOR ixx=1,nh-2 do BEGIN
		    ar23ave = (a2(ixx,iyy) + a3(ixx,iyy))/2.0

		    ardiff = abs(ar(ixx,iyy)-ar23ave)
		    if ardiff gt difmax then begin
			  i2(ixx,iyy) = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $
			    + ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1)        $
			    + ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0
			  nx1 = nx1 + 1L
		    endif
       	  ENDFOR
        ENDFOR
	  if nx1 gt 0 then print, 'Pass', nt, '  found', nx1, ' points to remove'
	  ar = i2
	  nt = nt + 1
	if nx1 lt nh*nv*0.001 then goto, endcos
	if nt gt 2 then goto, endcos
  goto, over
 endcos:
;
end
;
pro C1_WSMDIF,imn=imn,imset=imset,img=img
;
  if n_elements(imn) eq 0 then imn = '1201'
  if n_elements(imset) eq 0 then imset = 'WSMC1_'
  if n_elements(img) eq 0 then img = 5955
;
namex = 'D:\Lasco Data\C1 Data\'
name='D:\LASCO Data\C1 Differences\'
  nami  = '5954'
  nami2 = '5955'
  nami3 = '5956'
  fini1 = '.fts'
  fini2 = '.dif'
  ifstim = 5955
  ilstim = 5963
;
  nameb = namex+imn+nami+fini1
;
  base = lasco_readfits(nameb,hdrb)
  bse = readfits(nameb,hdrbe)
  exptb = fxpar(hdrbe,'EXPTIME')
  base = float(base)/exptb
;
  nameb2 = namex+imn+nami2+fini1
  a2 = lasco_readfits(nameb2,hdrb3)
  bse = readfits(nameb2,hdrbe2)
  exptb2 = fxpar(hdrbe2,'EXPTIME')
  a2 = float(a2)/exptb2
;
  nameb3 = namex+imn+nami3+fini1
  a3 = lasco_readfits(nameb3,hdrb3)
  bse = readfits(nameb3,hdrbe3)
  exptb3 = fxpar(hdrbe3,'EXPTIME')
  a3 = float(a3)/exptb3
;
  hs  = hdrb.NAXIS1
  vs  = hdrb.NAXIS2
  sun = getsuncen(hdrb)
  hc  = sun.xcen
  vc  = sun.ycen
  print, 'image set base ', hs,vs,hc,vc,exptb
;
  basen0 = base
  name0  = name+imset+'0'+fini2
  basen0 = fix(basen0*exptb)
  writefits, name0,basen0,hdrbe
;
; Plot out base image
;
  tv, basen0
;
; Remove cosmic rays from the base image
;
  print, 'Total number of points =', hs*vs
  difmax1 = 3.
  cosrm = cosrmv(base,a2,a3,hs,vs,difmax1,base)
;
; Plotthe base image
;
  basen0 = base
  name0  = name+imset+'0'+fini2
  basen0 = fix(basen0*exptb)
;
; Plot out base image with cosmic rays removed
;
  tv, basen0
;
  difmax2 = 30.
  smo   =  smooc(base,hs,vs,1,1,1,hs,1,vs,difmax2,base)
;
; Plot and put the base image into a file
;
  basen0 = base
  name0  = name+imset+'0'+fini2
  basen0 = fix(basen0*exptb)
  writefits, name0,basen0,hdrbe
;
; Plot out base image with cosmic rays removed
;
  tv, basen0
;
 istop = 0
;  istop = 1
  if istop eq 1 then goto, end2
;
; Smoothed base
;
  base_sm = smooth(base,15)
  bases    = base_sm
  basessmo = base_sm
;
; Begin determination of flat field analysis
;
  bases    = 0.
  basessmo = 0.
  img1 = img+1
start1:
;
  namei = namex+imn+strcompress(img1,/rem)+fini1
  second = readfits(namei,hdr)
  expts  = fxpar(hdr,'EXPTIME')
  second = float(second)/expts
;
  namei2 = namex+imn+strcompress(img1-1,/rem)+fini1
  a2     = readfits(namei2,hdr2)
  expts2 = fxpar(hdr2,'EXPTIME')
  a2     = float(a2)/expts2
;
  if (img1+1) gt ilstim then begin
    namei3 = namex+imn+strcompress(img1-2,/rem)+fini1
  endif else begin
    namei3 = namex+imn+strcompress(img1+1,/rem)+fini1
  endelse
  a3     = readfits(namei3,hdr3)
  expts3 = fxpar(hdr3,'EXPTIME')
  a3     = float(a3)/expts3
;
  cosrm = cosrmv(second,a2,a3,hs,vs,difmax1,second)
  smo   =  smooc(second,hs,vs,1,1,1,hs,1,vs,difmax2,second)
;
;  second = second - base_sm
  basesmo = smooth(second,15)
  second  = second  - basesmo

  bases    = bases    + second
  basessmo = basessmo + basesmo
;
 img1 = img1 + 1
 if img1 gt ilstim then goto, end1
 goto, start1
;
 end1:
;
; Flat field image multiplier
;
 baseimgmult = 1. - bases/basessmo
;
 FOR iy=45,55 DO BEGIN
  FOR ix=45,55 DO BEGIN
   print, ' Sample image multiplier at',ix,iy, bases(ix,iy)/(img1-img-1), $
          base_sm(ix,iy), (img1-img-1), baseimgmult(ix,iy)
  ENDFOR
 ENDFOR
;
; Now remove smoothed base from individual flat-fielded images
;
 img2 = img
start2:

  namen = name+imset+strcompress(img2,/rem)+fini2
;
  namei = namex+imn+strcompress(img2,/rem)+fini1
  second = readfits(namei,hdr)
  expts  = fxpar(hdr,'EXPTIME')
  second = float(second)/expts
;
  if (img2-1) ge ifstim then begin
    namei2 = namex+imn+strcompress(img2-1,/rem)+fini1
  endif else begin
    namei2 = namex+imn+strcompress(img2+2,/rem)+fini1
  endelse
  a2     = readfits(namei2,hdr2)
  expts2 = fxpar(hdr2,'EXPTIME')
  a2     = float(a2)/expts2
;
  if (img2+1) gt ilstim then begin
    namei3 = namex+imn+strcompress(img2-2,/rem)+fini1
  endif else begin
    namei3 = namex+imn+strcompress(img2+1,/rem)+fini1
  endelse
  a3     = readfits(namei3,hdr3)
  expts3 = fxpar(hdr3,'EXPTIME')
  a3     = float(a3)/expts3
;
  cosrm = cosrmv(second,a2,a3,hs,vs,difmax1,second)
  smo   =  smooc(second,hs,vs,1,1,1,hs,1,vs,difmax2,second)
;
  second = second*baseimgmult - base_sm
;
  second = fix(second*expts)
  writefits, namen, second, hdr
;
  help, second
  tv, second
 img2 = img2 + 1
 if img2 gt ilstim then goto, end2
 goto, start2
 end2:
   print, 'End of analysis'
;
end