pro sp_make_cube_index, dstring,                    $ ;example: dstring = '20140220_233311'. required
                        wmap=wmap,                  $ ;output, written to fits directory.  An index of bad and interpolated columns.
                        nbad=nbad,                  $ ;passed through to sp_quick_cube.
                        outdir=outdir,              $ ;must be absolute path (not relative) otherwise movie-making will not work
                        quiet=quiet,                $ ;
                        fits_only=fits_only,        $ ; only generate the fits output products, no images or movies.
                        ffmpeg_path=ffmpeg_path,    $ ; path to ffmpeg executable.  If not set then look for ffmpeg in the execution directory.
                        filldata=filldata,          $ ; passes through to sp_quick_cube.pro: 0 or median switch to fill bad data values
                        loud=loud,                  $ ;
                        flat_dir_tree=flat_dir_tree   ; convenience output: don't use time-heirarchical tree structure
;
;+
;    Name:               sp_make_cube_index
;
;    Purpose:            Build multiple data cube products from spectropolarimeter data and save in iris FITS format.
;                        Also, make corresponding overview image products, i.e., images or movies, as appropriate.
;                        Current data cube products:
;                         1.)  blapp/
;                         2.)  btapp/'B_transverse'
;                         3.)  conti/'Continuum Intensity'
;                         4.)  dopp2/'Doppler Shift'
;                         5.)  magfield/'Magnetic Field Strength'
;                         6.)  magincl/'Magnetic Field Inclination'
;                         7.)  magazim/'Magnetic Field Azimuth'
;                         8.)  ptot/
;                         9.)  dopp1/
;                        10.)  rufaz
;
;    Input Parameters:   dstring -- the start time/date string of a hinode spectropolarimeter observation.  Required
;
;    Keyword Parameters:
;                        wmap   - returned from sp_quick-cube.pro
;                        nbad   - index number of bad/dropped column.  These 'bad' columns will be interpolated if they fall within interpolation rules.
;                        outdir - output top level directory for completed data cube products. (optional).  Default write is to the local directory.
;                         
;
;    Output:             fits files containing multiple hinode spectropolarimeter datacube products.
;
;    Calling Sequence:
;                        IDL> sp_make_cube_index, '20140220_233311' ; <-- simple example, output directory is the current execution directory
;                        more complex calling examples from helper function: "make_some_sp_cubes.pro":
;                        sp_make_cube_index, sp_filelist(i), outdir = outdir_str, quiet=quiet, fits_only=fits_only,filldata=filldata,loud=loud, flat_dir_tree=flat_dir_tree
;
;    History:            2016-08-16 Ted Tarbell (original, version2)
;                        2016-08-18 PGS         Clean-up, modified output fits format.
;                        2020-04-22 MLD Removed calls to stks2index and used a faster method
;                        2020-05-01 MLD Now works when the B_ANGLE, P_ANGLE, and SOLAR_RA keywords are missing from the L2 header
;                        2021-07-09 MLD Now interpolates Level 2 physical quantities rather than Level 1 estimates of such quantities
;                        2021-09-03 MLD Now inserts a HISTORY record into the FITS headers of the arrays indicating the version of this code being used
;                        2021-09-16 MLD Removed extra BTYPE and BUNIT keywords from main FITS headers
;                        2022-05-26 MLD upon encountering an "oddly empty scan", this routine now redefines outdir to be outdir_orig before exiting
;-
;

version_date='2021-09-16'
  
;some notes: 
;for testing purposes:
;20181204_183006 is a 2D image with some missing data.

;some questions:
;20190419

;  currently, the following resolve_routine calls are required (by PGS) to avoid name-space conflicts with other routines 
resolve_routine, 'image2movie2' ;right now, this has to be the local version because of the string array addition fix.
resolve_routine, 'sotsp_l2struct2data' ;
resolve_routine, 'strfind_urls',/IS_FUNCTION
resolve_routine, 'ssw_ffmpeg'   ;  right now, this has to be the local version because of the string array addition fix.

width_image  = -1  ;  initial value
height_image = -1  ;  initial value 

;----------------------------- Start of quick cube generation -------------------------;

;  output directory management.  Note: default output is to the local directory.
startdir = curdir()
if ~keyword_set(outdir) then outdir = curdir()  ;  needs to be absolute path, otherwise moviemaking won't work for some reason
outdir_orig = outdir 

;  get level 2 file path (/archive1/hinode/sot/level2hao/...)
maptime = file2time(dstring)
if ~keyword_set(quiet) then print,' sp_make_cube_index ',dstring,'  ',maptime
sp2=sotsp_time2level2(maptime)
if (n_elements(sp2) gt 1) then begin
  w=where(strpos(sp2(*),dstring) ne -1,count)
  if (count gt 0) then sp2 = sp2(w)
end
if (max([strpos(sp2(*),dstring)]) eq -1) then begin  ;  return if no valid path
  if ~keyword_set(quiet) then print, dstring,'  file not found'
  return
endif
if ~file_exist(sp2) then begin  ;  return if no valid path is present
  print, 'file: ',sp2, ' not found.  Returning.'
  return
endif

;  determine level1d file path from level 2 file path
break_file,sp2,disk_log,dir,filnam,ext
dirlevel1d = str_replace(disk_log+dir,'level2hao','level1d')
if ~keyword_set(quiet) then print, dirlevel1d
lev1conti=dirlevel1d(0)+'sotsp_conti_'+filnam(0)+'.fits'  ;  path of level1d continuum FITS file
if ~file_test(lev1conti) then lev1conti=str_replace(lev1conti,'/archive','/archive1')
if ~file_test(lev1conti) then begin
  print,'  Level 1.5 continuum file for '+filnam+' not found, returning...'
  return
endif
lev1hdr=headfits(lev1conti)
index1=fitshead2struct(lev1hdr)
if ~keyword_set(quiet) then print,'successfully read level1d continuum FITS header'
if size(index1,/tname) ne 'STRUCT' then begin
  print, '  Level 1.5  continuum header for '+filnam+' is not a structure, returning...'
  return 
endif
index1 = index1(0)  ;  header to level1d continuum FITS file, mainly used to harvest keywords below

;  read in level 2 data
read_sotsp,sp2,index2,data,scan_info=scan_info,/silent,/xycoord
if ~keyword_set(quiet) then print,'successfully read level2 data file'
if n_dimensions(data) lt 3 then begin
  if ~keyword_set(quiet) then print, filnam+' ... number of dimensions of data lt 3 -- returning'
  return
endif

;  run through sotsp_rasterize
slitmap=sotsp_rasterize(scan_info)
nx=(size(slitmap,/dim))[0]
ny=index2[0].naxis2
if size(slitmap,/n_dim) eq 2 then nt=(size(slitmap,/dim))[1] else nt=1l

;  read in level 1 spectra (why?)
dirlevel1 = str_replace(disk_log+dir,'level2hao','level1hao')
l1images = file_list(dirlevel1,'*_stksimg.save')
stks=sotsp_stks2struct(l1images)

;xxxxxxxxxx

;filtering out duplicate images -- START
utimed = stks.utimed ;
if n_dimensions(utimed) eq 2 then begin
  mid_value_2nd_term = ((size(utimed))(2))/2 ;integer division -- for this purpose we want a columnar value near the middle.
  end_value_1st_term =  (size(utimed))(1)
endif else begin
  if ~keyword_set(quiet) then print, 'utimed not a 2 dimensional array.  Returning.'
  return
endelse

;here we try to remove duplicate obs by identifying and removing from the data array any that are too close together in time, which is the signature of a duplicate obs.
;the closeness-in-time threshold is defined below in DUP_OBS_DELTA_TIME_THRESH.  (Currently at .6 seconds.)
DUP_OBS_DELTA_TIME_THRESH = 0.6
;note, in a case of duplicate obs we assume the 2nd one is bad -- therefore the first element "0" of the list can never be bad, and we include it explicitely here.
;
;this is a vector of array positions of too-close-in-time observations (as found through a difference of offset-by-one times)
wbadt  = 1 + where(3600.*abs(utimed(1:(end_value_1st_term-1),mid_value_2nd_term) - utimed(0:(end_value_1st_term-2),mid_value_2nd_term)) le DUP_OBS_DELTA_TIME_THRESH, num_bad)
;this is a vector of non-duplict, i.e., not-too-close together in time) array positions as found through a different of offset-by-one times.
wgoodt = [0,1 + where(3600.*abs(utimed(1:(end_value_1st_term-1),mid_value_2nd_term) - utimed(0:(end_value_1st_term-2),mid_value_2nd_term)) gt DUP_OBS_DELTA_TIME_THRESH)]
;now, pull out only the non-nuplicate date to work with below:
data = data[wgoodt,*,*]; this should cause magfield, magincl, magazim and fillfac to be good.
;we are assuming that we only ever have 2 contiguous duplicate entries. Because of why/how they occur.  Does not handle more than two sequential dup obs. 

;  MLD notes: can add tags to existing structure using add_tag.pro
scan_info.slitpos  = scan_info.slitpos(wgoodt)
scan_info.times    = scan_info.times(wgoodt)
scan_info=add_tag(scan_info,index1.nslitpos,'nslitpos',/no_copy)
scan_info=add_tag(scan_info,index1.scn_sum,'scn_sum',/no_copy)
scan_info=add_tag(scan_info,index1.scn_rpt,'scn_rpt',/no_copy)
scan_info=add_tag(scan_info,index1.spmapctr,'spmapctr',/no_copy)
scan_info=add_tag(scan_info,index1.scn_step,'scn_step',/no_copy)
scan_info=add_tag(scan_info,index1.camssum,'camssum',/no_copy)

utimed = utimed[wgoodt,*]

