func dataread(x, name, h)
 ;x and h are returned, value of function is status
 ;1/7/2008 - add a $instrument_type flag
 ;4/13/2000 - tiff file support, not even baseline but at least some including
 ;the rare 16 bit variety
 ;4/28/99 set $multi_image_flag for appropiate FITS files, this is a hint
 ;to the calling routine to reprocess if it knows how
 ;4/22/99 load globals $data_tai, $data_xc, etc when possible
 ;6/10/97 add ippread for Jaffey's files, type 6 here
 ;5/31/95 more general, reads first bit of file to determine type
 ;4/10/95 based on some code in main browser from Zoe, reads fz and fits files
 ;at present
 h = ''		;in case the data type doesn't have a header
 ; added, LS 17apr00 (distance from Sun)
 ;added $data_focus for Solar B test data
 zeroifundefined, $data_focus, $last_type, $instrument_type
 zero, $data_time_tai, $data_xcen, $data_ycen, $data_xscale, -
 	$data_yscale, $multi_image_flag, $data_dist, $data_focus, $instrument_type		 
 ;special check for Nikon raw files, assume they have a .NEF
 if strpos(name, '.NEF') gt 0 or strpos(name, '.RAF') gt 0 then {
   ;use for both Nikon NEF files and Fuji RAF iles
   nikonrawread, x, name
   return, 1
 }
 if defined($last_type) eq 0 then $last_type = 0	;fz is our default guess
 ;we use $last_type for the first guess, 0=> fz, 1=> fits
 ;to avoid error messages, check if file exists first, should make all the
 ;the calls below less noisy
 if filesize(name) le 0 then return, -1
 ;if the guess is right, we save some time
 ;but ... because of some problems with read_jpeg for non-jpeg files on the
 ;PC's, verify that this is jpeg first by bypassing thei section
 if $last_type ne 2 then {
 ncase $last_type
   if fzreader(x, name, h) eq 1 then return, 1
   if rdfits_mdi_trace(x, name, h) eq 1 then return, 1
   if read_jpeg(x, name) eq 1 then return, 1
   if tjpeg12read(x, name) eq 1 then return, 1
   if gifread(x, name, $gif_ct) eq 1 then return, 1
   if p5read(x, name) eq 1 then return, 1
   if ippread(x, name) eq 1 then return, 1
   if tiffread(x, name) eq 1 then return, 1
 endcase
 }
 ;our guess didn't work or wasn't supported above, so now try them all
 lun = get_lun()
 if openr(lun, name) eq 0 then { ty,'no such file'  return, -1 }
 ;got something, read a bit of it
 fileread, lun,xq,0,20,0
 close, lun
 sq = smap(xq(0))
 ;ty,'first character = ', sq
 ;hex, xq
 xq0 = xq(0)
 ;do the possibilities in a row, each does a return
 if xq0 eq 0xaa or xq0 eq 0x55 then {
  ;the fz case
  if fzreader(x, name, h) eq 1 then { $last_type = 0  return, 1 } else return, 0
 }
 if sq eq 'S' then {
  ;the fits case
  ty,'the fits case'
  if rdfits_mdi_trace(x, name,h) eq 1 then
 	 { $last_type = 1  ty,'returning'   return, 1 } else return, 0
 }
 ;6/10/97, check for ipp file
 if sq eq 'n' then {
  ;we believe that these start with a ncaa
  if ippread(x,name,h) eq 1 then {$last_type = 6  return, 1 } else return, 0
  }
 ;the rest don't have headers/comments so make h a null string
 h=''
 if sq eq 'I' or sq eq 'M' then {
  ;probably a TIFF file, at least not anything else here
  if tiffread(x, name) eq 1 then { $last_type = 7  return, 1 }
  }
 if xq0 eq byte(255) then {
  ;could be a jpeg, 2 types
  if read_jpeg(x, name) eq 1 then { $last_type = 2  return, 1 }
  if tjpeg12read(x, name) eq 1 then { $last_type = 3  return, 1 } else return, 0
  }
 if sq eq 'G' then {
  ;the gif case
  if gifread(x, name, $gif_ct) eq 1 then
 	 { $last_type = 4  return, 1 } else return, 0
  }
 if sq eq 'P' then {
  ;a P5 file?
  if p5read(x, name) eq 1 then { $last_type = 5  return, 1 } else return, 0
  }
 ;nothing we know !, return a zero status
 x = 0
 return, 0
 endfunc
 ;========================================================
