;+
; PURPOSE:
;  This procedure stitches together individual binned skycells from
;  gpc1. The program assumes that the skycells form a 12x12 grid on
;  the sky, and does not make any attempt to correct for spherical
;  geometry projection effects. The 12x12 tesselation grid is created
;  by IPP commands like
;
;   skycells -mode LOCAL -scale 0.2 -nx 12 -ny 12 -size 4 4 -fix-ns
;     -center "racen" "deccen" -D CATDIR "catdir name"
;
; CATEGORY:
;  PS PR images
;
; CALLING SEQUENCE:
;  stitch, filter
;
; INPUTS:
;  filter: A single character denoting the filter to stitch 
;  (see usage notes below).
; 
; NOTES ON USAGE:
;  This script is highly specific to the way that I create
;  mosaics:
;    -First, I create 8x8 binned images of each skycell in mana,
;     naming the binned files X.YYY.b8.fits. Here, X is one of g,r, or i,
;     and indicates the filter name. 
;    -Next, I copy the files from each filter into a directories named
;     g, r, and i.
;    -I run stitch from the directory above each of these. The program
;     goes into the appropriate subdirectory, finds the files named
;     X.YYY.b8.fits, and stitches them together. It outputs a file
;     named X.merge.fits in the current directory.
;-
pro stitch, filter

  files = file_search(filter+'/*.fits', count = count)
  if count eq 0 then begin
     print, 'no fits files files found'
     return
  endif

  ;- get the image dimensions of a sample file
  im = mrdfits(files[0],0,h,/silent)
  nx = n_elements(im[*,0])
  ny = n_elements(im[*,1])
  grid = 12
  
  ;- create the mosaic array
  big = fltarr(1d * nx * grid, 1D * ny * grid)
  xind = rebin( indgen(nx), nx, ny)
  yind = rebin( 1#indgen(ny), nx, ny)

  xs = indgen(grid^2) mod grid
  ys = indgen(grid^2) / grid
  
  ;- a plot to show which files are present and which are missing
  plot, xs, ys, xra = [-1, 13], yra = [-1, 13], psym = 4, /xsty, /ysty
  missing = xs * 0

  red = '0000ff'xl
  white = 'ffffff'xl
  ;- read in and stitch each image
  for i = 1, grid^2-1, 1 do begin

     f=filter+'/'+filter+'.'+string(i,format='(i3.3)')+'.b8.fits'
     if ~file_test(f) then color = red else color = white
     xyouts, xs[i], ys[i], strtrim(i, 2), color = color, $
             charsize = 2
     
     if ~file_test(f) then continue
     
     xoff = ((i) mod grid) * nx
     yoff = ((i) / grid) * ny
     im = mrdfits(f,0,h,/silent)
    ;- the IPP doesn't seem to populate the outer 1 pixel border of
    ;- each skycell. Interpolate over these
     sz = size(im)
     im[0,*] = im[1,*]
     im[sz[1]-1, *] = im[sz[1]-2,*]
     im[*,0] = im[*,1]
     im[*, sz[2]-1] = im[*,sz[2]-2]
     big[xind + xoff, yind + yoff] = im
  endfor

  ;- save the mosaic
  print, 'writing fits'
  writefits, filter+'.merge.fits', big

  return
end