;prep data arrays of non-duplicate observations to be processed by sp_quick_cubes:
;a.)
blapp = stks.blapp                  ;pull into a local variable from the common data structure
blapp = blapp[wgoodt,*]             ;remove any bad==duplicate entries
;b.)
btapp = stks.btapp                  ;pull into a local variable from the common data structure
btapp = btapp[wgoodt,*]             ;remove any bad==duplicate entries 

;c.)
conti = stks.conti                  ;pull into a local variable from the common data structure
conti = conti[wgoodt,*]             ;remove any bad==duplicate entries

;d.)
dopp2 = reform(stks.ctrline(*,*,1)) ;pull into a local variable from the common data structure
dopp2 = dopp2[wgoodt, *]            ;remove any bad==duplicate entries
;
;e.)
;; slat = stks.slat                    ;pull into a local variable from the common data structure
;; slat = slat[wgoodt,*]               ;remove any bad==duplicate entries
;
;f.)
;; slong = stks.slong                  ;pull into a local variable from the common data structure
;; slong = slong[wgoodt,*]             ;remove any bad==duplicate entries
;
;g.)
dopp1 = reform(stks.ctrline(*,*,0)) ;pull into a local variable from the common data structure
dopp1 = dopp1[wgoodt,*]             ;remove any bad==duplicate entries
;
;h.)
;; rufaz = !radeg*stks.rufaz           ;pull into a local variable from the common data structure
;; rufaz = rufaz[wgoodt,*]             ;remove any bad==duplicate entries
;
;i.)
ptot = stks.ptot                    ;pull into a local variable from the common data structure
ptot = ptot[wgoodt,*]               ;remove any bad==duplicate entries
;
index2.naxis1 = (size(wgoodt))(1)   ;resize the number of columns in the level2 index naxis1 to correspond to the prepped "good/non-duplicate" data arrays.
;
l1size = size(blapp)
if (l1size(1) ne index2(0).naxis1 or l1size(2) ne index2(0).naxis2) then begin
  if ~keyword_set(quiet) then print,'Size mismatch L1: ',l1size(1:2),'  L2: ',index2(0).naxis1,index2(0).naxis2
  return
end

if ~keyword_set(nbad) then nbad = 0 

;Build Data Cube Image products -- START

;  we can interpolate the following quantities directly
;
;3.) conti
if ~keyword_set(quiet) then box_message, ['Sending conti array to sp_quick_cube']
wh=where(index2.ftype eq 'Continuum Intensity',nwh)  ;  usually 31
if nwh eq 0 then begin
  print,'  quantity not found in level 2 data'
  return
endif
conti = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud, /filldata)
conti_sp_qc = conti
;
;11.) dopp1
if ~keyword_set(quiet) then box_message, ['Sending dopp1 array to sp_quick_cube']
wh=where(index2.ftype eq 'Doppler Shift of 630.15nm',nwh)  ;  usually 3
if nwh eq 0 then begin
  print,'  quantity not found in level 2 data'
  return
endif
dopp1 = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,/filldata,loud=loud)
dopp1 = (dopp1-median(dopp1))
;dopp1 = (dopp1-median(dopp1))*0.021549/6301.5*2.997925E5  ;  what is this factor?
dopp1_sp_qc = dopp1
;
;4.) dopp2
if ~keyword_set(quiet) then box_message, ['Sending dopp2 array to sp_quick_cube']
wh=where(index2.ftype eq 'Doppler Shift of 630.25nm',nwh)  ;  usually 4
if nwh eq 0 then begin
  print,'  quantity not found in level 2 data'
  return
endif
dopp2 = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,/filldata,loud=loud)
dopp2 = (dopp2-median(dopp2))
;dopp2 = (dopp2-median(dopp2))*0.021549/6301.5*2.997925E5  ;  what is this factor?
dopp2_sp_qc = dopp2

;
;13.) ptot
if ~keyword_set(quiet) then box_message, ['Sending ptot array to sp_quick_cube']
wh=where(index2.ftype eq 'Polarization Degree',nwh)  ;  usually 32
if nwh eq 0 then begin
  print,'  quantity not found in level 2 data'
  return
endif
ptot = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,/filldata)
ptot_sp_qc = ptot

;  for the magnetic field, better to interpolate components and then
;  back out field strength, azimuth angle, and inclination angle
if ~keyword_set(quiet) then box_message, ['Sending magfield array to sp_quick_cube']
wh1=where(index2.ftype eq 'Magnetic Field Strength',nwh1)  ;  usually 0
wh2=where(index2.ftype eq 'Magnetic Field Inclination',nwh2)  ;  usually 1
wh3=where(index2.ftype eq 'Magnetic Field Azimuth',nwh3)  ;  usually 2
wh4=where(index2.ftype eq 'Stray Light Fraction',nwh4)  ;  usually 11
if (nwh1 eq 0) or (nwh2 eq 0) or (nwh3 eq 0) or (nwh4 eq 0) then begin
  print,'  quantity not found in level 2 data'
  return
endif
if ~keyword_set(quiet) then box_message, ['Sending B_z component to sp_quick_cube']
bzcomp = sp_quick_cube(data(*,*,wh1[0])*data(*,*,wh4[0])*cos(data(*,*,wh2[0])*!dtor),scan_info,wmap=wmap,/interp,nbad=nbad)
if ~keyword_set(quiet) then box_message, ['Sending B_x component to sp_quick_cube']
bxcomp = sp_quick_cube(data(*,*,wh1[0])*data(*,*,wh4[0])*sin(data(*,*,wh2[0])*!dtor)*cos(data(*,*,wh3[0])*!dtor),scan_info,wmap=wmap,/interp,nbad=nbad)
if ~keyword_set(quiet) then box_message, ['Sending B_y component to sp_quick_cube']
bycomp = sp_quick_cube(data(*,*,wh1[0])*data(*,*,wh4[0])*sin(data(*,*,wh2[0])*!dtor)*sin(data(*,*,wh3[0])*!dtor),scan_info,wmap=wmap,/interp,nbad=nbad)
if ~keyword_set(quiet) then box_message, ['Sending fill fac array to sp_quick_cube']
fillfac = sp_quick_cube(data(*,*,wh4[0]),scan_info,wmap=wmap,/interp,nbad=nbad) ;, /add_noise)

;  treat fillfac < 0.02 as non-magnetic
ww = where(fillfac lt 0.02, count)
if (count gt 0) then fillfac(ww) = 1.e7
fillfac_sp_qc = fillfac

;  calculate interpolated components
invfill=1.0/fillfac
magfield = sqrt(bzcomp^2 + bxcomp^2 + bycomp^2)*invfill
magincl = !radeg*atan(sqrt(bxcomp^2+bycomp^2)*invfill,bzcomp*invfill)
magazim = !radeg*atan(bycomp*invfill,bxcomp*invfill)
ww = where(magazim lt 0.,count)
if (count gt 0) then magazim(ww) = magazim(ww)+180.
magfield_sp_qc = magfield
magincl_sp_qc  = magincl
magazim_sp_qc  = magazim

;
;5.) magfield
;; if ~keyword_set(quiet) then box_message, ['Sending magfield array to sp_quick_cube']
;; wh=where(index2.ftype eq 'Magnetic Field Strength',nwh)  ;  usually 0
;; if nwh eq 0 then begin
;;   print,'  quantity not found in level 2 data'
;;   return
;; endif
;; magfield = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,interp=0,nbad=nbad,loud=loud)
;; magfield_sp_qc = magfield
;
;6.) magincl
;; if ~keyword_set(quiet) then box_message, ['Sending magincl array to sp_quick_cube']
;; wh=where(index2.ftype eq 'Magnetic Field Inclination',nwh)  ;  usually 1
;; if nwh eq 0 then begin
;;   print,'  quantity not found in level 2 data'
;;   return
;; endif
;; magincl = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud)
;; magincl_sp_qc = magincl
;
;7.) magazim
;; if ~keyword_set(quiet) then box_message, ['Sending magazim array to sp_quick_cube']
;; wh=where(index2.ftype eq 'Magnetic Field Azimuth',nwh)  ;  usually 2
;; if nwh eq 0 then begin
;;   print,'  quantity not found in level 2 data'
;;   return
;; endif
;; magazim = sp_quick_cube(data(*,*,wh[0]),scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud)
;; magazim_sp_qc = magazim
;
;8.) slat - this is an intermediate product and not written to file
;original: slat = sp_quick_cube(slat,scan_info,wmap=wmap,/interp,/xytinterp,loud=loud)
;; slat = sp_quick_cube(slat,scan_info,wmap=wmap,/xytinterp,loud=loud)  ;  removed /interp.  Verify with Marc that this is correct -- reference the copy of Ted's e-mail.
;; slat_sp_qc = slat
;
;9.) slong - this is an intermediate product and not written to file
;orig: slong = sp_quick_cube(slong,scan_info,wmap=wmap,/interp,/xytinterp,loud=loud)
;; slong = sp_quick_cube(slong,scan_info,wmap=wmap,/xytinterp,loud=loud)  ;  removed /interp.  Verify with Marc that this is correct -- reference the copy of Ted's e-mail.
;; slong_sp_qc = slong
;
;10.) utimed - this is an intermediate product and not written to file
;original utimed = sp_quick_cube(utimed,scan_info,wmap=wmap,/interp,/xytinterp,loud=loud)
;second original: utimed = sp_quick_cube(utimed,scan_info,wmap=wmap,/interp,/xytinterp,loud=loud); TODO: is there supposed to be an interp here???????????
;; utimed = sp_quick_cube(utimed,scan_info,wmap=wmap,/xytinterp,loud=loud); removed /interp.  Verify with Marc that this is correct -- reference the copy of Ted's e-mail.
;; utimed_sp_qc = utimed
;
;12.) rufaz
;; if ~keyword_set(quiet) then box_message, ['Sending rufaz array to sp_quick_cube']
;; rufaz = sp_quick_cube(rufaz,scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud)
;; rufaz_sp_qc = rufaz 
;; rfix_index = where(btapp lt 125)
;; if n_elements(rfix_index) gt 1 then rufaz(rfix_index) = 0 ;This is to catch the case when a single -1 == no elements is returned.  
;;                                                           ;Ambiguous, but if there is only one bad one, we should be ok.
;
; 1.) blapp
if ~keyword_set(quiet) then box_message, ['Sending blapp array to sp_quick_cube']
blapp = sp_quick_cube(blapp,scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud)
if isa(blapp,/int) then begin
  print, 'blapp is an int and not a struct.  returning.'
  return ; it is supposed to be a cube.  presumably we only need to test once, and we can test at the first one.