func fzreader(x, name, h)
 ;just a front end for fzread that also checks the header for a La Palma
 ;style date
 if fzread(x, name, h) eq 0 then return, 0
 ;check h, if we get in trouble just bail out with a success return but
 ;no time set
 iq = strpos(h,'_im')
 if iq lt 0 then return, 1
 ;this means we probably have a La Plama file
 iq = iq+3
 hh = h(iq:*)
 day = strtol(hh)
 iq = iq + 1
 if strpos('0123456789', h(iq)) ge 0 then iq = iq + 1
   
   mm = upcase(h(iq:(iq+2))
   case
     mm eq 'JAN':im=1
     mm eq 'FEB':im=2
     mm eq 'MAR':im=3
     mm eq 'APR':im=4
     mm eq 'MAY':im=5
     mm eq 'JUN':im=6
     mm eq 'JUL':im=7
     mm eq 'AUG':im=8
     mm eq 'SEP':im=9
     mm eq 'OCT':im=10
     mm eq 'NOV':im=11
     mm eq 'DEC':im=12
     else im = 0
   endcase
 if im eq 0 then return, 1
 ;now the year
 iq = iq + 3
 year = strtol(h(iq:*))
 if year lt 100 then year = year + 1900
 if year lt 1985 then year = year + 100
 ;get doy
 dy = doy(year, im, day)
 ;now the time, look for a colon
 hh = h(iq:*)
 iq = strpos(hh,':')
 if iq lt 0 then return, 1
 ;assume a standard format (make more fancy later in C code)
 hour = strtol(hh((iq-2):iq))
 m = strtol(hh((iq+1):(iq+2))
 s = float(hh( (iq+4):*)
 $data_time_tai = tai_from_date(year, dy, hour, m, s)
 return, 1
 endfunc
 ;========================================================
subr dataread, x, name, h
 ;works as an interface to the function version, just complains if it can't
 ;read the file
 h = ''		;because sometimes passing an undefined string doesn't work
 iq = dataread(x, name, h)
 if iq eq 1 then return
 x = 0
 ty,'dataread failure, bad or non-exsistent file'
 return
 endsubr
 ;========================================================
func datawrite(type_in, x, name, h)
 ;just writes a file of the specified type and passes along status
 type = type_in
 ;all tiff cases handled by array characteristics now
 if type eq 9 or type eq 11 then type = 5
 ;also for JPEG
 if type eq 10 then type = 4
 if defined(h) eq 0 then h=''
 nd = num_dim(x)
 ;;ty,'type in datawrite', type
 ncase type
  iq = fzwrite(x, name, h)
  iq = fcwrite(x, name, h)
  iq = fcrunwrite(x, name, h)
  iq = wrfits(x, name, h)
  {
  if defined($jpeg_quality) eq 0 then $jpeg_quality=75
  iq = write_jpeg(x, name, $jpeg_quality)
  }
  iq = tiffwrite(x, name)
  {
   if nd eq 3 and dimen(x,0) eq 3 then {
   ;a rgb triplet, don't have a general solution yet so ...
   xq = rgb2gray(x)
   if isarray($gif_color) then iq = gifwrite(xq, name, $gif_color) else
  	iq = gifwrite(xq, name)
   }  else {
   if isarray($gif_color) then iq = gifwrite(x, name, $gif_color) else
  	iq = gifwrite( x, name)  }
  }
  iq = pgmwrite(x, name)
  {
  if nd eq 2 then iq = wrt_sgi(x, name) else {
  	if nd eq 3 and dimen(x,0) eq 3 then {
  	imr = x(0,*,*)   img = x(1,*,*)   imb = x(2,*,*)
	iq = wrt_sgi_rgb(imr, img, imb, name)
	imr = 0   img = 0   imb = 0
	} else return, -1
  }
  }
 endcase
 return, iq
 endfunc
 ;========================================================
func data_type(name, header, params)
 ;not very clever but returns probable data type by reading a bit of the file
 ;read a bit of the file to figure out type
 ;5/6/96 added optional header and params which we often need as well
 narg = !narg
 header =''
 params = zero(lonarr(11))-1
 lun = get_lun()
 if openr(lun, name) eq 0 then { ty,'no such file'  return, -1 }
 ;got something, read a bit of it
 fileread, lun,xq,0,20,0
 close, lun
 sq = smap(xq(0))
 xq0 = xq(0)
 ;do the possibilities in a row
 if xq0 eq 0xaa or xq0 eq 0x55 then {
   ;probably an fz file
   $last_type = 0
   if narg gt 1 then fzinspect,name, params, header	;get both anyhow
   return, $last_type }

 if sq eq 'S' then {
   ;the fits case (probably)
   $last_type = 1
   if narg gt 1 then fitsinspect, name, header, params
   return, $last_type }

 if sq eq 'P' then {
   ;a P5 file?
   if narg gt 1 then p5inspect,name, params
   $last_type = 5  return, $last_type
   }
 ;6/10/97, check for ipp file
 if sq eq 'n' then {
   if narg gt 1 then ippinspect, name, header, params 
   $last_type = 6  return, $last_type
   }
 
 ;the above all had inspectors, for the remaining at least get the file size
 params(0) = filesize(name)

 if xq0 eq byte(255) then {
   ;could be a jpeg, 2 types, but we'll just guess the 8 bit type here
   $last_type = 2  return, $last_type
   }
 if sq eq 'G' then {
   ;the gif case
   $last_type = 4  return, $last_type
   }
 if sq eq 'I' or sq eq 'M' then {
  ;probably a TIFF file, at least not anything else here
  $last_type = 7  return, $last_type }
 ;nothing we know !, return a zero status
 return, -1
 endfunc
 ;========================================================
subr p5inspect(name, params)
 lun = get_lun()
 params = zero(lonarr(11)
 params(0) = -1
 if openr(lun ,name) eq 0 then { ty,'p5read: no such file'  return, 0 }
 fileread, lun, xq, 0, 2, 0
 if smap(xq) ne 'P5' then { ty,'not a P5 file'  close, lun  return, 0 }
 rewindf, lun
 s=''
 readf, lun, s
 nx = 0		ny = 0
 readf, lun, nx, ny
 dmax = 0
 readf, lun, dmax
 params(2) = 2
 params(3) = nx
 params(4) = ny
 if dmax le 255 then { type = 0 } else if dmax le 65535 then { type = 1 } else
 	type = 2
 if dmax lt 0 then type = 2	;?
 params(1) = type
 close, lun
 params(0) = filesize(name)
 endsubr
 ;========================================================
func p5read(x, name)
 lun = get_lun()
 if openr(lun ,name) eq 0 then { ty,'p5read: no such file'  return, 0 }
 fileread, lun, xq, 0, 2, 0
 if smap(xq) ne 'P5' then { ty,'not a P5 file'  close, lun  return, 0 }
 rewindf, lun
 s=''
 readf, lun, s
 nx = 0		ny = 0
 readf, lun, nx, ny
 dmax = 0
 readf, lun, dmax
 ;ty,s
 ;ty,nx,ny
 ;ty,dmax
 if dmax le 255 then { type = 0 } else if dmax le 65535 then { type = 1 } else
 	type = 2
 if dmax lt 0 then type = 2	;?
 fileread, lun, x, -1, nx*ny, type
 close, lun
 redim, x, nx, ny
 if bigendian() eq 1 then {
   ;7/26/2004 - use the builtin endian routine rather than swapb, swapl
   endian, x
   ;;if symdtype(x) eq 1 then swapb, x
   ;;if symdtype(x) eq 2 then swapl, x
 }
 return, 1
 endfunc
 ;========================================================
func pgmwrite(x, name)
 ;writes an I*1 pgm binary file, returns 1 if OK
 if symdtype(x) ne 0 or num_dim(x) ne 2 then {
 	ty,'PGMWRITE: input must be a 2-D I*1 array'  return, 0 }
 lun = get_lun()
 if openw(lun ,name) eq 0 then {
 	ty,'PGMWRITE: can''t open output file'  return, 0 }
 printf, lun, 'P5'
 nx = dimen(x,0)	ny = dimen(x,1)
 s=istring(nx,1)+' '+istring(ny,1)	;just to keep compact
 printf, lun, s
 printf, lun, '255'			;this is the max value
 filewrite, lun, x
 close, lun
 return, 1
 endfunc
 ;========================================================
func wrt_sgi(im, file)
 ;this writes out an image in SGI format
 nx=dimen(im,0)	 ny=dimen(im,1)
 lun =get_lun()
 if openw(lun ,file) eq 0 then {
	 ty,'WRT_SGI: can''t open output file'  return, 0 }
 hdr=zero(intarr(256)
 hdr(0)  = 474 ; magic number
 hdr(1)  = 1   ; type
 
 hdr(2)  = 2
 
 hdr(3)  = nx    ; xdim
 hdr(4)  = ny    ; ydim
 hdr(5)  = 1	; zdim
 
 hdr(9)  = max(im) ; range
 hdr(7)  = min(im)
 
 filewrite, lun, hdr
 filewrite, lun, im
 close, lun
 
 return, 1
 end
 ;==============================================================
func wrt_sgi_rgb(imr, img, imb, file)
 ;this writes out an image in SGI format
 nx=dimen(imr,0)	 ny=dimen(imr,1)
 lun =get_lun()
 if openw(lun ,file) eq 0 then {
	 ty,'WRT_SGI: can''t open output file'  return, 0 }
 hdr=zero(intarr(256)
 hdr(0)  = 474 ; magic number
 hdr(1)  = 1   ; type
 
 hdr(2)  = 3
 
 hdr(3)  = nx    ; xdim
 hdr(4)  = ny    ; ydim
 hdr(5)  = 3	; zdim
 
 hdr(9)  = 255 ; range
 hdr(7)  = 0
 if $little_endian then swab, hdr
 filewrite, lun, hdr
 filewrite, lun, imr
 filewrite, lun, img
 filewrite, lun, imb
 close, lun
 return, 1
 endfunc
 ;==============================================================
func yuvwrite_gray(x, name)
 ;writes YUV files for MPEG, returns 1 if OK
 ;these are 3 files, Y is the luminance, etc; we have b/w so we
 ;just write the data in Y and fill U and V with 128, note that
 ;the U and V files are sampled every other pixel so they are x4
 ;smaller; this means that nx and ny are forced even
 if symdtype(x) ne 0 or num_dim(x) ne 2 then {
 	ty,'YUVWRITE: input must be a 2-D I*1 array'  return, 0 }
 lun = get_lun()
 if openw(lun ,name+'.Y') eq 0 then {
 	ty,'PGMWRITE: can''t open output Y file'  return, 0 }
 nx = dimen(x,0)	ny = dimen(x,1)
 if nx%2 ne 0 or ny%2 ne 0 then {
 ;have to take an even extraction
 filewrite, lun, x(0:(nx-2), 0:(ny-2)) } else { filewrite, lun, x }
 close, lun
 ;now create "fake" U and V files
 xq = zero(bytarr(nx/2, ny/2)) + byte(128)
 lun = get_lun()
 if openw(lun ,name+'.U') eq 0 then {
 	ty,'PGMWRITE: can''t open output U file'  return, 0 }
 filewrite, lun, xq
 close, lun
 lun = get_lun()
 if openw(lun ,name+'.V') eq 0 then {
 	ty,'PGMWRITE: can''t open output V file'  return, 0 }
 filewrite, lun, xq
 close, lun
 return, 1
 endfunc
 ;========================================================
func yuvwrite(x, name, redmap, greenmap, bluemap)
 ;writes YUV files for MPEG, returns 1 if OK
 ;note that x could be modified if there are odd valued dimensions
 ;this is the false color version
 ;redmap, greenmap, and bluemap must be 256 arrays
 ;these are 3 files, Y is the luminance, note that
 ;the U and V files are sampled every other pixel so they are x4
 ;smaller; this means that nx and ny are forced even
 if symdtype(x) ne 0 or num_dim(x) ne 2 then {
 	ty,'YUVWRITE: input must be a 2-D I*1 array'  return, 0 }
 lun = get_lun()
 if openw(lun ,name+'.Y') eq 0 then {
 	ty,'PGMWRITE: can''t open output Y file'  return, 0 }
 nx = dimen(x,0)	ny = dimen(x,1)
 if nx%2 ne 0 or ny%2 ne 0 then {
 ;have to take an even extraction
 x = x(0:(nx-2), 0:(ny-2))
 }
 y = byte(.299 * redmap(x) + .587 * greenmap(x) + .114 * bluemap(x) )

 filewrite, lun, y
 close, lun
 ;now create real U and V files, first make a compressed version of x
 
 xc = compress(x,2)
 
 lun = get_lun()
 if openw(lun ,name+'.U') eq 0 then {
 	ty,'PGMWRITE: can''t open output U file'  return, 0 }
 y = byte(128+word(-.1687*redmap(xc)-.3313*greenmap(xc)+.5*bluemap(xc)))
 filewrite, lun, y	;write U image
 close, lun
 lun = get_lun()
 if openw(lun ,name+'.V') eq 0 then {
 	ty,'PGMWRITE: can''t open output V file'  return, 0 }
 y = byte(128+word(0.5*redmap(xc)-.4187*greenmap(xc)-.0813*bluemap(xc)))
 filewrite, lun, y	;write V image
 close, lun
 return, 1
 endfunc
 ;========================================================
subr ippinspect(name, h, params)
 lun = get_lun()
 params = zero(lonarr(11)
 params(0) = -1
 if openr(lun ,name) eq 0 then { ty,'ippread: no such file'  return, 0 }
 ;this is to be safe, just read 4 bytes
 fileread, lun, xq, 0, 4, 0
 if smap(xq) ne 'ncaa' then { ty,'not a ipp file'  close, lun  return, 0 }
 ;OK, rewind it and do cr terminated reads
 rewindf, lun
 s=''
 ;be careful, read at most 100 lines
 nc = 100
 gotem = 0
 h = ''
 while 1 {
 nc = nc - 1
 if nc le 0 then { ty,'problem finding end of header of ipp file'
 	close, lun  return, 0 }
 readf, lun, s
 h = h + s + smap(byte(10))	;clumsy way to get the header
 if num_elem(s) le 0 then { ty,'bad line in ipp header' close, lun return, 0 }
 if s(0) eq smap(byte(12)) then break	;we are done
 ty,s
 if strpos(s,'size') eq 0 then {
 ;found the size parameter
 ss = strskp(s,'=')
 nx = fix(strtok(ss,' '))
 ny = fix(strtok('','; '))
 gotem = gotem + 1
 }
 if strpos(s,'bits') eq 0 then {
 ;found the bits parameter
 ss = strskp(s,'=')
 nb = fix(strtok(ss,' '))
 ty,'nb =', nb
 if nb ne 8 and nb ne 16 and nb ne 32 then { ty,'bad bits parameters in ipp file =', nb
 	close, lun return, 0 }
 gotem = gotem + 1
 }
 }
 ty,'gotem =', gotem
 if gotem lt 2 then { ty,'could not find enough info in ipp header'
 	close, lun return, 0 } 
 params(2) = 2
 params(3) = nx
 params(4) = ny
 if nb eq 8 then { type = 0 } else if nb eq 16 then { type = 1 } else
    if nb eq 32 then type = 3
 params(1) = type
 close, lun
 params(0) = filesize(name)
 endsubr
 ;========================================================
func ippread(x, name, h)
 lun = get_lun()
 if openr(lun ,name) eq 0 then { ty,'ippread: no such file'  return, 0 }
 ;this is to be safe, just read 4 bytes
 fileread, lun, xq, 0, 4, 0
 if smap(xq) ne 'ncaa' then { ty,'not a ipp file'  close, lun  return, 0 }
 ;OK, just read the whole thing here (different approach than in
 ;ippinspect
 n = filesize(name)
 if n gt 32000000 then { ty,'this looks too big for an ipp file, size =', n
 	close, lun	return, 0 }
 fileread, lun, xx, 0, n, 0
 close, lun
 ;find the FF, assume in first 1024 bytes (easier in C)
 iq = where(xx(0:1023) eq byte(12))
 if isscalar(iq) then { ty,'can''t find end of ipp header' return, 0 }
 iq = iq(0)
 h = smap(xx(0: (iq+1) ) )
 x = xx((iq+2) : *)
 xx = 0
 #nl = smap(byte(10))
 ;decode the header
 s = strskp(h, #nl+'size=')	;need to avoid rawsize
 ;found the size parameter
 nx = fix(strtok(s,' '))
 ny = fix(strtok('','; '))
 s = strskp(h, #nl+'bits=')	;get bits
 nb = fix(strtok(s,' ;'))
 if nb ne 8 and nb ne 16 and nb ne 32 then { ty,'bad bits parameters in ipp file =', nb
 	return, 0 }
 if nb eq 16 then x = wmap(x)
 if nb eq 32 then x = fmap(x)
 redim, x, nx, ny
 ;assume that these are always big endian files and swab if machine is not
 if bigendian() eq 0 then {
   ;7/26/2004 - use the builtin endian routine rather than swapb, swapl
   endian, x
   ;;if symdtype(x) eq 1 then swapb, x
   ;;if symdtype(x) eq 2 then swapl, x
   ;;if symdtype(x) eq 3 then swapl, x
 }
 return, 1
 endfunc
 ;========================================================
func issolarbsot(x, name, h, mode)
 ;checks a fits header and returns 1 if we think this is a Solar B SOT/FPP file
 ;6/14/2010 - start work to support cfitsio files also
 nlines = dimen(h,0)
 for i = 0, nlines - 1 do {
   ;9/12/2006 - look for Solar B keys, only FG and SP have extensions at present
   if strpos(h(i),'INSTRUME= ''SOT') ge 0 then {
     ;;ty,'got a solar B extension'
     ;we assume we have a solar B file with an image header extension
     x2 = 0
     ;might be a compressed file so use rdfits if no extension case
     if mode eq 1 then iq = fits_read(x, name, h, x2, h2) else iq = rdfits(x, name, h)
     ;;iq = rdfitsex(x, x2, name, h, h2)
     ;;if iq ne 1 then return, iq
     solarboptions, x, h, x2
     $instrument_type = 3	;to be compatible with $view_data_type
     return, 1		;return 1 here even if fits_read did not
   }
 }
 return, 0
 endfunc
 ;========================================================
func rdfits_mdi_trace(x, name, h)
 ;9/12/2006 - check for image header extensions
 ;7/26/2004 - start adding checks for Solar B FITS files
 ;4/28/99 tries to detect the strange multi-image files and sets
 ;$multi_image_flag if it seems appropiate, in this case we return x=0
 ;specialized to recognize some cropped MDI data
 ;added the test trace files 10/20/97
 ;3/14/98 (pi day)
 ;because of the TRACE 3-D hourly files, it is now dangerous to read a fits file
 ;without checking first, even the extension can be rather large
 iq = fits_header(name, h, h2)
 if iq ne 1 then return, 0   ;we have to have a header
 
 ;if no extension, we will assume a civilized FITS file
 if isscalar(h2) eq 0 then {
   ;has an extension, look for some suspicious keys
   ;6/10/2010 - for compressed fits files ala cfitsio style, these are done as a binary table
   ;extension, we look for an NAXIS = 0 in h as the first hint
   iq = fits_key(h,'NAXIS')
   if iq le 0 then {
     ;6/21/2010 - this is not completely general yet, we don't look for
     ;a second extension and assume the header we want is in h2, we also don't
     ;handle Solar B or MDI, we assume it must be AIA or HMI or unknown type
     h = h2
     iq = 0
     iq = cfitsio_read(x, name)
     if iq ne 1 then { ty,'error reading compressed file'  x = 0 return, 0 }
     ;try to get times and positions, this also gets BSCALE and BZERO if they exist
     fits_time_decode, h
     return, 1
   }
   ;should run through the keys just once with a more efficient code
   if issolarbsot(x, name, h, 1) then return, 1
   nlines = dimen(h,0)
   for i = 0, nlines - 1 do {
     if strpos(h(i),'TR_WRT_FITS') gt 0 or strpos(h(i),'TR_REFORMAT') gt 0 then {
       $multi_image_flag = 1    x=0   return, 1 }
   }
   ;has an extension but apparently not a known type of multi-image file
   ;the old single image TRACE files have extensions so check first for that
   iq = rdfitsex(x, x2, name, h, h2)
   if iq ne 1 then return, iq
   if isarray(x2) then {
     ;could be, probably not a cropped MDI at least so we'll just return from here
     ;if we don't need to decode a compressed TRACE file, note that trace_decode
     ;will modify the x array if the data was compressed
     return, trace_decode(x, x2, name)
   }
 }

 ;assuming ordinary FITS file (no extension)
 ;check for a Solar B file again to catch ones w/o an extension
 if issolarbsot(x, name, h, 0) then return, 1
 $multi_image_flag = 0
 iq = rdfits(x, name, h)
 ;read in the data as a simple FITS files
 ;5/16/2010 for early AIA level 1 files, need to take out 0x8000000 entries
 if symdtype(x) eq 2 then zeronondata, x
 ;try to get times and positions, this also gets BSCALE and BZERO if they exist
 fits_time_decode, h
 ;check if x is 1-D
 if num_dim(x) eq 1 then {
   ;is 1-D, check if there is a DPC field, likely a cropped MDI image

   if ($data_dpc and 0xffff) eq 0x0fc0 then {
     ;a good case to de-crop
     x = uncrop_mdi(x)
     return, iq
   }
   ;no DPC, assume really just a 1-D array
 }
 ;if a fdmg file, check if we are running with widgets and perform some
 ;ops on the fdmg if certain options set
 if $data_fdmg_flag then {
   ;this has to be safe to run without the browser
   if defined($fdmg_rotate) then {
     if xmtogglegetstate($fdmg_scale_to_gauss(1)) then {
       if $bscale eq 0 then {
	 ;check if bunit already set to gauss (3/5/2010), if so, assume lack of bscale means in gauss
	 if strpos(upcase($bunit),'GAUSS') lt 0 then {
           errormess,'can''t scale to gauss\nNo BSCALE available\nor not a FITS file'
	 }
       } else {
         x = x * $bscale + $bzero
       }
   }
     if xmtogglegetstate($fdmg_rotate(1)) then {
      if $data_p_angle ne 0.0 then {
       ;I think we need to rotate about $data_x0, $data_y0
       ;a nearest neighbor rotate is done to avoid getting a range
       ;of "non-valid data" values from interpolation, this confuses
       ;the scaling code which will drop the lowest single value if it is
       ;less than -10000 but doesn't work if there is a spread
       x = rotate(x,-$data_p_angle, $data_x0,$data_y0)
      }
     }
     ;the mask is done after a possible rotate or else we have zeroes in
     ;the corners
     if xmtogglegetstate($fdmg_mask_limb(1)) then {
      ;we assume dataread has already gotten the parameters
      if $data_r_sun gt 0 then {
       xin = (indgen(fltarr(dimen(x,0))) - $data_x0)^2
       yin = (indgen(fltarr(dimen(x,1))) - $data_y0)^2
       redim, xin, dimen(xin, 0), 1
       redim, yin, 1, dimen(yin, 0)
       r2 = xin + yin
       xin = 0   yin = 0
       mask = sieve( r2 gt ($data_r_sun + 1)^2 )
       x(mask) = min(x)  ;make sure it is black
      }
     }
   }
 } else {
   ;for generic fits images, apply $bscale and $bzero if present, this will
   ;also convert to F*4
   if $bscale ne 0 then  { if $bscale ne 1 or $bzero ne 0 then x = x * $bscale + $bzero }
 }
 return,iq
 endfunc
 ;=======================================================================
subr get_table_ids(xx, qt, ht)
 ;get the qt and ht numbers from the ext array
 sai = xx(413)
 if $little_endian then { swab, sai }
 ncase sai
  xq = xx(294+8)
  xq = xx(294+9)
  xq = xx(294+12)
  else: ty, 'illegal subarea, sai =', sai
 endcase
 if $little_endian then { swab, xq }
 qt=xq and 0x3f   ht = (xq and 0xc0)/64
 endsubr
 ;=========================================================================
function trace_decode(x, x2, fname)
 xx = wmap(x2)
 ;if x is 1-D it should be compressed but also check trace_type
 trace_type = xx(2)
 if trace_type ne 1 then {
 ;probably uncompressed and already unpacked if it was packed
 if num_dim(x) ne 2 then {
 	ty,'possible error, data is 1-D but not compressed type?'
 }
 return, 1	;return the data as is
 }
 ;here is compressed
 get_table_ids, xx, qt, ht
 ty,'qt, ht =', qt, ht
 ;qt =' /hosts/shimmer/usr/people/shine/trace/qtables/qt'+string(qt)
 ;ht =' /hosts/shimmer/usr/people/shine/trace/hufpack1'
 path = getenv('TOPS_IDS_PATH')
 if path eq '' then {
  ;a default if the TOPS_IDS_PATH is not defined
  qt =' /tsw/dbase/tables/qt'+string(qt)
  ht =' /tsw/dbase/tables/ht0'
  } else {
  qt =' ' + path + '/tables/qt'+string(qt)
  ht =' ' + path + 'tables/ht0'
 }
 sq='/tsw/ana/trace_decode '+fname+' 0 junk_decode '+ht+qt+' > decode.record'
 ;ty,sq
 spawn, sq
 ;now read it back from junk
 iq = rdfits(x, 'junk_decode', h)
 ;and remove junk
 spawn, 'rm junk_decode'
 return, iq
 endfunc
 ;========================================================
subr solarboptions, x, h, x2
 ;this is now called by dataread, so it has to work standalone
 ;parse some keywords that we are likely to need and stash in globals
 fits_time_decode, h
 ;also prepare a text entry, needed in browser
 solarbkeystrs, h, sobs	;this also sets $sbtype
 $solarbstrs = '\n' + sobs
 if $sbtype le 1 then {
   ;keep originals for C code approach
   $fg_ix0 = fix(fits_key(h, 'FGCCDIX0'))
   $sb_ix1 = 4095 - $fg_ix0 ;note reversal of ix0/1
   $fg_ix1 = fix(fits_key(h, 'FGCCDIX1'))
   $sb_ix0 = 4095 - $fg_ix1
   $fg_iy1 = fix(fits_key(h, 'FGCCDIY1'))
   $sb_iy0 = 2047 - $fg_iy1
   $fg_iy0 = fix(fits_key(h, 'FGCCDIY0'))
   $sb_iy1 = 2047 - $fg_iy0
   $sb_genid = fix(fits_key(h, 'GEN_ID'))
   $sb_sum = fix(fits_key(h, 'CAMSSUM'))
   $sb_bin = fix(fits_key(h, 'FGBINX'))
   ;figure ix0, ix1 for summed or binned pixels
   ;(not sure why I do this, it simplifies some things but makes others more awkward)
   if $sb_sum gt 1 or $sb_bin gt 1 then {
     iq = $sb_sum * $sb_bin
     $sb_ix0 = $sb_ix0/iq
     $sb_ix1 = ($sb_ix1+iq-1)/iq -1
     $sb_iy0 = $sb_iy0/iq
     $sb_iy1 = ($sb_iy1+iq-1)/iq -1
   }
   $sb_roistop = fix(fits_key(h, 'ROISTOP'))
 }
 ;for FG's, some image processing related to pedestal and central gap
 if notdefined($sbhasimage) then {
   compile_file, getenv('ANA_WLIB') + '/solarbvariables.ana'
 }
 if $sbtype le 1 then {
   zeroifundefined, $solarbpedestalflag, $solarbmidfixflag, $sbmagpatchtype, $solarbmagfixflag ,$motif_running, $roiloopflag
   ;however, we might be reading a "processed" solar B file, so abort if we have any history
   ;entries, obviously clumsy
   ;also assume processed if not I*2
   if symdtype(x) eq 1 then {
     if fits_key_string(h, 'HISTORY') eq 'undefined' then {
       ;;ty,'running patches'
       if $solarbmidfixflag then solarbpatchgap, x
       ;do the solarbpatchgap first since it will convert I's to FP assuming unsigned I*2, otherwise
       ;the pedestal would add for data over 2^15
       if $solarbpedestalflag then solarbpedremove, x
       if $solarbmagfixflag then solarbmagpatch, x	;target will determine if a V to fix
     }
   }
 }
 ;put the image header into a global after a few checks
 if num_dim(x2) eq 1 then { exchange, x2, $sbimagehdr } else {
   $sbimagehdr = 0
   ;;ty, 'something wrong with Solar B image header? (missing extension?)'
   ;;errormess,'something wrong with\nSolar B image header?\n(missing extension?)'
 }
 ;FG or SP or something else?
 if $sbtype eq 0 then $sbtype = 1
 ncase $sbtype-1
  { 
   ;a FG
   sq = getsbwave(h)
   ;put wavelength on next line since obs type is sometimes rather long
   $solarbstrs += ', ' + sq
   sq = fits_key(h, 'WAVEOFF')
   if isscalar(sq) then $solarbstrs = $solarbstrs + sprintf('\n%dmA', sq)
   ;if actual exposure time is available, use it, otherwise
   ;default to a commanded one (there are some older fits files
   ;with the MEXPOSE keyword, these should age off eventually but
   ;support for a while)
   sq = fits_key(h, 'EXP0')
   if isscalar(sq) eq 0 then sq = fits_key(h, 'EXPTIME')
   if isscalar(sq) eq 0 then sq = fits_key(h, 'MEXPOSE')
   if isscalar(sq) then $sb_exptime = sq else { $sb_exptime = 0.0  sq = 0.0 }
   ;but don't use if shutterless
   ;note, already have gen ID from dataread
   ;;genid = fits_key(h, 'GEN_ID')
   if isscalar($sb_genid) then { if $sb_genid lt 32 then
      $solarbstrs += sprintf(',%5.2fs', sq) }
   ;and the focus
   sq = fits_key(h, 'FOCUS')
   if isscalar(sq) then { sq = sq - 2048
   $solarbstrs += sprintf(', focus:%4d\n', sq) }
   ;also do the image compression type (e.g., JPEG) and the bits/pixel
   getcompstrings, h, c1str, c2str
   $solarbstrs += c1str + ' ' + c2str
   ;;ty, 'our string = ', $solarbstrs
  }
  {
   ;a SP
   ;;$lastsp = $list_item_position
   $lastsphead = h
   sq = fits_key(h, 'SLITPOS')
   if isscalar(sq) then $solarbstrs = $solarbstrs + sprintf('\nslit %d', fix(sq))
   sq = fits_key(h, 'FOCUS')
   if isscalar(sq) then { sq = fix(sq) - 2048  $solarbstrs += sprintf(', focus:%4d', fix(sq)) }
   getcompstrings, h, c1str, c2str
   $solarbstrs += '\n' + c1str + ' ' + c2str
  }
 endcase

 endsubr
 ;========================================================
subr solarbkeystrs, h, sobs
 ;gets some keyword stuff that is used in several places
 sobs = fits_key(h,'OBS_TYPE')
 if strcontains(upcase(sobs), 'UNDEFINED') then {
   ;might be an SP w/o the OBS_TYPE keyword so try PRODUCT
   sobs = fits_key(h, 'PRODUCT')
 }
 ;also add camera summing and/or binning info here, assume square for FG
 iq = fits_key(h, 'CAMSSUM')
 if isscalar(iq) then {
   if iq ge 2 then sobs = sobs + sprintf('[S%d]', iq)
 }
 iq = fits_key(h, 'FGBINX')
 if isscalar(iq) then {
   if iq ge 2 then sobs = sobs + sprintf('[B%d]', iq)
 }
 if strcontains(sobs,'FG') then $sbtype = 1 else
   if strcontains(sobs,'SP') then $sbtype = 2 else
     if strcontains(sobs,'AC') then $sbtype = 3 else
       if  strcontains(sobs,'CT') then $sbtype = 4 else
         if  strcontains(sobs,'jitter') then $sbtype = 5 else
           if  strcontains(sobs,'analog') then $sbtype = 6 else
             if  strcontains(sobs,'Scan') then $sbtype = 7 else
         $sbtype = 0
 ;the INSTRUME check is only needed for older data, once these are
 ;re-re-formatted, we can drop
 sq = fits_key(h, 'INSTRUME')
 if strcontains(sq,'Filtergraph') then $sbtype = 1 else
   if strcontains(sq,'Spectro') then $sbtype = 2 else
     if strcontains(sq,'Analog') then $sbtype = 3 else
       if strcontains(sq,'Correlation') then $sbtype = 4
 ;;ty,'in solarbkeystrs, $sbtype,sobs =', $sbtype,'  ', sobs
 endsubr
 ;========================================================
func getsbwave(h)
 ;check if a dark, if so the wavelength isn't too important
 sq = fits_key(h, 'DARKFLAG')
 ;also keep in a global since we need to know later
 $darkflag = 0
 if isscalar(sq) then { if sq eq 1 then {
   sq = 'DARK'
   $darkflag = 1
   return, sq
   }
 }
 
 sq = fits_key(h, 'WAVE   ')
 ;if this is a no move, try to get filter wheel position and decode
 if strcontains(sq, 'no m') then {
   iok = 0
   if strcontains(sq, 'NFI') then {
     iq = fits_key(h, 'NBFW')
     if isscalar(iq) then {
       for i=0,5 do { if abs(iq-$sbnbfw(i,1)) lt 3 then { iok = 1  break } }
       if iok then sq = sprintf('NB, at %d', $sbnbfw(i,0)) else
       	sq =sprintf('NB, wheel=%d?', iq)
     }
   } else {
     iq = fits_key(h, 'WBFW')
     if isscalar(iq) then {
       for i=0,5 do { if abs(iq-$sbwbfw(i,1)) lt 3 then { iok = 1  break } }
       if iok then sq = sprintf('WB, at %d', $sbwbfw(i,0)) else
       	sq =sprintf('WB, wheel=%d?', iq)
       
     }
   }
 }
 ;some of these are rather long, remove some stuff here
 sq = strreplace(sq,'TF ', '')
 sq = strreplace(sq,'cont ', '')
 sq = strreplace(sq,'bandhead ', '')
 return, sq
 endfunc
 ;========================================================
func strcontains(s, test)
 ;returns 1 if s contains test and 0 otherwise
 if isstring(s) and isstring(test) then {
   if strpos(s, test) ge 0 then return, 1
 } 
 return, 0
 endfunc
 ;========================================================
func strcontainstype(s, stest)
 ;returns index of match if a path s contains one of the strings in stest
 ;as a path of the path and -1 otherwise
 if isstring(s) ne 1 then return, -1
 if isstrarr(stest) then {
   n = num_elem(stest)
   for i=0,n-1 do {
     if strpos(s, '/'+stest(i)+'/') ge 0 then return, i
   }
   ;no match (or we wouldn't get here) so return -1
   return, -1
 } else if isstring(stest) then {
   if strpos(s, '/'+stest+'/') ge 0 then return, 0 else return, -1
 } else return, -1
 endfunc
 ;========================================================
func strcontainsone(s, stest)
 ;returns index of first match if s contains one of the strings in stest
 ;and -1 otherwise, works best if no strings are subsets of others 
 if isstring(s) ne 1 then return, -1
 if isstrarr(stest) then {
   n = num_elem(stest)
   for i=0,n-1 do {
     if strpos(s, stest(i)) ge 0 then return, i
   }
   ;no match (or we wouldn't get here) so return -1
   return, -1
 } else if isstring(stest) then {
   if strpos(s, test) ge 0 then return, 0 else return, -1
 } else return, -1
 endfunc
 ;========================================================
subr getcompstrings(h, s1, s2)
 ;decode the compression parameters from a Solar B FITS header
 if defined($sbqvalue) eq 0 then $sbqvalue = ['Q98','Q90','Q75','Q50','Q95','Q92','Q85','Q65']
 s1 = ''	s2 = ''
 ;;sq = fits_key(h, 'BITCOMP1')
 ;;if isscalar(sq) then s1 = sprintf('L%d', sq)
 ;;sq = fits_key(h, 'IMGCOMP1=')
 ;;if isscalar(sq) then s1 = s1 + $sbctype(sq)
 sq = fits_key(h, 'QTABLE1')
 if isscalar(sq) then s1 = s1 + $sbqvalue(sq)
 sq = fits_key(h, 'BITSPPI')
 if isscalar(sq) then s1 = s1 + sprintf(' %4.1fbpp ', sq)
 sq = fits_key(h, 'QTABLE2')
 if isscalar(sq) then s2 = s2 + $sbqvalue(sq)
 sq = fits_key(h, 'BITSPPQ')
 if isscalar(sq) then s2 = s2 + sprintf(' %4.1fbpp ', sq)
 endsubr
 ;========================================================
subr solarbpedremove, x
 ;this doesn't work properly for sums, only single exposures
 ;remove pedestal from a FG image, first get some keys
 ;2/17/2009 - modified to avoid V's and get I's in 3-D products
 ;note that fits_key is not safe for all keys
 ;we only can remove from an I image, if there is one, we assume first in stack
 if $sb_genid ge 32 then iflag = 1 else iflag = $sbhasimage($sb_genid < 18)
 if iflag eq 0 then return   
 nx = dimen(x,0)
 if $sb_sum eq 1 then {
   ps = $sbpedestal($sb_ix0:$sb_ix1,*)
 } else if $sb_sum eq 2 {
   ps = $sbpedestals2($sb_ix0:$sb_ix1,*)
 } else if $sb_sum eq 4 then {
   ty,'not handling 4x4 summing yet in solarbpedremove'
   return
 }
 ps = word(ps)
 redim, ps, dimen(ps,0), 1
 ;following added to handle processing cases for compound images outside of browser
 if defined($knownnumberofi) then ps = ps * $knownnumberofi
 ;we now have a vector with the pedestal varying in x direction only
 if num_dim(x) eq 3 then {
   ;we already checked for existence of I
   x(0,0,0) = x(*,*,0) - ps
 } else  if num_dim(x) eq 2 then {
   x = x - ps
 } else { ty, 'unexpected dimension for data?' }
 return
 endsubr
 ;========================================================
subr solarbpatchgap, x
 if symdtype(x) ne 1 then x = short(x)
 ;patch the central gap if appropiate for this image (must be a FG)
 ;2/17/2009 - extend to non-I cases, might be some problems still, start using fgccdshift
 ;;if num_dim(x) gt 2 then { ty,'patch CCD center: we don''t do 3-D cases yet' return }
 if $sb_sum eq 4 then { ty,'not handling 4x4 summing yet in "patch center"'  return }
 if $sb_bin gt 1 then { ty,'not handling binning summing yet in "patch center"'  return }
 ;we also have to give up if this is a shutterless with ROI looping
 if $roiloopflag then {
   ;check if it might have ROI looping, it must be shutterless of course
   if $sb_genid ge 32 then {
     ;worth the expense of a fits keyword check
     roiloop = fits_key($sbhead,'ROILOOP')
     if isscalar(roiloop) eq 0 then {
	errormess, 'no ROILOOP keyword in FITS file?\npatch center disabled'
	return }
     if roiloop ge 1 then { ty,'not doing "patch center" for ROI loops yet'  return }
   }
 }
 if num_dim(x) eq 2 then {
   fgccdshift, x, $sb_roistop, $fg_ix1, $fg_iy1, $sb_sum, $sb_bin, $sb_genid
 } else if num_dim(x) eq 3 then {
   ;;until we make fgccdshift smarter, we have to extract and reinsert each 2-D image
   nd = dimen(x, 2)
   for id = 0, nd-1 do {
     xq = x(*,*,id)
     fgccdshift, xq, $sb_roistop, $fg_ix1, $fg_iy1, $sb_sum, $sb_bin, $sb_genid
     x(0,0,id) = xq
   }
   ;if we have an I image in a 3-D, we have to float now to make any I signal a
   ;unsigned quantity, because the 3-D keeps getting stuff inserted, we have to
   ;do it before this can cause unintended casts between short and float in operations
   ;that follow, needs to be done more efficiently someday
   if $sb_genid ge 32 then iflag = 1 else iflag = $sbhasimage($sb_genid < 18)
   if iflag then {
     x = float(x)
     y = x(*,*,0)
     y = y + 65536. * (y lt 0)
     x(0,0,0) = y
   }
 }
   
//      }
//      ;insert interpolated columns, special cases abound but we have already
//      ;returned if there are no middle columns so we know we have at least one
//      ;and if there is one and not the other, their weight sum ws has already been
//      ;adjusted
//      x1 = rfix((2.0*xa + xb) * ws)
//      x2 = rfix((2.0*xb + xa) * ws)
//      redim, x1, 1, ny
//      redim, x2, 1, ny
//      if im ge 0 and im le nxm1 then insert, x, x1, im
//      if imp1 ge 0 and imp1 le nxm1 then insert, x, x2, imp1
//    }
//  return		;end of the $sb_roistop lt 1025 case
//  }
//  
 ;as a temporary fix for 2x2 partial reads, back out 1/2 pixel of horizontal shift
 ;this should be modified to avoid interpolation (causing a bias to one side instead)
 ;and integrated in fgccdshift
 if $sb_roistop lt 1025 then {
   if $sb_sum eq 2 then {
     mid = 1023 - $sb_ix0
     if mid gt 0 and mid lt $sb_ix1 then {
       if num_dim(x) eq 2 then {
	 z1 = x(0:mid,*)
	 z2 = z1
	 insert, z2, z1(0:(mid-1),*),1,0
	 z1 = 0.5*(z1+z2)
	 z2 = x((mid+1):*,*)
	 insert, x, z1, 0, 0
	 z1 = z2
	 insert, z1, z2(1:*, *),0,0
	 z1 = 0.5*(z1+z2)
	 insert, x, z1, mid+1, 0
       } else if num_dim(x) eq 3 then {
         ;do the stack, keep it simple for now, optimize when included in fgccdshift
            for id = 0, nd-1 do {
	      xq = x(*,*,id)
	      z1 = xq(0:mid,*)
	      z2 = z1
	      insert, z2, z1(0:(mid-1),*),1,0
	      z1 = 0.5*(z1+z2)
	      z2 = xq((mid+1):*,*)
	      insert, xq, z1, 0, 0
	      z1 = z2
	      insert, z1, z2(1:*, *),0,0
	      z1 = 0.5*(z1+z2)
	      insert, xq, z1, mid+1, 0
	      x(0,0,id) = xq
	  }
       }

     }
   }
 }
 if $sb_roistop eq 1025 then {
   ;if not a partial read there are 2 bright central columns that we
   ;are approximately patching for some special cases
   ;we only can remove from an I image, if there is one, we assume first in stack
   if $sb_genid ge 32 then iflag = 1 else iflag = $sbhasimage($sb_genid < 17)
   if iflag then {
     if num_dim(x) eq 3 then {
       ;we already checked for existence of I
       exchange, x, xsave
       x = xsave(*,*,0)		;x is now just the I part
     } else xsave = 0

     if $sb_sum eq 2 then {
       ;check if the 1023 column is in the image
       if $sb_ix0 le 1023 and $sb_ix1 ge 1023 then {
	 inx = 1023 - $sb_ix0
	 c1 = x(inx,*)*0.912
	 redim, c1, 1, num_elem(c1)
	 insert, x, c1, inx, 0
       }
       ;check if the 1024 column is in the image
       if $sb_ix0 le 1024 and $sb_ix1 ge 1024 then {
	 inx = 1024 - $sb_ix0
	 c1 = x(inx,*)*0.921
	 redim, c1, 1, num_elem(c1)
	 insert, x, c1, inx, 0
       }
     }


     if $sb_sum eq 1 then {
       ;check if the 2047 column is in the image
       if $sb_ix0 le 2047 and $sb_ix1 ge 2047 then {
	 inx = 2047 - $sb_ix0
	 c1 = x(inx,*)*0.84
	 redim, c1, 1, num_elem(c1)
	 insert, x, c1, inx, 0
       }
       ;check if the 2048 column is in the image
       if $sb_ix0 le 2048 and $sb_ix1 ge 2048 then {
	 inx = 2048 - $sb_ix0
	 c1 = x(inx,*)*0.852
	 redim, c1, 1, num_elem(c1)
	 insert, x, c1, inx, 0
       }
     }
     ;need to check if we just corrected a 3-D
     if num_dim(xsave) eq 3 then {
       x = short(x)
       xsave(0,0,0) = x
       exchange, x, xsave
     } 
   }
 }
 endsubr
 ;========================================================
subr solarbmagpatch, x, forceflag
 ;check if a V image in here unless we have the forceflag argument, then we just assume
 ;a 2-D array and apply (checks if 2-D)
 if !narg lt 2 then {
   if $sb_genid gt 40 then vflag = 0 else vflag = $sbhasv($sb_genid)
 } else {
   ;must be 2-D
   if num_dim(x) ne 2 then { ty,'solarbmagpatch must have 2-D array if force flag is set, no action'  return }
   vflag = 1
 }
 if vflag then {
   isum = 1
   if $sb_sum gt 1 or $sb_bin gt 1 then { isum = $sb_sum * $sb_bin }
   if isum gt 4 then { ty,'not handling total sum of ', isum, ' yet in "mag fix"'  return }
   vflag = vflag - 1
   ;;ty, 'have a V at index:', vflag
   if isum eq 2 then {
     ;check if the 1023 column is in the image
      if $sb_ix0 le 1023 and $sb_ix1 ge 1023 then {
        inx = 1023 - $sb_ix0
	w = $magdecay
	} else return
   } else if isum eq 1 then {
     ;check if the 2047 column is in the image
     if $sb_ix0 le 2047 and $sb_ix1 ge 2047 then {
        inx = 2047 - $sb_ix0
	w = 2.0 * $magdecay
	} else return
   } else if isum eq 4 then {
     ;check if the 511 column is in the image
     if $sb_ix0 le 511 and $sb_ix1 ge 511 then {
        inx = 511 - $sb_ix0
	w = 4.0 * $magdecay
	} else return
   }
   inx = inx + 1
   nx = dimen(x,0)
   f = zero(fltarr(nx))
   if num_dim(x) eq 3 then {
     ;we already checked for existence of V
     exchange, x, xsave
     x = xsave(*,*,vflag)		;x is now just the V part
   } else xsave = 0
   xx = indgen(fltarr(inx)) - inx
   ;estimate h, this may not be the final method
   h = mean(x((inx-5):(inx-1), *)) - mean(x(inx:(inx+4), *))
   ;;ty,'jump =', h
   s1 = x((inx-5):(inx-1), *)
   s2 = x(inx:(inx+4), *)
   ;for each, take a sdev
   inq = sieve(abs(s1) lt 2.0*sdev(s1))
   if isarray(inq) then s1 = mean(s1(inq)) else s1 = mean(s1)
   inq = sieve(abs(s2) lt 2.0*sdev(s2))
   if isarray(inq) then s2 = mean(s2(inq)) else s2 = mean(s2)
   hw3 = s1 - s2
   ;;ty,'hw3 =', hw3
   ;if $sbmagpatchtype eq 1 then h = hw3
   ;just take hw3 for now
   h = hw3
   f(0) = h * (1.0+$magbounce) * exp(xx/w) - h * $magbounce
   redim, f, nx, 1
   x = x - f
   if num_dim(xsave) eq 3 then {
     x = short(x)
     xsave(0,0,vflag) = x
     exchange, x, xsave
   } 
 }
 endsubr
 ;=============================================================================
subr zeronondata, x
 ;used to zero the pesky 0x80000000's in HMI data
 ix = sieve(x eq 0x80000000)
 if isarray(ix) then {
   x(ix) = 0
   ix = 0
 }
 endsubr
 ;========================================================