endif
blapp_sp_qc = blapp

;
;2.) btapp
if ~keyword_set(quiet) then box_message, ['Sending btapp array to sp_quick_cube']
btapp = sp_quick_cube(btapp,scan_info,wmap=wmap,/interp,nbad=nbad,loud=loud)
btapp_sp_qc = btapp

; Now convert to magnetic field components and interpolate those, per
; KD Leka's recommendation
; avoids mixing angles like 0 & 180 degrees which are equivalent in magazim

;; fillfac = data(*,*,11)

;; if ~keyword_set(quiet) then box_message, ['Sending B_z component to sp_quick_cube']
;; bzcomp = sp_quick_cube(data(*,*,0)*cos(data(*,*,1)*!dtor)*fillfac,scan_info,wmap=wmap,/interp,nbad=nbad)

;; if ~keyword_set(quiet) then box_message, ['Sending B_x component to sp_quick_cube']
;; bxcomp = sp_quick_cube(data(*,*,0)*sin(data(*,*,1)*!dtor)*fillfac*cos(data(*,*,2)*!dtor), $
;;                        scan_info,wmap=wmap,/interp,nbad=nbad)

;; if ~keyword_set(quiet) then box_message, ['Sending B_y component to sp_quick_cube']
;; bycomp = sp_quick_cube(data(*,*,0)*sin(data(*,*,1)*!dtor)*fillfac*sin(data(*,*,2)*!dtor), $
;;                        scan_info,wmap=wmap,/interp,nbad=nbad)

;; ;14) fillfac
;; if ~keyword_set(quiet) then box_message, ['Sending fill fac array to sp_quick_cube']
;; fillfac = sp_quick_cube(data(*,*,11),scan_info,wmap=wmap,/interp,nbad=nbad) ;, /add_noise)
;; fillfac_sp_qc = fillfac

;; ; Recover magfield, etc. from interpolated components
;; fill2 = fillfac
;; ; treat fillfac < 0.02 (TBC) as non-magnetic
;; ww = where(fill2 lt 0.02, count)
;; if (count gt 0) then fill2(ww) = 1.e7
;; fill2 = 1./fill2
;; magfield = sqrt(bzcomp^2 + bxcomp^2 + bycomp^2)*fill2
;; magincl = !radeg*atan(sqrt(bxcomp^2+bycomp^2)*fill2,bzcomp*fill2)

;; w_magincl_good = where(wmap gt 0) ;exclude all initially bad one, even if interpolated 
;; w_magincl_bad  = where(wmap lt 0)

;; magazim = !radeg*atan(bycomp*fill2,bxcomp*fill2)

;; ww = where(magazim lt 0.,count)
;; if (count gt 0) then magazim(ww) = magazim(ww)+180.
;; print,'  mag field change after interpolation (stdev) ',stdev(magfield-magfield_sp_qc)
;; print,'inclination change after interpolation (stdev) ',stdev(magincl-magincl_sp_qc)
;; print,'    azimuth change after interpolation (stdev) ',stdev(magazim-magazim_sp_qc)
;; ;help,magfield
;; ww = where(wmap eq -3, count3)
;; ww = where(wmap eq -2, count2)
;; ww = where(wmap eq -1, count1)

;; print,'Missing columns wmap = -3: ',count3
;; print,'Missing columns wmap = -2: ',count2
;; print,'Missing columns wmap = -1: ',count1

;; magfield_sp_qc = magfield
;; magincl_sp_qc  = magincl
;; magazim_sp_qc  = magazim

outdir = string([outdir+strtrim('/',2)])

;----------------------------- End of quick cube generation -------------------------;

;Dimensional information is used in building the output path, as well as in image consruction.
height_image = (size(dopp1))(2)
width_image  = (size(dopp1))(1)

;test only start
minimum_allowable_width =20.
minimum_allowable_height = 128.
if width_image le minimum_allowable_width then begin
  print, 'minimum allowable image width = ', strtrim(string(minimum_allowable_width),2)
  print, 'original image width = ', strtrim(string(width_image),2), ' original image height = ', strtrim(string(height_image),2)
  new_width_factor = minimum_allowable_width/width_image
  width_image = fix(new_width_factor*width_image) + fix(new_width_factor*width_image mod 2)
  if ( fix(new_width_factor*height_image) + fix(new_width_factor*height_image mod 2) le 8192 ) and (width_image gt 20) then begin ; width_image gt 20 <<-- kind of a hack.
    height_image= fix(new_width_factor*height_image) + fix(new_width_factor*height_image mod 2)
  endif 
  print, 'new image width = ', strtrim(string(width_image),2), ' new image height = ', strtrim(string(height_image),2)
endif

num_dims   = size(dopp1,/n_dimensions)
old_style_image_markup = 0; currently default, make this keyword_driven.
new_style_image_markup = 1
if num_dims eq 3 then num_images =  (size(dopp1))(3) else num_images = 1

;Building the base, movie, image, and fits directories. -- START 
if num_dims le 2 then begin
  if strmid(outdir,0,1) ne "/" then begin 
    if keyword_set (flat_dir_tree) then begin
      outdir = curdir() + "/" + dstring + "/"
    endif else begin 
      outdir = curdir() + "/2D/" + anytim2dir(dstring) + "/" + strtrim(dstring,2) +"/"
    endelse
    single_img_dir = outdir +"image/"
  endif else begin
    if keyword_set (flat_dir_tree) then begin
      outdir = outdir + "/" + strtrim(dstring,2) +"/"
    endif else begin
      outdir = outdir + "/2D/"+"/"+anytim2dir(dstring)+"/"+dstring + "/"
    endelse
    single_img_dir = outdir +"image/"
  endelse
  mk_dir, outdir
  if ~keyword_set(fits_only) then mk_dir, single_img_dir
endif else begin ;should be the case where num_dims eq 3
  if strmid(outdir,0,1) ne "/" then begin
    if keyword_set(flat_dir_tree) then begin
      outdir = outdir+ "/" + strtrim(dstring,2) + "/" 
    endif else begin
      outdir =  outdir + '/3D/' + anytim2dir(dstring) +"/"+ strtrim(dstring,2)+"/" ;concatanating dstring is new -- need to also be done in the 2d case TODO
    endelse 
    movie_dir = outdir + "movie"
  endif else begin
    if keyword_set(flat_dir_tree) then begin
      outdir = outdir + "/" + strtrim(dstring,2)+"/"
    endif else begin
      outdir =  outdir +'/3D/' + anytim2dir(dstring) + "/"+strtrim(dstring,2)+"/" ;
    endelse
    movie_dir = outdir+"movie/"
  endelse
  mk_dir, outdir
  if ~keyword_set(fits_only) then mk_dir, movie_dir
endelse
outdir_fits = outdir[0] + '/fits/' ;will it work from here?
mk_dir, outdir_fits ;always
savegenx, wmap, file= outdir_fits+'wmap_'+dstring,/overwrite
;------------------ Start of FITS writing Section -------------------;

;Build Index Structure for Data Cubes -- START
csize = size(blapp_sp_qc) ;all "_qc" data cubes are presumed to be the same size. So, we can query an _qc data product to determine csize.
nx = csize(1) 
ny = csize(2)
;; if (csize(1) ne nx) or (csize(2) ne ny) then begin  ;  MLD notes that eventually we want to use output from sotsp_rasteriize to set nx and ny
;;   print,'  dimension mismatch with slitmap vs quick cubes'
;;   stop
;;   return
;; endif

;Using the blapp_sp_qc for num images is correct because we have to make a fits header for each fits image
if n_dims(blapp_sp_qc) eq 3 then begin 
  ;; nt = csize(3)
  if (nt - num_images) eq 1 then nt = num_images
endif else begin
  if n_dims(blapp_sp_qc) eq 2 then nt = 1 ;using blapp here is wrong, because it might have had its last image removed (if more than 50 percent of its columns were uninterpolated) and these are for the fits files.
endelse

cubeindex=0

;note: for compatibility with iris raster data readers, the following keywords are
;a.) Required:     NAXIS, NAXISn, CDELTn, CRPIXn, CRVALn, CROTAn or SAT-ROT or pc-matrix,
;                  DATE_OBS
;b.) Very useful:  CTYPEn, CUNITn, BTYPE, BUNIT, EXPTIME
;c.) nice to have: TELESCOP, INSTRUME, LVL_NUM, ORIGIN, BLD_VERS, CAMERA, DSUN_OBS, 
;                  WAVE_STR, WAVELNTH, KEYWDDOC, HISTORY
;Missing tags from lists a, b, c are:
;Missing "a": None. Note: using CROTAn
;Missing "b": None. 
;Missing "c": ORIGIN, BLD_VERS, CAMERA, DSUN_OBS, note: we have 'wave', but not 'wave_str' or 'wavelnth', 
;             KEYWDDOC, HISTORY
;
;add some information to the primary header
taglist=[          'simple',   'bitpix',  'naxis',  'naxis1',   'naxis2',   'telescop', 'instrume' ] 
taglist=[taglist , 'date_obs', 'cdelt1',  'cdelt2', 'crpix1',   'crpix2',   'xcen',     'ycen'     ]
taglist=[taglist,  'crval1',   'crval2',  'crota1', 'crota2',   'ctype1',   'ctype2',   'cunit1'   ]
taglist=[taglist,  'cunit2',   'exptime', 'wave',   'prog_ver', 'seqn_ver', 'parm_ver', 'nslitpos' ]
taglist=[taglist,  'scn_sum',  'scn_rpt', 'spmapctr','scn_step', 'spccdix0', 'spccdix1', 'num_side' ]
taglist=[taglist,  'camssum', 'sp_extid', 'macroid' ]
;
; index1 is a level 1 index
;; sotsp_stks2index,dstring,index1;  ;  MLD note: isn't index1 determined above?
if (index1[0].nslitpos eq 0) && (index1.scn_sum gt 1) then stop, 'Odd value for nslitpos, verify'
if index1[0].scn_sum eq 2 then index1.nslitpos = index1.nslitpos/2 
;also, what do we do if nslitpos is an odd number?  Is that a possibility? 

;will sot_fix_pointing be integrated into this code? (or, can this be removed?)
;index1_temp = index1
;index1 = sot_fix_pointing(index1, oindex=oindex1,/sp,/verbose) ;doesn't always have to be verbose.

is_index1_valid = isa(index1)
if 0 eq is_index1_valid then begin 
  if ~keyword_set(quiet) then print, 'no valid return from sotsp_stks2index so index1 is undefined -- Returning'
  return 
endif else if ~keyword_set(quiet) then print, 'Valid index1 -- Continuing'

newindex=str_subset(index1(0),taglist,/include)

;From the level 2 header
ephem=pb0r(index2(0).tstart)
newindex=add_tag(newindex, index2(0).tstart,'tstart' )
newindex=add_tag(newindex, index2(0).tend  ,'tend'   )
newindex=add_tag(newindex, index2(0).mapid ,'mapid'  )
newindex=add_tag(newindex, 'B_LOS'         ,'btype'  )
newindex=add_tag(newindex, 'Gauss'         ,'bunit'  )
newindex=add_tag(newindex, 1.5             ,'lvl_num')
if tag_exist(index2(0),'p_angle') then newindex=add_tag(newindex, index2(0).p_angle ,'p_angle') else newindex=add_tag(newindex, ephem(0) ,'p_angle')
if tag_exist(index2(0),'b_angle') then newindex=add_tag(newindex, index2(0).b_angle ,'b_angle') else newindex=add_tag(newindex, ephem(1) ,'b_angle')
if tag_exist(index2(0),'solar_ra') then newindex=add_tag(newindex, index2(0).solar_ra ,'solar_ra') else newindex=add_tag(newindex, 60*ephem(2) ,'solar_ra')

;
;Build Index Structure for Data Cubes -- END
nx = index1[0].nslitpos

;Populate Data Cube Index Common Values - START
newindex.naxis    = 3
newindex.naxis1   = nx
newindex.naxis2   = ny
newindex.instrume = 'SOT-SP'
newindex.ctype1   = 'Solar-X'
newindex.ctype2   = 'Solar-Y'
;newindex.cdelt1  = index2(0).cdelt1
newindex.cdelt1   = index1(0).scn_step*0.148569*index1(0).scn_sum
;newindex.cdelt2  = index2(0).cdelt2
newindex.cdelt2   = index1(0).camssum*0.1599
newindex.cunit1   = 'arcseconds'
newindex.cunit2   = 'arcseconds'
newindex.scn_sum  = index1(0).scn_sum
newindex.scn_step = index1(0).scn_step
newindex.scn_rpt  = index1(0).scn_rpt
newindex.spmapctr = index1(0).spmapctr
newindex.camssum  = index1(0).camssum
newindex.num_side = index1(0).num_side
newindex.spccdix0 = index1(0).spccdix0
newindex.spccdix1 = index1(0).spccdix1
    
if (newindex.scn_sum gt 1) then newindex.exptime = 2.*newindex.exptime

;Make a cube of index values -- one per image
cubeindex_base = replicate(newindex,nt)

;  determine FITS header keywords for each individual raster

;  shimizu-san's offset table gives bad values after 2015
;; offset = get_shimizu(index2(0).date_obs)
;; if ~keyword_set(quiet) then print,'Shimizu offset ',offset,' for ',index2(0).date_obs

for i=0,nt-1 do begin

  ;  determine DATE_OBS for this single raster (assume linear across scan)
  slitix=where(slitmap(*,i) ge 0,nindices)
  if nindices eq 0 then begin
    box_message,'Issue in sp_make_cube_index.pro with '+dstring+': oddly empty scan?'
    outdir=outdir_orig
    return
  endif
  if nx eq 1 then cubeindex_base(i).date_obs = anytim(scan_info.times(slitmap(slitix,i)),/vms) else begin
    fitresult=linfit(slitix,anytim(scan_info.times(slitmap(slitix,i))))
    cubeindex_base(i).date_obs = anytim(fitresult[0]+fitresult[1]*(0.5*(nx-1)),/vms)
  endelse

  ;  determine XCEN for this single raster (assume linear across scan)
  wh=where(index2.ftype eq 'Map X Coordinates',nwh)  ;  usually 36
  if nwh eq 0 then begin
    print,'  quantity not found in level 2 data'
    return
  endif
  if nx eq 1 then cubeindex_base(i).xcen = data(slitmap(slitix,i),0,wh[0]) else begin
    fitresult=linfit(slitix,data(slitmap(slitix,i),ny/2,wh[0]))   ;  integer division
    cubeindex_base(i).xcen = fitresult[0]+fitresult[1]*(0.5*(nx-1)) ;+offset[0]
  endelse

  ;  determine XCEN for this single raster (take median value from all valid slit positions)
  wh=where(index2.ftype eq 'Map Y Coordinates',nwh)  ;  usually 37
  if nwh eq 0 then begin
    print,'  quantity not found in level 2 data'
    return
  endif
  cubeindex_base(i).ycen = median(data(slitmap(slitix,i),ny/2,wh[0]));+offset[1]  ;  use median Y value to eliminate outliers

  ;  set CRVALi = XCEN or YCEN
  cubeindex_base(i).crval1 = cubeindex_base(i).xcen
  cubeindex_base(i).crval2 = cubeindex_base(i).ycen

  ;  set CRPIXi = central pixel of this single raster
  cubeindex_base(i).crpix1 = (FLOAT(cubeindex_base(i).NAXIS1)+1.)/2.
  cubeindex_base(i).crpix2 = (FLOAT(cubeindex_base(i).NAXIS2)+1.)/2.

  ;  set CROTAi
  cubeindex_base(i).crota1 = 0.41197900                                    ; per Ted (should be this now, but am forcing it to correspond with Ted's definition.)
  cubeindex_base(i).crota2 =  0.41197900+0.26                              ;per Ted

endfor  

;  define output file
if ~keyword_set(outdir) then outdir = './' else outdir = outdir
if ~keyword_set(filename) then filename = 'SPcubes_'+dstring

;  make minimal FITS header
mkhdr, mainheader, 0, /extend ;make initial main header.
cubeindex_main = cubeindex_base
cubeindex_main = rem_tag(cubeindex_main, 'SIMPLE')
cubeindex_main = rem_tag(cubeindex_main, 'BITPIX')
cubeindex_main = rem_tag(cubeindex_main, 'NAXIS')
cubeindex_main = rem_tag(cubeindex_main, 'NAXIS1')
cubeindex_main = rem_tag(cubeindex_main, 'NAXIS2')
cubeindex_main = rem_tag(cubeindex_main, 'DATE')
cubeindex_main = rem_tag(cubeindex_main, 'END')

cube_index_single = cubeindex_main[0]
tags = tag_names(cube_index_single)
n_tags = n_elements(tags)
for i1=0, n_tags-1 do begin
  output=execute(string('tag_val = cubeindex_main[', strtrim(string(0),2),  '].',tags(i1)))
  sxaddpar, mainheader, tags[i1], tag_val
endfor
sxdelpar, mainheader, 'BTYPE'
sxdelpar, mainheader, 'BUNIT'
sxaddhist,'sp_make_cube_index VERSION:'+version_date,mainheader

file_extension = '.fits'
;
;1.) blapp
blapp_str = '_blapp_index'
filename_blapp = strtrim(string(outdir_fits+ filename+ blapp_str+ file_extension),2)
writefits, filename_blapp,  0, mainheader ;top level fits header for 'blapp' fits file .
box_message, ['writing fits files to: ', outdir_fits]
for i_blapp=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_blapp = cubeindex_main[i_blapp] ;this is the index that has had tags removed from it.
  mkhdr, subheader_blapp, blapp[*,*,i_blapp],/image,/extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_blapp, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_blapp.' + string(strtrim(tags(i),2))))
    sxaddpar, subheader_blapp, tags[i], tag_val
  endfor
  writefits, filename_blapp, reform(blapp_sp_qc[*,*,i_blapp]), subheader_blapp, /append ;attempt to write out the correct data array "the sp_qc file is the original file"
endfor
if ~keyword_set(debug) then delvar, blapp, filename_blapp, i_blapp, blapp_str, subheader_blapp, blapp_str, cubeindex_blapp
;
;2.) btapp -- B-Transverse
btapp_str = '_btapp_index'
filename_btapp = strtrim(string(outdir_fits+ filename+ btapp_str+ file_extension),2) 
writefits, filename_btapp,  0, mainheader ;top level fits header for 'btapp' fits file.
for i_btapp=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_btapp       = cubeindex_main[i_btapp] 
  cubeindex_btapp.btype = 'B_transverse'
  mkhdr, subheader_btapp, btapp[*,*,i_btapp], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_btapp, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_btapp.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_btapp, tags[i], tag_val
  endfor
  writefits, filename_btapp, reform(btapp_sp_qc[*,*,i_btapp]), subheader_btapp, /append ; attempt to write out the correct array. the "sp_qc" file is the orignal file
endfor
if ~keyword_set(debug) then delvar, btapp, filename_btapp, i_btapp, subheader_btapp, btapp_str, cubeindex_btapp
;
;3.) conti
conti_str = '_conti_index'
filename_conti = strtrim(string(outdir_fits+ filename+ conti_str+ file_extension),2)
writefits, filename_conti,  0, mainheader ;top level fits header for 'conti' fits file.
for i_conti=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_conti       = cubeindex_main[i_conti]
  cubeindex_conti.btype = 'Continuum Intensity'
  cubeindex_conti.bunit = 'DN'
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  mkhdr, subheader_conti, conti[*,*,i_conti], /image, /extend
  sxaddpar, subheader_conti, 'XTENSION', 'IMAGE'
  for i=0, n_tags - 1 do begin
    output=execute(string('tag_val = cubeindex_conti.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_conti, tags[i], tag_val
  endfor
  writefits, filename_conti, reform(conti_sp_qc[*,*,i_conti]), subheader_conti, /append ;attempt to write out the correct array. The sp_qc file is the original file.
endfor
if ~keyword_set(debug) then delvar,conti, filename_contin, i_conti, subheader_conti, conti_str, cubeindex_conti
;
;4.) dopp2
dopp2_str = '_dopp2_index'
filename_dopp2 = strtrim(string(outdir_fits+ filename+ dopp2_str+ file_extension),2)
writefits, filename_dopp2, 0, mainheader ;top level fits header for 'dopp2' fits file.
for i_dopp2=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_dopp2       = cubeindex_main[i_dopp2]
  cubeindex_dopp2.btype = 'Doppler Shift of 6302.5 Line'
  cubeindex_dopp2.bunit = 'km/s'
  mkhdr, subheader_dopp2, dopp2[*,*,i_dopp2], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_dopp2, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_dopp2.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_dopp2, tags[i], tag_val
  endfor
  writefits, filename_dopp2, reform(dopp2_sp_qc[*,*,i_dopp2]), subheader_dopp2, /append ;attempt to write out the correct array the sp_qc file is the original file.
endfor
if ~keyword_set(debug) then delvar, dopp2, filename_dopp2, i_dopp2, subheader_dopp2, dopp2_str, cubeindex_dopp2
;
;5.) magfield
magfield_str = '_magfield_index'
filename_magfield = strtrim(string(outdir_fits+ filename+ magfield_str+ file_extension),2)
writefits, filename_magfield, 0, mainheader ;top level fits header for 'magfield' fits file.
for i_magfield=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_magfield         = cubeindex_main[i_magfield]
  cubeindex_magfield.btype   = 'Magnetic Field Strength'
  cubeindex_magfield.bunit   = 'Gauss'
  cubeindex_magfield.lvl_num = 2.0
  mkhdr, subheader_magfield, magfield[*,*,i_magfield], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_magfield, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_magfield.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_magfield, tags[i], tag_val 
  endfor 
  writefits, filename_magfield, reform(magfield_sp_qc[*,*,i_magfield]), subheader_magfield, /append ;attempt to write out the correct array the sp_qc file is the original file.
endfor
if ~keyword_set(debug) then delvar, magfield, filename_magfield, i_magfield, subheader_magfield, magfield_str, cubeindex_magfield
;
;6.) magincl
magincl_str = '_magincl_index'
filename_magincl = strtrim(string(outdir_fits+ filename+ magincl_str+ file_extension),2)
writefits, filename_magincl, 0, mainheader ;top level fits header for 'magfield' fits file.
for i_magincl=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_magincl       = cubeindex_main[i_magincl]
  cubeindex_magincl.btype = 'Magnetic Field Inclination'
  cubeindex_magincl.bunit = 'deg'
  mkhdr, subheader_magincl, magincl[*,*,i_magincl], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_magincl, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin 
    output=execute(string('tag_val = cubeindex_magincl.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_magincl, tags[i], tag_val 
  endfor 
  writefits, filename_magincl, reform(magincl_sp_qc[*,*,i_magincl]), subheader_magincl, /append ;attempt to write out the correct array the sp_qc file is the original file.
endfor
if ~keyword_set(debug) then delvar, magincl, filename_magincl, i_magincl, subheader_magincl, magincl_str, cubeindex_magincl
;
;7.) magazim
magazim_str = '_magazim_index'
filename_magazim = strtrim(string(outdir_fits+ filename+ magazim_str+ file_extension),2)
writefits, filename_magazim, 0, mainheader 
for i_magazim=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_magazim       = cubeindex_main[i_magazim]
  cubeindex_magazim.btype = 'Magnetic Field Azimuth'
  cubeindex_magazim.bunit = 'deg'
  mkhdr, subheader_magazim, magazim[*,*,i_magazim], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_magazim, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_magazim.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_magazim, tags[i], tag_val
  endfor 
  writefits, filename_magazim, reform(magazim_sp_qc[*,*,i_magazim]), subheader_magazim, /append ;attempt to write out the correct array the sp_qc file is the original file
endfor
if ~keyword_set(debug) then delvar, magazim, filename_magazim, i_magazim, subheader_magazim, magazim_str, cubeindex_magazim
;
;8.) dopp1
dopp1_str = '_dopp1_index'
filename_dopp1 = strtrim(string(outdir_fits+ filename+ dopp1_str+ file_extension),2)
writefits, filename_dopp1, 0, mainheader
for i_dopp1=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_dopp1       = cubeindex_main[i_dopp1]
  cubeindex_dopp1.btype = 'Doppler Shift of 6301.5 Line'
  cubeindex_dopp1.bunit = 'km/s' 
  mkhdr, subheader_dopp1, dopp1[*,*,i_dopp1], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_dopp1, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_dopp1.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_dopp1, tags[i], tag_val
  endfor
  writefits, filename_dopp1, reform(dopp1_sp_qc[*,*,i_dopp1]), subheader_dopp1, /append ; attempt to write out correct data to fits file. the sp_qc file is the original file
endfor
if ~keyword_set(debug) then delvar, dopp1, filename_dopp1, i_dopp1, subheader_dopp1, dopp1_str, cubeindex_dopp1
;
;9.) ptot
ptot_str = '_ptot_index'
filename_ptot = strtrim(string(outdir_fits+ filename+ ptot_str+ file_extension),2)
writefits, filename_ptot, 0, mainheader
for i_ptot=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_ptot       = cubeindex_main[i_ptot]
  cubeindex_ptot.btype = 'Total Polarization'
  cubeindex_ptot.bunit = 'unitless number between 0 and 1'
  mkhdr, subheader_ptot, ptot[*,*,i_ptot], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_ptot, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_ptot.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_ptot, tags[i], tag_val
  endfor
  writefits, filename_ptot, reform(ptot_sp_qc[*,*,i_ptot]), subheader_ptot, /append ;attempt to write out correct data to fits file. the sp_qc file is the original file
endfor
if ~keyword_set(debug) then delvar, ptot, filename_ptot, i_ptot, subheader_ptot, ptot_str, cubeindex_ptot
;
;10.) rufaz
;; rufaz_str = '_rufaz_index'
;; filename_rufaz = strtrim(string(outdir_fits+ filename+ rufaz_str+ file_extension),2)
;; writefits, filename_rufaz, 0, mainheader
;; for i_rufaz=0, n_elements(cubeindex_base)-1 do begin
;;   cubeindex_rufaz       = cubeindex_main[i_rufaz]
;;   cubeindex_rufaz.btype = 'Rough Magnetic Field Azimuth'
;;   cubeindex_rufaz.bunit = 'Degrees counterclockwise from Solar West'
;;   mkhdr, subheader_rufaz, rufaz[*,*,i_rufaz], /image, /extend
;;   tags = tag_names(cubeindex_main)
;;   n_tags = n_elements(tags)
;;   sxaddpar, subheader_rufaz, 'XTENSION', 'IMAGE'
;;   for i=0, n_tags-1 do begin
;;     output=execute(string('tag_val = cubeindex_rufaz.'+ string(strtrim(tags(i),2))))
;;     sxaddpar, subheader_rufaz, tags[i], tag_val
;;   endfor
;;   writefits, filename_rufaz, reform(rufaz_sp_qc[*,*,i_rufaz]), subheader_rufaz, /append ; the sp_qc file is the original file
;; endfor
;; if ~keyword_set(debug) then delvar, rufaz, filename_rufaz, i_rufaz, subheader_rufaz, rufaz_str, cubeindex_rufaz
;
;11.) fillfac
fillfac_str = '_fillfac_index'
filename_fillfac = strtrim(string(outdir_fits+ filename+ fillfac_str+ file_extension),2)
writefits, filename_fillfac, 0, mainheader
for i_fillfac=0, n_elements(cubeindex_base)-1 do begin
  cubeindex_fillfac       = cubeindex_main[i_fillfac]
  cubeindex_fillfac.btype = 'Magnetic Filling Factor'
  cubeindex_fillfac.bunit = 'unitless number between 0 and 1'
  mkhdr, subheader_fillfac, fillfac[*,*,i_fillfac], /image, /extend
  tags = tag_names(cubeindex_main)
  n_tags = n_elements(tags)
  sxaddpar, subheader_fillfac, 'XTENSION', 'IMAGE'
  for i=0, n_tags-1 do begin
    output=execute(string('tag_val = cubeindex_fillfac.'+ string(strtrim(tags(i),2))))
    sxaddpar, subheader_fillfac, tags[i], tag_val
  endfor
  writefits, filename_fillfac, reform(fillfac_sp_qc[*,*,i_fillfac]), subheader_fillfac, /append ; the sp_qc file is the original file
endfor
if ~keyword_set(debug) then delvar, fillfac, filename_fillfac, i_fillfac, subheader_fillfac, fillfac_str, cubeindex_fillfac
if ~keyword_set(debug) then delvar, mainheader, cubeindex_main

cd, startdir
outdir = outdir_orig
;------------------ End of FITS writing -------------------;
;
;------------------ Start of image/movie code -------------;

if ~keyword_set(fits_only) then begin

  modified_width  = 0
  modified_height = 0
  if 0 ne (width_image mod 2) then begin
    width_image = width_image + 1
    modified_width = 1
  endif
  if 0 ne (height_image mod 2) then begin
    height_image = height_image +1
    modified_height = 1
  endif
  if (0 ne modified_width mod 2) || (0 ne modified_height mod 2) then begin
    blapp    = congrid(blapp,    width_image, height_image, num_images,/interp)
    btapp    = congrid(btapp,    width_image, height_image, num_images,/interp)
    conti    = congrid(conti,    width_image, height_image, num_images,/interp)
    dopp1    = congrid(dopp1,    width_image, height_image, num_images,/interp)
    magfield = congrid(magfield, width_image, height_image, num_images,/interp)
    magincl  = congrid(magincl,  width_image, height_image, num_images,/interp)
    magazim    = congrid(magazim,    width_image, height_image, num_images,/interp)
  endif

  ;image out size adjustment -- verify that these rules are ok., i.e., that they don't puff-up the images too much.  Also, does this add appreciable time to movie generation?
  if width_image le 60 then begin
    min_width_image = 64.
    max_width_image = floor(8190/7.) ;= 1170 (clamped at 8192) verify that width_image is a pre-concatentated image.  That is the reason we are diving by 7
    min_height_image = 128.
    max_height_image = 1024 ; <-- just a guess 
    new_dim_factor_x = 1.0
    new_dim_factor_y = 1.0
    new_dim_factor = 1

    if width_image  gt max_width_image then begin 
      print, 'image width ge max_width_image, Returning.'
      return
    endif
    if height_image gt max_height_image then begin
      print, 'image height ge max_height_image.  Returning.'
      return
    endif

    if ~keyword_set(quiet) then box_message, ['Old x,y dimensions were: ', string(width_image), string(height_image)]
    ;here we build candidate dim_factors in the case that a dimension is too 'small'.
    if width_image  lt min_width_image  then new_dim_factor_x    = max([ceil((min_width_image)/width_image)  , 1.])     ;this new_dim_factor_x will yield the maximum x possible
    if height_image lt min_height_image then new_dim_factor_y    = max([ceil((min_height_image)/height_image), 1.])     ;this new_dim_factor_y will yield the maximum y possible

    if ((new_dim_factor_x*width_image le max_width_image) and (new_dim_factor_y*height_image le max_height_image)) $
      then begin ;now we want to find the maximum multiplier to get final x, y dimensions that preserve the aspect ratio
      dim_factor_x_test = new_dim_factor_x
      dim_factor_y_test = new_dim_factor_y
      for ii = floor(min([new_dim_factor_x, new_dim_factor_y])), max([(floor(max_width_image/width_image)), floor(max(max_height_image/height_image))])-1 do begin
        if (((dim_factor_x_test+1.)*width_image le max_width_image) and ((dim_factor_y_test+1.)*height_image le max_height_image)) then begin
          dim_factor_x_test = dim_factor_x_test + 1
          dim_factor_y_test = dim_factor_y_test + 1
        endif
      endfor
      new_dim_factor_x = dim_factor_x_test       ; these might remain the same
      new_dim_factor_y = dim_factor_y_test       ; these might remain the same
    endif
    new_dim_factor = min([new_dim_factor_x, new_dim_factor_y])
    
    if 0 ne (new_dim_factor_x*width_image)  mod 2 then width_image  = width_image -1     ;make divisible by 2
    if 0 ne (new_dim_factor_y*height_image) mod 2 then height_image = height_image-1     ;make divisible by 2
    new_dim_factor_x = new_dim_factor
    new_dim_factor_y = new_dim_factor
    if width_image lt 50 and new_dim_factor_x eq 1.0 then new_dim_factor_x = 2. ;tryint it our for case: '20160517_121728' pgs 20180403 
    if 10./(new_dim_factor_x*width_image) gt 1 then new_dim_factor_x = ceil(min_width_image/width_image)
    if ~keyword_set(quiet) then box_message, ['Multiplying individual (x,y) images by factors of ' + strtrim(string(new_dim_factor_x),2) + ' ' +strtrim(string(new_dim_factor_y,2))]
    width_image  = new_dim_factor_x*width_image
    height_image = new_dim_factor_y*height_image
    if ~keyword_set(quiet) then box_message, ['New x,y dimensions and dim_factor are: ', string(width_image), string(height_image), string(new_dim_factor)]
    blapp    = congrid(blapp,    width_image, height_image, num_images,/interp)
    btapp    = congrid(btapp,    width_image, height_image, num_images,/interp)
    conti    = congrid(conti,    width_image, height_image, num_images,/interp)
    dopp1    = congrid(dopp1,    width_image, height_image, num_images,/interp)
    magfield = congrid(magfield, width_image, height_image, num_images,/interp)
    magincl  = congrid(magincl,  width_image, height_image, num_images,/interp)
    magazim  = congrid(magazim,  width_image, height_image, num_images,/interp)
  endif

  off_limb_mask=conti*0. + 1.
  obs_xy_radial_dist_from_center =  sqrt(index2[0].xcen^2.0 + index2[0].ycen^2.0) ; NEED TO RECOVER OR CALCULATE THIS TODO
  solar_radius = newindex.solar_ra                                                ;
  fract_radial_dist = obs_xy_radial_dist_from_center/solar_radius
  if ~keyword_set(quiet) then box_message, ['fractional distance of obs to limb is: ', strtrim(string(fract_radial_dist),2)]
  if fract_radial_dist gt .90 then begin ;this is just to the center of the FOV.  If we know the FOV, which we should, then we should be able to calculate if there is any limbing in the image or not
                                         ;and replace this with an estimate
    LIMB_THRESHOLD_DISCRIMINATOR = 0.35
    where_dn_low = where(conti le LIMB_THRESHOLD_DISCRIMINATOR*median(conti)) ;
    off_limb_mask[where_dn_low] = 0.
    if ~keyword_set(quiet) then box_message, ['Near limb so using bit_mask to set all values <= 0.35*median(conti) to black']
  endif

  ;------------------------------------- 2D Case start -------------------------------;
  if num_dims le 2 then begin
    existing_2d_image = file_exist(single_img_dir+'/'+dstring+'.gif')
    single_sp_image = 1b 
    agg_img        = bytarr(width_image*7, height_image)
    agg_img_fits_t = fltarr(width_image*7, height_image) ;new
    if not (max(blapp) eq 0 and min(blapp) eq 0) then begin
      blapp_t        = off_limb_mask*bytscl(magscl(off_limb_mask*blapp, 25, 1000))
      if max(blapp_t) ne 0 then  blapp = blapp_t else blapp = bytscl(blapp*off_limb_mask)
    endif
    if not (max(btapp) eq 0 and min(btapp) eq 0) then begin
      btapp_t        = off_limb_mask*bytscl(hist_equal(btapp,    percent=1))
      if max(btapp_t)  then  btapp = btapp_t else btapp = bytscl(btapp*off_limb_mask)
    endif
    if not( max(conti) eq 0 and min(conti) eq 0) then begin
      conti_t        = off_limb_mask*bytscl(hist_equal(conti,    percent=2))
      if max(conti_t) ne 0 then  conti = conti_t else cont = bytscl(conti*off_limb_mask)
    endif
    if not (max(dopp1) eq 0 and min(dopp1) eq 0) then begin
      dopp1_t        = off_limb_mask*bytscl(   255.-hist_equal(dopp1*off_limb_mask,    percent=1)*1  ) ;this reverses 'dopp1' -- need to do this for 2d as well.
      if max(dopp1_t) ne 0 ne 0 then  dopp1 = dopp1_t else dopp1 = off_limb_mask*(255-bytscl(dopp1)) 
    endif
    if not (max(magfield) eq 0 and min(magfield) eq 0) then begin
      magfield_t     = off_limb_mask*bytscl(hist_equal(magfield*off_limb_mask, percent=.3))
      if max(magfield) ne 0 then  magfield = magfield_t else magfield = bytscl(magfield*off_limb_mask)
    endif
    if not (max(magincl) eq 0 and min (magincl) eq 0) then begin
      magincl_t      = off_limb_mask*bytscl(hist_equal(magincl,  percent=1))
      if max(magincl_t) ne 0  then  magincl = magincl_t else magincl = bytscl(magincl*off_limb_mask)
    endif
    if not (max(magazim) eq 0 and min(magazim) eq 0) then begin
      magazim_t        = off_limb_mask*bytscl(hist_equal(magazim,    percent=1))
      if max(magazim_t) ne 0 then  magazim = magazim_t else magazim =  bytscl(off_limb_mask*magazim)
    endif
    sot_fov = 0            ; find a proc to determine the fov
    pointing_check = 1     ; find or write a proc that estimate/determines whether the field of view includes the solar limb.

    if (old_style_image_markup eq 1) then begin
      ;-------add image type labels before aggregating them.  START
      cor = 0
      blapp    = align_label_cube('B_lon', blapp,    /lc )
      btapp    = align_label_cube('B_tv',  btapp,    /lc )
      conti    = align_label_cube('Ic',    conti,    /lc )
      dopp1    = align_label_cube('Vel',   dopp1,    /lc )
      magfield = align_label_cube('Mag',   magfield, /lc )
      magincl  = align_label_cube('Incl',  magincl,  /lc )
      magazim  = align_label_cube('Az',    magazim,    /lc )
      ;-------add image type labels before aggregating them.  END
    endif

    if (new_style_image_markup eq 1) then begin
      if width_image le 60 then tag_fontsize=10 else tag_fontsize=20
      x_text_fract_offset = .1
      y_text_fract_offset = .1*width_image/height_image
      font_size_multiplier = width_image/128
      
      ;blapp
      blapp_font_size = tag_fontsize
      w_blapp         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_blapp         = image(blapp,/overplot,/clip)
      i_blapp         = text(x_text_fract_offset, y_text_fract_offset,'B$_\parallel$',font_size=blapp_font_size, font_color='yellow')
      blapp_check     = i_blapp.CopyWindow()
      blapp1          = reform(blapp_check[0,*,*], width_image, height_image)
      blapp           =blapp1
      w_blapp.close
      heap_free, w_blapp
      w_blapp = 0
      heap_free, i_blapp

      ;btapp
      btapp_font_size = tag_fontsize
      w_btapp         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_btapp         = image(btapp,/overplot,/clip)
      i_btapp         = text(x_text_fract_offset, y_text_fract_offset,'$B_\perp$',font_size=btapp_font_size, font_color='yellow')
      btapp_check     = i_btapp.CopyWindow()
      btapp1          = reform(btapp_check[0,*,*], width_image, height_image)
      btapp=btapp1    
      w_btapp.close   
      heap_free, w_btapp
      w_btapp = 0
      heap_free, i_btapp

      ;conti
      conti_font_size = tag_fontsize
      w_conti         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_conti         = image(conti,/overplot,/clip)
      i_conti         = text(x_text_fract_offset, y_text_fract_offset,'$I_{\rm c}$',font_size=conti_font_size, font_color='yellow')
      conti_check     = i_conti.CopyWindow()
      conti1          = reform(conti_check[0,*,*], width_image, height_image)
      conti=conti1
      w_conti.close
      heap_free, w_conti
      w_conti = 0
      heap_free, i_conti

      ;dopp1
      dopp1_font_size = tag_fontsize
      w_dopp1         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_dopp1         = image(dopp1,/overplot,/clip)
      i_dopp1         = text(x_text_fract_offset, y_text_fract_offset,'$v_{\rm LOS}$',font_size=dopp1_font_size, font_color='yellow')
      dopp1_check     = i_dopp1.CopyWindow()
      dopp11          = reform(dopp1_check[0,*,*], width_image, height_image)
      dopp1=dopp11
      w_dopp1.close
      heap_free, w_dopp1
      w_dopp1 = 0
      heap_free, i_dopp1

      ;magfield
      magfield_font_size = tag_fontsize
      w_magfield         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_magfield         = image(magfield,/overplot,/clip)
      i_magfield         = text(x_text_fract_offset, y_text_fract_offset,'$B$',font_size=magfield_font_size, font_color='yellow')
      magfield_check     = i_magfield.CopyWindow()
      magfield1          = reform(magfield_check[0,*,*], width_image, height_image)
      magfield=magfield1
      w_magfield.close
      heap_free, w_magfield
      w_magfield = 0
      heap_free, i_magfield

      ;magincl
      magincl_font_size = tag_fontsize
      w_magincl         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_magincl         = image(magincl,/overplot,/clip)
      i_magincl         = text(x_text_fract_offset, y_text_fract_offset,'$\gamma$',font_size=magincl_font_size, font_color='yellow')
      magincl_check     = i_magincl.CopyWindow()
      magincl1          = reform(magincl_check[0,*,*], width_image, height_image)
      magincl=magincl1
      w_magincl.close
      heap_free, w_magincl
      w_magincl = 0
      heap_free, i_magincl

      ;magazim
      magazim_font_size = tag_fontsize
      w_magazim         = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
      i_magazim         = image(magazim,/overplot,/clip)
      i_magazim         = text(x_text_fract_offset, y_text_fract_offset,'$\chi$',font_size=magazim_font_size, font_color='yellow')
      magazim_check     = i_magazim.CopyWindow()
      magazim1          = reform(magazim_check[0,*,*], width_image, height_image)
      magazim=magazim1
      w_magazim.close
      heap_free, w_magazim
      w_magazim = 0
      heap_free, i_magazim
      
    endif

    agg_img=[blapp, btapp, conti, dopp1, magfield, magincl, magazim]


    dims = size(agg_img,/dimensions)
    n_dims = size(agg_img,/n_dimensions)
    if dims[1] lt 256 and n_dims eq 2 then begin
      height_adjustment_factor = ceil(256./dims[1])
      agg_img = congrid(agg_img, dims[0]*height_adjustment_factor, dims[1]*height_adjustment_factor,/interp)
    endif

    ;new test start 2D case
    ;-------------------------------------------- draw clock ----------------------------------- start
    temp_time_array = cubeindex_base.date_obs
    for i=0, (num_images-1) do begin
      drawclock, cubeindex_base[0], agg_img[*,*,i],agg_img_temp, pos=4,/shadow
      agg_img[*,*,i] = agg_img_temp
    endfor
    ;-------------------------------------------- draw clock ----------------------------------- end
    ;
    ;-------------------------------------------- add date string ------------------------------ start
    time_labels = strmid(cubeindex_base[0].date_obs,0,11)
    agg_img = align_label_cube(time_labels, agg_img, /ur)
    ;-------------------------------------------- add date string ------------------------------ end
    ;new test end 2D case

    savegenx, wmap, file= outdir_fits+'/wmap_'+dstring,/overwrite  ;do we want to save this inependently of whether an image or a movie is made?
    single_img_file = single_img_dir+dstring+".gif"
  endif else begin

    ;------------------------------------------ 3D case start ----------------------------------

    single_sp_image = 0b
    blapp    = off_limb_mask*bytscl(magscl(off_limb_mask*blapp, 25, 1000))
    btapp    = off_limb_mask*bytscl(hist_equal(btapp,   percent=1))
    conti    = off_limb_mask*bytscl(hist_equal(conti,   percent=2))
    dopp1    = off_limb_mask*bytscl(   255.-hist_equal(dopp1*off_limb_mask, percent=1 )*1) ;this reverses 'conti' -- verify that this was done for 2d, as well.
    magfield = off_limb_mask*bytscl(hist_equal(magfield*off_limb_mask, percent=.3))
    magincl  = off_limb_mask*bytscl(hist_equal(magincl, percent=1))
    magazim    = off_limb_mask*bytscl(hist_equal(magazim,   percent=1))
 
    if (old_style_image_markup eq 1) then begin
      ;-------add image type labels before aggregating them.  START
      cor = 0
      for agg_img_j = 0, num_images-1 do begin
        blapp[*,*,agg_img_j]    = align_label_cube(' B_lon ', blapp[*,*,agg_img_j],    /lc)
        btapp[*,*,agg_img_j]    = align_label_cube(' B_tv ' , btapp[*,*,agg_img_j],    /lc)
        conti[*,*,agg_img_j]    = align_label_cube(' Ic '   , conti[*,*,agg_img_j],    /lc)
        dopp1[*,*,agg_img_j]    = align_label_cube(' Vel '  , dopp1[*,*,agg_img_j],    /lc)
        magfield[*,*,agg_img_j] = align_label_cube(' Mag '  , magfield[*,*,agg_img_j], /lc)
        magincl[*,*,agg_img_j]  = align_label_cube(' Incl ' , magincl[*,*,agg_img_j],  /lc)
        magazim[*,*,agg_img_j]    = align_label_cube(' Az '   , magazim[*,*,agg_img_j],    /lc)
      endfor
      ;-------add image type labels before aggregating them.  END
    endif
 
    agg_img = bytarr(width_image*7, height_image, num_images)
    print, 'Processing ', strtrim(string(num_images),2),' images: '
    for agg_i=0,num_images-1 do begin
      print, string(13b), (agg_i+1), Format='(A, I, $)'
      if agg_i eq (num_images - 1) then print, ''
      if (new_style_image_markup eq 1) then begin
        if width_image le 60 then tag_fontsize=10 else tag_fontsize=20
        font_size_string1 = 'image width = ' + strtrim(string(width_image),2)
        font_size_string2 = 'font_size in image = ' + strtrim(string(tag_fontsize),2)
        if ~keyword_set(quiet) then box_message = [font_size_string1, font_size_string2]
        x_text_fract_offset = .1
        y_text_fract_offset = .1*width_image/height_image
        font_size_multiplier = width_image/128

        ;
        ;blapp
        blapp_font_size = tag_fontsize
        w_blapp = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_blapp = image(blapp[*,*,agg_i],/overplot,/clip)
        i_blapp = text(x_text_fract_offset, y_text_fract_offset,'B$_\parallel$',font_size=blapp_font_size, font_color='yellow') ; test only
        blapp_check      = i_blapp.CopyWindow()                                                                                 ;test
        blapp1           = reform(blapp_check[0,*,*], width_image, height_image)
        blapp[*,*,agg_i] = blapp1
        w_blapp.close
        heap_free, w_blapp
        w_blapp = 0
        heap_free, i_blapp
        delvar, blapp_check
        delvar, blapp1

        ;btapp
        btapp_font_size = tag_fontsize
        w_btapp = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_btapp = image(btapp[*,*,agg_i],/overplot,/clip)
        i_btapp = text(x_text_fract_offset, y_text_fract_offset,'$B_\perp$',font_size=btapp_font_size, font_color='yellow')
        btapp_check      = i_btapp.CopyWindow()
        btapp1           = reform(btapp_check[0,*,*], width_image, height_image)
        btapp[*,*,agg_i] = btapp1
        w_btapp.close
        heap_free, w_btapp
        w_btapp = 0
        heap_free, i_btapp
        delvar, btapp_check
        delvar, btapp1
 
        ;conti
        conti_font_size = tag_fontsize
        w_conti = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_conti = image(conti[*,*,agg_i],/overplot,/clip)
        i_conti = text(x_text_fract_offset, y_text_fract_offset,'$I_{\rm c}$',font_size=conti_font_size, font_color='yellow')
        conti_check       = i_conti.CopyWindow()
        conti1            = reform(conti_check[0,*,*], width_image, height_image)
        conti[*,*,agg_i] = conti1
        w_conti.close
        heap_free, w_conti
        w_conti = 0
        heap_free, i_conti
        delvar, conti_check
        delvar, conti1
      
        ;dopp1
        dopp1_font_size = tag_fontsize
        w_dopp1 = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_dopp1 = image(dopp1[*,*,agg_i],/overplot,/clip)
        i_dopp1 = text(x_text_fract_offset, y_text_fract_offset,'$v_{\rm LOS}$',font_size=dopp1_font_size, font_color='yellow')
        dopp1_check      = i_dopp1.CopyWindow()
        dopp11           = reform(dopp1_check[0,*,*], width_image, height_image)
        dopp1[*,*,agg_i] = dopp11
        w_dopp1.close
        heap_free, w_dopp1
        w_dopp1 = 0
        heap_free, i_dopp1
        delvar, dopp1_check
        delvar, dopp11
 
        ;magfield
        magfield_font_size = tag_fontsize
        w_magfield = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_magfield = image(magfield[*,*,agg_i],/overplot,/clip)
        i_magfield          = text(x_text_fract_offset, y_text_fract_offset,'$B$',font_size=magfield_font_size, font_color='yellow')
        magfield_check      = i_magfield.CopyWindow()
        magfield1           = reform(magfield_check[0,*,*], width_image, height_image)
        magfield[*,*,agg_i] = magfield1
        w_magfield.close
        heap_free, w_magfield
        w_magfield = 0
        heap_free, i_magfield
        delvar, magfield_check
        delvar, magfield1

        ;magincl
        magincl_font_size = tag_fontsize
        w_magincl = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_magincl = image(magincl[*,*,agg_i],/overplot,/clip)
        i_magincl = text(x_text_fract_offset, y_text_fract_offset,'$\gamma$',font_size=magincl_font_size, font_color='yellow')
        magincl_check      = i_magincl.CopyWindow()
        magincl1           = reform(magincl_check[0,*,*], width_image, height_image)
        magincl[*,*,agg_i] = magincl1
        w_magincl.close
        heap_free, w_magincl
        w_magincl = 0
        heap_free, i_magincl 
        delvar, magincl_check
        delvar, magincl1
 
        ;magazim
        magazim_font_size = tag_fontsize
        w_magazim = window(dimensions=[width_image, height_image],/buffer,/no_toolbar)
        i_magazim = image(magazim[*,*,agg_i],/overplot,/clip)
        i_magazim = text(x_text_fract_offset, y_text_fract_offset,'$\chi$',font_size=magazim_font_size, font_color='yellow')
        magazim_check      = i_magazim.CopyWindow()
        magazim1           = reform(magazim_check[0,*,*], width_image, height_image)
        magazim[*,*,agg_i] = magazim1
        w_magazim.close
        heap_free, w_magazim 
        w_magazim = 0
        heap_free, i_magazim 
        delvar, magazim_check
        delvar, magazim1

      endif 
 
      agg_img[*,*,agg_i] = [blapp[*,*,agg_i]   , $
                            btapp[*,*,agg_i]   , $
                            conti[*,*,agg_i]   , $
                            dopp1[*,*,agg_i]   , $
                            magfield[*,*,agg_i], $
                            magincl[*,*,agg_i] , $
                            magazim[*,*,agg_i]]
 
    endfor ;end of for loop for all of the images
 
    ;Remove last image in a movie if more than 50% is bad.'  
    FRACT_THRESH_REMOVE_LAST_IMAGE = 0.50
    if num_images gt 2 then begin  
      num_elements_wmap_last = n_elements(wmap[*,num_images-1])
      where_last_image_lt_0  =where(wmap[*,(num_images-1)] lt 0 ,count_last_frame_lt_0)
      if (1.*count_last_frame_lt_0)/num_elements_wmap_last ge FRACT_THRESH_REMOVE_LAST_IMAGE then begin
        if ~keyword_set(quiet) then print, 'More than 50 percent of columns in the last frame are interpolated -- removing the last frame from movie'
        num_images = num_images-1
        agg_img = agg_img[*,*,0:(num_images-1)] ;num_images was reduced by 'one' in the prior line, making this work correctly.
      endif
    endif

    ;if we assume that this is always a movie, and therefore has 3 dimensions:
    ;this is new for short images, I don't know if I will have to duplicate this in the 3d case.
    dims = size(agg_img,/dimensions)
    n_dims = size(agg_img,/n_dimensions)
    if dims[1] lt 256 and n_dims eq 3 then begin
      height_adjustment_factor = ceil(256./dims[1])
      agg_img = congrid(agg_img, dims[0]*height_adjustment_factor, dims[1]*height_adjustment_factor, dims[2],/interp)
    endif

    ;-------------------------------------------- draw clock ----------------------------------- start
    temp_time_array = cubeindex_base.date_obs
    for i=0, (num_images-1) do begin
      drawclock, cubeindex_base[i], agg_img[*,*,i],agg_img_temp, pos=4,/shadow
      agg_img[*,*,i] = agg_img_temp
    endfor
    ;-------------------------------------------- draw clock ----------------------------------- end
    ;
    ;-------------------------------------------- add date string ------------------------------ start
    time_labels = strmid(cubeindex_base[0:(num_images-1)].date_obs,0,11)
    agg_img = align_label_cube(time_labels, agg_img, /ur)
    ;-------------------------------------------- add date string ------------------------------ end

    movie_name_prefix = "SP_"
    movie_name_suffix = "__B_lon__B_tv__Ic__Vel__Mag__Incl__Az_"
    movie_name = strjoin([movie_name_prefix,  dstring,  movie_name_suffix],/single)

    htmldoc=movie_dir+movie_name+'1'+'.html' ; Top level html doc name
  endelse  ;movie generation flag set.

  if single_sp_image then ssw_write_gif, single_img_file, bytscl(agg_img) else begin
    if ~keyword_set(quiet) then print, 'width image = '+ strtrim(string(width_image),2) + ' height_image = ' +strtrim(string(height_image),2)
    ;original special_movie2, temp_time_array, agg_img,  movie_dir=movie_dir, movie_name=movie_name,/nolabel;
    movie_dir2 = movie_dir[0]
    special_movie2, temp_time_array, agg_img,  movie_dir=movie_dir2, movie_name=movie_name,/nolabel,/no_gifanimate,/thumbnail ;

      ;; movie_path = movie_dir+movie_name+'_m.mp4'
      ;; ;
      ;; if ~keyword_set(quiet) then box_message, ['debug: Before first movie spawn.']
      ;; ;
      ;; num_mp4_movie_repeats = 1 ; This code is from an earlier design. There are more efficient ways to do this.  No longer needed for ffmpeg.  Am leaving in place.
      ;; openw, lun, movie_dir + '/mylist.txt', /get_lun
      ;; for jj = 0, num_mp4_movie_repeats-1 do begin ;make into switch above
      ;;   movie_concat_line = 'file ' + "'" +movie_dir + movie_name + 'A.mp4'+"'"
      ;;   printf, lun, movie_concat_line
      ;; endfor
      ;; close,lun
      ;; free_lun,lun

      ;; ;FIRST ffmpeg movie spawn
      ;; if n_elements(ffmpeg_path) eq 0 then ffmpeg_path = './ffmpeg'
      ;; if ~file_exist(ffmpeg_path) then begin
      ;;   print, 'no valid ffmpeg_path file associated with ffmpeg_path keyword (or, if not passed in, in the execution directory for sp_make_cube_index.pro.'
      ;;   return
      ;; endif
      ;; wait, 20 ; is a wait really necessary, and, if so, can it be a bit more fine grained?
      ;; box_message, ['Current ffmpeg_path is: '+ffmpeg_path, 'If ffmpeg not found, supply a path to a recent version of ffmpeg.'] 
      ;; x_platform_movie_str = ffmpeg_path + ' -i ' + movie_path + ' -c:v libx264 -profile:v high -level 4.1 -pix_fmt yuv420p -movflags  +faststart ' +  movie_dir+movie_name+'A.mp4'
      ;; spawn, x_platform_movie_str, listing
      ;; wait, 10 ;do we really need a wait here? 
      
      ;; if ~keyword_set(quiet) then print, 'Before second movie spawn.'

      ;; movie_multi_rpt_str = ffmpeg_path + ' -f concat -safe 0 -i '+movie_dir+'mylist.txt' + ' -codec copy ' + movie_dir+movie_name +'_m2.mp4'
      ;; spawn, movie_multi_rpt_str, listing2
      ;; wait, 20 ;do we really need a wait here for the spawn to complete?
      ;; if (file_exist(movie_dir+movie_name+'A.mp4')) then file_delete, movie_dir+movie_name+'A.mp4'
      ;; if (file_exist(movie_dir+movie_name+'_m.mp4')) then file_delete, movie_dir+movie_name+'_m.mp4'
      ;; if (file_exist(movie_dir+movie_name+'_m2.mp4')) then file_move, movie_dir+movie_name+'_m2.mp4', movie_path
  endelse 
endif

end