;+
; NAME:
;   mdrp_reduceoneifu
;
; PURPOSE:
;   Reduce a single IFU through the 3d MaNGA pipeline.
;   This program is called once by mdrp_reduce3d for each IFU.
;
; CALLING SEQUENCE:
;   mdrp_reduceoneifu
;
; INPUTS:
;   framenames: String array containing full filepaths to the input frames
;   plan: mgPlan3d structure describing the exposures to combine
;   ifuDesign: The one ifu to combine data for (e.g., 12701).
;
; OPTIONAL INPUTS:
;   outdir: Output directory, by default $MANGA_SPECTRO_REDUX/DRPVER/PLATE/stack/
;   outnameroot: Root output name, by default manga-$plateid-$ifudesign-
;         Program will scan this string looking for '$' symbols.  If none found,
;         it simply uses this string as the root.  If found, it replaces these
;         elements with proper variable values.  Requires variables to be '-' delimited.
;   cubescale: Output spaxel size in cube in arcsec (default 0.5 arcsec/spaxel)
;   extast_log: Output logfile for the extended astrometry module
;
; OPTIONAL KEYWORDS:
;   /noeam: Extended astrometry module (EAM) broadband registration of
;         frames is on by default if the reference files can be found.
;         This switch turns it off.
;
; OUTPUT:
;
; COMMENTS:
;
; EXAMPLES:
;
; BUGS:
;   This routine spawns the Unix command 'mkdir'.
;
; PROCEDURES CALLED:
;
; INTERNAL SUPPORT ROUTINES:
;
; REVISION HISTORY:
;   October-2013  Written by David Law (Dunlap Institute; drlaw@di.utoronto.ca)
;   19-Mar-2014 Updated for production run (D. Law)
;   02-Oct-2014 Make default spaxel scale 0.5 arcsec (D. Law)
;   12-Nov-2014 Check/reject for dropped IFUs in a set (D. Law; dlaw@stsci.edu)
;   08-Dec-2014 Use EAM by default, change keywords (D. Law)
;   05-Mar-2015 Add srvymode keyword to identify bright-time plates (D. Law)
;   31-Mar-2015 Add foreground star flagging (D. Law)
;   24-Jun-2015 Add computation of LINEAR format files (D. Law)
;   15-Jul-2015 Add FLUXEXT, ERREXT, MASKEXT, ERRTYPE keyword headers (D. Law)
;   15-Jul-2015 Bug fix to linear headers, remove X/Y arr from cubes (D. Law)
;   16-Jul-2015 Eliminate multi-component elements in obsinfo (D. Law)
;   27-Aug-2015 Eliminate u-band image, psf calculations (D. Law)
;   14-Mar-2016 Fix major bug in psfparam creating reference PSF (D. Law)
;-
;------------------------------------------------------------------------------

pro mdrp_reduceoneifu, framenames, plan, ifuDesign, outdir=outdir, outnameroot=outnameroot1, noeam=noeam, logonly=logonly, linonly=linonly, dominimal=dominimal, cubescale=cubescale, extast_log=extast_log, srvymode=srvymode

; Default error behavior is to stop immediately
on_error,0
; Use 32-bit integers by default, and strict array usage
compile_opt idl2

; Initialize drp3qual to 0
drp3qual=0L

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Setup- define number of exposures, output directory, and
; partial slitmap for the first exposure (i.e., slitmap cropped
; to only the ifuDesign value)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

nframes=n_elements(framenames)

; Default output spaxel size is 0.5 arcsec
if (keyword_set(cubescale)) then cubescale=cubescale $
else cubescale=0.5

; Check that all frames exist
if (total(file_test(framenames)) ne nframes) then $
  message,'Frames not found!'

; Check that the number of exposures matches the plan file
if (nframes ne n_elements(plan)) then $
  message, 'Error: Number of exposures doesnt match plan file!'

splog,'nframes = '+strcompress(string(nframes),/remove_all)

; Check that ifuDesign has either a single value or 1 value for each exposure
ndesign=n_elements(ifuDesign)
if ((ndesign ne 1)and(ndesign ne nframes)) then $
  message, 'ERROR: Too many or too few values of ifuDesign'

; Read in slitmap for first exposure
ml_mgframeread,framenames[0],slitmap=slitmap
; partialslit is the slitmap entry for the 1st exposure for this IFU only
partialslit=slitmap[where(slitmap.ifuDesign eq ifuDesign[0])]

; Set output directory using provided filepath, otherwise construct it
; based on MaNGA data model
if (keyword_set(outdir)) then outdir=outdir $
else begin
  outdir=concat_dir(ml_getenv('MANGA_SPECTRO_REDUX'),mangadrp_version(/simple))
  outdir=concat_dir(outdir,strcompress(string(plan[0].plate),/remove_all))
  outdir=concat_dir(outdir,'stack/')
endelse

splog,'Output directory set to: '+outdir
; If output directory doesn't exist, create it and any necessary parent directories
spawn, 'mkdir -p ' + outdir

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Setup- define output filename root
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   
; Set the output filename root (manga-$plateid-$ifudesign-)
if (keyword_set(outnameroot1)) then outnameroot=outnameroot1 $
else outnameroot='manga-$plateid-$ifudesign-'

; If outnameroot has any variable names in it (indicated by $) replace them
; now with actual values.  NOTE!  This needs to be explicitly done for
; each possible variable.  Requires variables to be '-' delimited.
pieces=strsplit(outnameroot,'-',/extract)
npiece=n_elements(pieces)
for i=0,npiece-1 do begin
  ; If this piece starts with a $ symbol, remove the symbol to get variable name
  if (strpos(pieces[i],'$') eq 0) then begin
    varname=strmid(pieces[i],1)
    ; Depending on varname, change to the correct variable value
    if (varname eq 'plateid') then pieces[i]=strcompress(string(plan[0].plate),/remove_all)
    if (varname eq 'ifudesign') then pieces[i]=strcompress(string(ifuDesign[0]),/remove_all)
    if (varname eq 'mangaid') then pieces[i]=strcompress(string(partialslit[0].mangaID),/remove_all)
  endif
endfor
; Reassemble the pieces into the final output filename root
outnameroot=''
for i=0,npiece-1 do outnameroot+=pieces[i]+'-'

splog,'Outname root set to: '+outnameroot

; If not provided with an extended astrometry filename, set it here
if (~keyword_set(extast_log)) then extast_log=outdir+'eam.out'

; Set directory for qa plots
qaplotdir=concat_dir(outdir,'qa/')
; Create directory if needed
spawn, 'mkdir -p ' + qaplotdir
; Define qa filename
qaplotfile=qaplotdir+outnameroot+'qa.ps'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Setup- define total number of fibers, preimaging directory
; Make sure to reject sets in which one or more exposures had
; a dropped ifu or a bogey.  Overall set QA is done later.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Read in the list of bogeys
bogey=yanny_readone(concat_dir(ml_getenv('MANGACORE_DIR'),'apocomplete/bogey.par'),'BOGEY')

; Set a rejection vector.  If an exposure is to be rejected b/c
; of a dropped ifu, it will be set here.
rejectexp=intarr(nframes)

; Determine total number of fibers across all exposures by
; reading in each of the slitmaps.  In principle, there's no
; reason why this can't handle different physical ifus, with a
; different number of fibers, for each exposure.

; Loop over exposures the first time to determine accept/reject
nfiber=intarr(nframes)
splog,'Counting fibers...'
for i=0,nframes-1 do begin
  ; Allow for multiple values of ifuDesign
  if (ndesign eq 1) then thisdesign=ifuDesign $
  else thisdesign=ifuDesign[i]

  ; Read in the slitmap and FITS header for this exposure
  ml_mgframeread,framenames[i],slitmap=slitmap,hdr=hdr
  ; Crop out lines appropriate for this IFU
  partialslit=slitmap[where(slitmap.ifuDesign eq thisdesign)]
  nfiber[i]=n_elements(partialslit)

  if nfiber[i] eq 0 then $
    message, 'ERROR- ifuDesign= '+strcompress(string(ifuDesign),/remove_all)+' not found in '+framenames[i]

  ; Check for bogeys against the master list
  thisexp=sxpar(hdr,'EXPOSURE')
  index=where((bogey.exposure eq thisexp)and(bogey.ifudesign eq long(thisdesign)),nbogey)
  ; If there's a bogey, this exp and all of the other exposures in the set are bad (need uniform sets)
  if (nbogey ne 0) then begin
    index=where(plan.set eq plan[i].set)
    rejectexp[index]=1
    splog,'Bogey in frame '+string(i)+' rejecting!'
  endif

  ; Check for dropped IFUs using DRP2QUAL, 'BADIFU' bit
  drp2qual=long(sxpar(hdr,'DRP2QUAL'))
  ; Were there any bad ifus?
  if ((drp2qual and sdss_flagval('MANGA_DRP2QUAL','BADIFU')) ne 0) then begin
    ; Is THIS ifu bad?
    if (stregex(sxpar(hdr,'BADIFU'),strcompress(partialslit[0].harname,/remove_all),/fold_case) ge 0) then begin
      ; This exp and all of the other exposures in the set are bad (need uniform sets)
      index=where(plan.set eq plan[i].set)
      rejectexp[index]=1
    endif
  endif
endfor

; Loop over exposures again to do the actual fiber counting
nftotal=0
for i=0,nframes-1 do begin
  splog, 'Number of fibers in exposure '+strcompress(string(i+1),/remove_all)+' is '+strcompress(string(nfiber[i]),/remove_all)

  ; Accept fibers
  if (rejectexp[i] eq 0) then nftotal+=nfiber[i] $
  ; Reject fibers
  else begin
    splog,'REJECTING this exposure (likely dropped ifu)'
  endelse
endfor

; Trim things to only accepted exposures
index=where(rejectexp eq 0,ngood)
; Failsafe
if (ngood eq 0) then begin
  splog,'No good exposures, quit!'
  return
endif else begin
  nframes=ngood
  nfiber=nfiber[index]
  finalframenames=framenames[index]
  finalplan=plan[index]
  if (ndesign gt 1) then finalifuDesign=ifuDesign[index] $
  else finalifuDesign=ifuDesign
endelse

splog,strcompress(string(nftotal),/remove_all)+' fibers total.'

; Define the preimaging directory.
preim_dir=concat_dir(ml_getenv('MANGAPREIM_DIR'),'data')
; If it's not defined, then turn off EAM
if (preim_dir eq '') then begin
  splog,'MANGAPREIM_DIR not defined, turning off EAM.'
endif

; Define the blank observing info structure.  This is the .par structure
; that consolidates key information about what went into each cube
; in a manner that can be read by the data simulator.
; We'll add one line to this structure later for each exposure
; combined.
obsinfo_line=ml_createobsstruc()

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Now we've got multiple sets of data to process.
; We have nframes worth of data for 3 data flavors:
; LINEAR: Flux calibrated data on a linear wave grid
; LOG: Flux calibrated data on a log wave grid
;
; Do the LOG loop first.  We'll work out astrometry here and
; apply it to the LINEAR loop later.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Data has the MaNGA log wavelength solution
wave=ml_setwcalib()
nwave=n_elements(wave)
; Figure out indices where the wavelength solution best matches
; SDSS g,r,i,z broadband effective lambda
ml_sdssfilters,gpivot=gpivot,rpivot=rpivot,ipivot=ipivot,zpivot=zpivot
bbwave=[gpivot,rpivot,ipivot,zpivot]
bbindex=intarr(n_elements(bbwave))
for j=0,n_elements(bbwave)-1 do begin
  temp=min(abs(wave-bbwave[j]),minat)
  bbindex[j]=minat
endfor

; Set up arrays for the data
flux_arr=fltarr(nwave,nftotal)
ivar_arr=fltarr(nwave,nftotal)
mask_arr=lonarr(nwave,nftotal)
disp_arr=fltarr(nwave,nftotal)
x_arr=fltarr(nwave,nftotal)
y_arr=fltarr(nwave,nftotal)

; Set up arrays for the psf reference
psf_flux=fltarr(n_elements(bbwave),nftotal)
psf_xarr=fltarr(n_elements(bbwave),nftotal)
psf_yarr=fltarr(n_elements(bbwave),nftotal)

; Call the PSF code for the first frame only to determine
; dimensionality of the psfplane array returned by fcpsf1d
if (ndesign eq 1) then thisdesign=finalifuDesign $
else thisdesign=finalifuDesign[0]
ml_mgframeread,finalframenames[0],slitmap=slitmap,hdr=hdr
temp=fileandpath(finalframenames[0],path=cogdir)
psfpar=psfparams(long(fxpar(hdr,'PLATEID')),long(fxpar(hdr,'MJD')),long(fxpar(hdr,'EXPOSURE')),reduxdir=cogdir)
psffac=float(fxpar(hdr,'PSFFAC'))
psffac=(psffac > 0.5) < 1.0
index=where(slitmap.ifuDesign eq thisdesign)
ml_mgframeread,finalframenames[0],index,slitmap=partialslit
fcpsf1d,psfpar.params,partialslit[0].xfocal,partialslit[0].yfocal,psfplane,wavearr=bbwave,xarr=xpsfplane,fluxcorr=psfpar.fluxcorr,psffactor=psffac
; Define a couple of arrays to hold psfplane and xpsfplane
; for each exposure
allpsfplane=fltarr([size(psfplane,/dim),nframes])
allxpsfplane=fltarr([size(xpsfplane,/dim),nframes])

; Running index for which fiber in nftotal a frame starts on
fiberstart=0L

; Loop over individual exposures
for i=0,nframes-1 do begin
  splog,''
  splog,strcompress('Processing exposure number '+string(i+1))
  splog,'File = '+finalframenames[i]

  ; Allow for multiple values of finalifuDesign
  if (ndesign eq 1) then thisdesign=finalifuDesign $
  else thisdesign=finalifuDesign[i]

  ; Read slitmap from file first
  ml_mgframeread,finalframenames[i],slitmap=slitmap,slithdr=slithdr

  ; index is the lines from slitmap appropriate for this IFU
  index=where(slitmap.ifuDesign eq thisdesign)

  ; Read data from file for only these ifu lines only
  ml_mgframeread,finalframenames[i],index,objflux=objflux,objivar=objivar, mask=objmask, $
    dispimg=dispimg, slitmap=partialslit, hdr=hdr

  ; If this is the first frame, save the FITS header
  if (i eq 0) then hdrsave=hdr

  ; Number of fibers read in from this exposure
  nfiber=n_elements(partialslit)

  ; Check that loglam (the wavelength solution from the data file)
  ; has nwave elements
  if ((size(objflux))[1] ne nwave) then $
    message,'ERROR- Incorrect number of wavelength samples.'
  ; DRL- check that it matches default wave solution too???
  ; Should this be generalized?

  ; Figure out the simplistic on-sky PSF fwhm for this IFU at u,g,r,i,z bands
  ; taking into account focal offset as a function of wavelength and plate location
  ; First get the basic psf from the guider images
  ; Directory for the cogimg guider images
  temp=fileandpath(finalframenames[i],path=cogdir)
  psfpar=psfparams(long(fxpar(hdr,'PLATEID')),long(fxpar(hdr,'MJD')),long(fxpar(hdr,'EXPOSURE')),reduxdir=cogdir)
  ; Compute the effective fwhm at these wavelengths, for this plate position
  ; Shrink by some percentage from the guider values, which are
  ; always high.  Use PSFFAC derived from minibundles in flux cal
  ; routine to determine appropriate percentage
  psffac=float(fxpar(hdr,'PSFFAC'))
  ; Look for catastrophic psffac failures, require in range 0.5 - 1
  if ((psffac lt 0.5)or(psffac gt 1.0)) then $
    drp3qual=drp3qual OR sdss_flagval('MANGA_DRP3QUAL','BADPSF')
  psffac=(psffac > 0.5) < 1.0
  fcpsf1d,psfpar.params,partialslit[0].xfocal,partialslit[0].yfocal,psfplane,wavearr=bbwave,xarr=xpsfplane,fwhmout=fwhm_eff,/dofwhm,fluxcorr=psfpar.fluxcorr,psffactor=psffac
  ; Save the PSF information into allpsfplane and allxpsfplane
  ; (will be used again in the second run of the EAM)
  allpsfplane[*,*,i]=psfplane
  allxpsfplane[*,i]=xpsfplane

  ; Populate this line of obsinfo with basic relevant information
  ml_addobsinfo,obsinfo_line,partialslit,slithdr,hdr,pf_fwhm=fwhm_eff
  ; Add this line to the obsinfo file
  if (~keyword_set(obsinfo)) then obsinfo=obsinfo_line $
  else obsinfo=[obsinfo,obsinfo_line]
  ; Define the set membership
  obsinfo[i].set=finalplan[i].set
  ; Record psffac
  obsinfo[i].psffac=psffac

  ; Combine all previous 2d maskbits into the mask array
  newmask=objmask

  ; Decide what 2d maskbits are important enough that these pixels will
  ; be rejected when constructing the cube.  For easy reference, add
  ; the '3DREJECT' bit for each of these.
  ; Look for where the MANGA_DRP2PIXMASK bits 'FULLREJECT', 'COMBINEREJ',
  ; 'LOWFLAT', or 'NODATA' are set.  Set these pixels to zero inverse variance
  ; and activate the '3DREJECT' bit
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','COMBINEREJ')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','FULLREJECT')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','LOWFLAT')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','NODATA')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif

  ; Figure out where the broken fibers were using the slitmap.gbu flag
  ; and add this information to the MANGA_DRP2PIXMASK bitmask
  index=where(partialslit.gbu ne 0,nindex)
  if (nindex ne 0) then begin
    tempmask=replicate(0L,nfiber)
    tempmask[index]=sdss_flagval('MANGA_DRP2PIXMASK','DEADFIBER')
    newmask = newmask OR rebin(reform(tempmask,1,nfiber),nwave,nfiber)
  endif

  ; Look for where the inverse variance is zero, and zero out the
  ; flux at these locations too so that they don't confuse anyone
  index=where(objivar eq 0.,nindex)
  if (nindex gt 0) then objflux[index]=0.

  ; Stuff the flux, inverse variance, and new bitmask into the big arrays
  flux_arr[*,fiberstart:fiberstart+nfiber-1]=objflux[*,*]
  ivar_arr[*,fiberstart:fiberstart+nfiber-1]=objivar[*,*]
  mask_arr[*,fiberstart:fiberstart+nfiber-1]=newmask[*,*]
  disp_arr[*,fiberstart:fiberstart+nfiber-1]=dispimg[*,*]

  ; Work out the basic astrometry for these fiber positions and stick them in
  ; the vectors xfiber and yfiber on-sky.  If this is the first exposure,
  ; use the /setbase keyword to mark this as the zeropoint against which
  ; DC shifts should be measured
  if (i eq 0) then begin
    mdrp_astrometry,xfiber,yfiber,partialslit,slithdr,hdr,wave
  endif else begin
    mdrp_astrometry,xfiber,yfiber,partialslit,slithdr,hdr,wave
  endelse

  ; Run the extended astrometry module to tweak up positions, so long
  ; as the preimaging directory is defined and we aren't in 
  ; APOGEE-led mode (i.e., looking at stars)
  if ((~keyword_set(noeam))and(keyword_set(preim_dir))and(srvymode ne 'APOGEE')) then begin
    ; Convert from x,y fiber positions in arcsec into RA,DEC
    rafiber=partialslit[0].ra-xfiber/3600.D/cos(partialslit[0].dec*!DPI/180.)
    decfiber=partialslit[0].dec+yfiber/3600.D

    ; Construct filepath to the reference image
    refimagepath=concat_dir(preim_dir,ml_plategroup(fxpar(hdr,'DESIGNID'),/design))+'/'+strcompress(string(fxpar(hdr,'DESIGNID')),/remove_all)+'/'
    refimagefile=lookforgzip(refimagepath+'preimage-'+strcompress(string(partialslit[0].mangaID),/remove_all)+'.fits')

    ; Can we find this file?  If so, do the astrometric tweaking
    if (file_test(refimagefile)) then begin
      ; Set up the output qa plot path
      frametemp=fileandpath(finalframenames[i],path=pathtemp)
      qadir=concat_dir(pathtemp,'qa/')
      ; Create qa directory if it doesn't exist
      if file_test(qadir,/directory) eq 0 then spawn, '\mkdir -p '+qadir
      ; Construct full path to output qa file
      outplot=concat_dir(qadir,(strsplit(frametemp, '.fits', /regex, /extract))[0]+'_'+thisdesign+'_eam')

      ; Run the module.
      ; rafiber, decfiber, objflux, objivar, objmask all have dimension [nwave, nfiber]
      ; wave is a vector of size [nwave]
       mdrp_eam, rafiber, decfiber, objflux, objivar, objmask, wave, [obsinfo[i].pf_fwhm_g,obsinfo[i].pf_fwhm_r,obsinfo[i].pf_fwhm_i,obsinfo[i].pf_fwhm_z], refimagefile, raout=raout, decout=decout, fluxout=fluxout, errorout=errorout, plan=finalplan[i], ifudesign=thisdesign, outputfile=extast_log, /plot, outplot=outplot, outpar=outpar

      ; Save derived values for dTHETA0 (theta when it was free)
      obsinfo[i].EAMfit_theta0=outpar[2]
      obsinfo[i].EAMfit_theta0err=outpar[7]
    endif
  endif

  ; Stuff the x,y fiber relative position information into the big arrays
  x_arr[*,fiberstart:fiberstart+nfiber-1]=xfiber[*,*]
  y_arr[*,fiberstart:fiberstart+nfiber-1]=yfiber[*,*]

  ; Work out the psf information.  Using the actual xfiber,yfiber info,
  ; figure out the flux that would have gone into each from a point
  ; source convolved with true on-sky profile (pre-fiber).
  ; Do this for 4 wavelengths g,r,i,z
  for j=0,n_elements(bbindex)-1 do begin
    psf_xarr[j,fiberstart:fiberstart+nfiber-1]=xfiber[bbindex[j],*]
    psf_yarr[j,fiberstart:fiberstart+nfiber-1]=yfiber[bbindex[j],*]

    roff=reform(sqrt(xfiber[bbindex[j],*]^2+yfiber[bbindex[j],*]^2),nfiber)
    ftemp=replicate(0.,nfiber)
    ; Only calculate flux in fibers within 6'' of center
    index=where(roff lt 6.,nindex)
    for k=0,nindex-1 do $
      ; Linear interpolation
      ftemp[index[k]]=interpol(psfplane[*,j],xpsfplane,roff[index[k]])

    psf_flux[j,fiberstart:fiberstart+nfiber-1]=ftemp[*]
  endfor

  ; Increment the fiber running index
  fiberstart+=nfiber
endfor

; Loop over all pluggings to calculate the bundle rotation theta
; on average for this plugging.
theplugs=obsinfo[uniq(obsinfo.slitfile)].slitfile
nplug=n_elements(theplugs)
; Define a new vector to contain the mean theta values
newtheta=fltarr(nframes)
for i=0,nplug-1 do begin
  theseexp=where(obsinfo.slitfile eq theplugs[i])
  newtheta[theseexp]=ml_meanclip(obsinfo[theseexp].EAMfit_theta0)
endfor

; Loop over all exposures and run the EAM again, this time with
; fixed theta values and actually apply the results to our
; astrometry for the fibers.
if ((~keyword_set(noeam))and(keyword_set(preim_dir))and(srvymode ne 'APOGEE')) then begin
  fiberstart=0L
  splog,'Rerunning EAM with fixed theta if applicable'
  for i=0,nframes-1 do begin

    ; Allow for multiple values of finalifuDesign
    if (ndesign eq 1) then thisdesign=finalifuDesign $
    else thisdesign=finalifuDesign[i]
    ; Read slitmap from file first
    ml_mgframeread,finalframenames[i],slitmap=slitmap,slithdr=slithdr
    ; index is the lines from slitmap appropriate for this IFU
    index=where(slitmap.ifuDesign eq thisdesign,nfiber)
    partialslit=slitmap[index]

    ; We can read information out of the big arrays instead of having
    ; to re-read it from disk and repeat things from before
    objflux=flux_arr[*,fiberstart:fiberstart+nfiber-1]
    objivar=ivar_arr[*,fiberstart:fiberstart+nfiber-1]
    objmask=mask_arr[*,fiberstart:fiberstart+nfiber-1]
    xfiber=x_arr[*,fiberstart:fiberstart+nfiber-1]
    yfiber=y_arr[*,fiberstart:fiberstart+nfiber-1]
    ; Convert from x,y fiber positions in arcsec into RA,DEC
    rafiber=partialslit[0].ra-xfiber/3600.D/cos(partialslit[0].dec*!DPI/180.)
    decfiber=partialslit[0].dec+yfiber/3600.D

    ; Construct filepath to the reference image
    refimagepath=concat_dir(preim_dir,ml_plategroup(fxpar(hdr,'DESIGNID'),/design))+'/'+strcompress(string(fxpar(hdr,'DESIGNID')),/remove_all)+'/'
    refimagefile=lookforgzip(refimagepath+'preimage-'+strcompress(string(partialslit[0].mangaID),/remove_all)+'.fits')
    ; Can we find this file?  If so, do the astrometric tweaking
      if (file_test(refimagefile)) then begin
        ; Make sure to call this with /plot so that it doesn't try to
        ; send plots to xdisplay at Utah, but we dont' need them kept
        ; so don't specific output plot name
        mdrp_eam, rafiber, decfiber, objflux, objivar, objmask, wave, [obsinfo[i].pf_fwhm_g,obsinfo[i].pf_fwhm_r,obsinfo[i].pf_fwhm_i,obsinfo[i].pf_fwhm_z], refimagefile, raout=raout, decout=decout, fluxout=fluxout, errorout=errorout, plan=finalplan[i], ifudesign=thisdesign, outputfile=extast_log, outpar=outpar, FIXROT=newtheta[i],/plot

        ; Convert from RA,DEC back to x,y fiber offsets in arcsec relative to fiducial pointing
        xfiber=-(raout-partialslit[0].ra)*3600.D*cos(partialslit[0].dec*!DPI/180.)
        yfiber=(decout-partialslit[0].dec)*3600.D
        ; Stuff the x,y fiber relative position information into the big arrays
        x_arr[*,fiberstart:fiberstart+nfiber-1]=xfiber[*,*]
        y_arr[*,fiberstart:fiberstart+nfiber-1]=yfiber[*,*]

      ; Save derived values for dRA, dDEC, dTHETA, A, B (mean combined across bands)
      ; Note that A,B are flux scaling and offset: while these aren't used
      ; they are useful for QA.
      obsinfo[i].EAMfit_ra=outpar[0]
      obsinfo[i].EAMfit_dec=outpar[1]
      obsinfo[i].EAMfit_theta=outpar[2]
      obsinfo[i].EAMfit_a=outpar[3]
      obsinfo[i].EAMfit_b=outpar[4]
      obsinfo[i].EAMfit_raerr=outpar[5]
      obsinfo[i].EAMfit_decerr=outpar[6]
      obsinfo[i].EAMfit_thetaerr=outpar[7]
      obsinfo[i].EAMfit_aerr=outpar[8]
      obsinfo[i].EAMfit_berr=outpar[9]
      endif

    ; Restore psfplane information from allpsfplane array
    psfplane=allpsfplane[*,*,i]
    xpsfplane=allxpsfplane[*,i]
    ; Revise PSF information based on updated positions
    for j=0,n_elements(bbindex)-1 do begin
      psf_xarr[j,fiberstart:fiberstart+nfiber-1]=xfiber[bbindex[j],*]
      psf_yarr[j,fiberstart:fiberstart+nfiber-1]=yfiber[bbindex[j],*]

      roff=reform(sqrt(xfiber[bbindex[j],*]^2+yfiber[bbindex[j],*]^2),nfiber)
      ftemp=replicate(0.,nfiber)
      ; Only calculate flux in fibers within 6'' of center
      index=where(roff lt 6.,nindex)
      for k=0,nindex-1 do $
        ; Linear interpolation
        ftemp[index[k]]=interpol(psfplane[*,j],xpsfplane,roff[index[k]])

      psf_flux[j,fiberstart:fiberstart+nfiber-1]=ftemp[*]
    endfor

    ; Increment the fiber running index
    fiberstart+=nfiber
  endfor
endif

; Loop over sets to calculate Omega
if (srvymode ne 'APOGEE') then begin
  thesets=obsinfo[uniq(obsinfo.set)].set
  nsets=n_elements(thesets)
  for i=0,nsets-1 do begin
    theseexp=where(obsinfo.set eq thesets[i])
    allomega=mdrp_omega(obsinfo[theseexp])
    obsinfo[theseexp].omegaset_u=allomega[0]
    obsinfo[theseexp].omegaset_g=allomega[1]
    obsinfo[theseexp].omegaset_r=allomega[2]
    obsinfo[theseexp].omegaset_i=allomega[3]
    obsinfo[theseexp].omegaset_z=allomega[4]
  endfor
endif

; Quality control plots and flagging based on the obsinfo structure.
; Don't bother if srvymode=APOGEE cuz this will look bad since not dithered
if (srvymode ne 'APOGEE') then $
  ml_qaobsinfo,obsinfo,filename=qaplotfile,drp3qual=drp3qual

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Create data cube in log wavelength grid
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Cube size in arcsec based on relative fiber positions
; across all exposures
cube_zsize=nwave
; Require square format output
cube_xysize=fix(max(x_arr)-min(x_arr)) > fix(max(y_arr)-min(y_arr))

; For consistency make Odd integers true
; Convert to output pixel scale and add some padding
cube_xysize = fix(cube_xysize/cubescale) +10
; For consistency make the number of pixels even
; (Note that odd integers count as TRUE)
if (cube_xysize) then cube_xysize+=1
; Assign to both x and y
cube_xsize = cube_xysize
cube_ysize = cube_xysize

; Create data cube and inverse variance structure
cube_flux=fltarr(cube_xsize,cube_ysize,cube_zsize)
cube_ivar=fltarr(cube_xsize,cube_ysize,cube_zsize)
cube_mask=lonarr(cube_xsize,cube_ysize,cube_zsize)

; Define 1d vectors of x,y,flux,ivar,mask to be used
; at each wavelength slice
xvec=fltarr(nftotal)
yvec=fltarr(nftotal)
fvec=fltarr(nftotal)
ivec=fltarr(nftotal)
mvec=lonarr(nftotal)

; Define output psf cube
psfcube_xsize=cube_xsize*1.
psfcube_ysize=cube_ysize*1.
psfcube_scale=cubescale/1.
psfcube=fltarr(psfcube_xsize,psfcube_ysize,n_elements(bbwave))
; FWHM of single-gaussian fit to resulting cube
cubefwhm=fltarr(n_elements(bbwave))

; Figure out reconstructed psf in each of the 4 wavebands
for i=0,n_elements(bbwave)-1 do begin
  xvec=psf_xarr[i,*]/psfcube_scale+psfcube_xsize/2.
  yvec=psf_yarr[i,*]/psfcube_scale+psfcube_ysize/2.
  fvec=psf_flux[i,*]
  scale=psfcube_scale*psfcube_scale/!PI
  psfcube[*,*,i]=ml_griddata(xvec,yvec,fvec,[psfcube_xsize,psfcube_ysize],1.6/psfcube_scale,0.7/psfcube_scale,scale=scale)
  temp=gauss2dfit(psfcube[*,*,i],coef)
  cubefwhm[i]=mean(coef[2:3])*2.355*psfcube_scale
  ; Safety case just in the event that gauss2dfit returned bogus values
  ; Call <0'' or > 10'' bogus and set them to 0''
  if ((cubefwhm[i] lt 0.)or(cubefwhm[i] gt 10.)) then cubefwhm[i]=0.
endfor

; Loop over all wavelengths figuring out the cube at that wavelength
; (unless we're in bright-time mode, in which case don't do this loop;
; cubes will be just empty placeholders.  This saves ENORMOUS cpu time.)
if (srvymode ne 'APOGEE') then begin
  for i=0,cube_zsize-1 do begin
    ; Print a message to splog every 10% of the loop
    if (i mod fix(cube_zsize/10) eq 0) then $
      splog,'Constructing cube: '+strcompress(string(round(i*100./cube_zsize)),/remove_all)+'% complete'

    ; Create input vectors for this wavelength slice
    xvec[*]=x_arr[i,*]/cubescale+cube_xsize/2.
    yvec[*]=y_arr[i,*]/cubescale+cube_ysize/2.
    fvec[*]=flux_arr[i,*]
    ivec[*]=ivar_arr[i,*]
    mvec[*]=mask_arr[i,*]

    ; Scale correction factor is the ratio of area between a PI area fiber
    ; (in arcsec^2) and the output pixel size in arcsec^2
    ; The result means that the cube will be in calibrated units/pixel
    scale=cubescale*cubescale/!PI

;if (i eq 2545) then stop
;if (i eq 1312) then stop

    cube_flux[*,*,i]=ml_griddata(xvec,yvec,fvec,[cube_xsize,cube_ysize],1.6/cubescale,0.7/cubescale,ivar=ivec,invarimg=invarimg,maskvec=mvec,maskimg=maskimg,scale=scale)
    ; Populate inverse variance and mask arrays
    cube_ivar[*,*,i]=invarimg
    cube_mask[*,*,i]=maskimg
  endfor
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Work out the median spectral resolving power
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Median across all fibers in this IFU
disp_med=median(disp_arr,dimension=2)

; Bspline fits that will be used in place of any zero values
; (e.g., if NODATA for some reason)
; Weight anything with disp_med=0 with 0 weight
tmpwght=replicate(1.,nwave)
index=where(disp_med eq 0.,nindex)
if (nindex gt 0) then tmpwght[index]=0.
nord=3
everyn=100
bkpt=0
sset=bspline_iterfit(wave,disp_med,nord=nord,everyn=everyn,bkpt=bkpt,maxrej=0,invvar=tmpwght,yfit=yfit)
if (nindex gt 0) then disp_med[index]=yfit[index]


; Figure out variability at each lambda between
; fibers and between exposures (i.e., variability
; along the slit, and from exp to exposure)
specvar=fltarr(nwave)
for i=0,nwave-1 do begin
  ; Sigma clip to take out lone bad values
  temp=ml_meanclip(disp_arr[i,*],themean,thesig,ignore=0.)
  ; Fractional RMS variability of the LSF
  specvar[i]=thesig/themean
endfor

; Bspline fits that will be used in place of any zero values
; (e.g., if NODATA for some reason)
; Weight anything with disp_med=0 with 0 weight
tmpwght=replicate(1.,nwave)
index=where(specvar eq 0.,nindex)
if (nindex gt 0) then tmpwght[index]=0.
nord=3
everyn=100
bkpt=0
sset=bspline_iterfit(wave,specvar,nord=nord,everyn=everyn,bkpt=bkpt,maxrej=0,invvar=tmpwght,yfit=yfit)
if (nindex gt 0) then specvar[index]=yfit[index]

; Mean spectral resolution R=lambda/(FWHM)
specres=wave/(disp_med*2.35)
; Standard deviation of spectral resolution about this mean
specres_stddev=specres*specvar

; QA flag VARIABLELSF only in true failure cases where
; the LSF varies by more than 50% across 5% or more of the
; wavelength range
temp=where(specvar gt 0.5,nhighvar)
if (float(nhighvar)/nwave gt 0.05) then $
  drp3qual=drp3qual OR sdss_flagval('MANGA_DRP3QUAL','VARIABLELSF')

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Make the basic cube and RSS file headers
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

mkhdr,rsshdr,flux_arr,/image
mkhdr,cubehdr,cube_flux,/image
; Add MaNGA keycards
mdrp_drp3dhead,drp3qual=drp3qual,cubehdr=cubehdr,rsshdr=rsshdr,obsinfo=obsinfo,mgcf_hdr=hdrsave,wave=wave,cubescale=cubescale,cubefwhm=cubefwhm,cube0hdr=cube0hdr,rss0hdr=rss0hdr

; Make headers for the ivar and mask extensions
mkhdr,rss_ivar_hdr,ivar_arr,/image
mkhdr,rss_mask_hdr,mask_arr,/image
mkhdr,cube_ivar_hdr,cube_ivar,/image
mkhdr,cube_mask_hdr,cube_mask,/image

; Add info to the global headers to indicate flux, error, mask extensions
mdrp_hduclass,cubehdr,errhdr=cube_ivar_hdr,maskhdr=cube_mask_hdr,type='CUBE'
mdrp_hduclass,rsshdr,errhdr=rss_ivar_hdr,maskhdr=rss_mask_hdr,type='IMAGE'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Quality check- look for known foreground stars in the data cube
; so that we can flag them
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

thismangaid=strcompress(fxpar(cubehdr,'MANGAID'),/remove_all)
knownstars=yanny_readone(concat_dir(ml_getenv('MANGACORE_DIR'),'platedesign/foregroundstars/foregroundstars.par'),'FORESTAR')
index=where(strcompress(knownstars.mangaid,/remove_all) eq thismangaid,nindex)
starmask=intarr(cube_xsize,cube_ysize)

for k=0,nindex-1 do begin
  maskra=knownstars[index[k]].raj2000
  maskdec=knownstars[index[k]].dej2000
  maskradius=knownstars[index[k]].radius
  starmask[*,*]=0L
  for i=0,cube_xsize-1 do begin
    for j=0,cube_ysize-1 do begin
      xyad,cubehdr,i,j,thisra,thisdec
      deltara=3600.*(maskra-thisra)*cos(thisdec*!PI/180.)
      deltadec=3600.*(maskdec-thisdec)
      if (sqrt((deltara*deltara)+(deltadec*deltadec)) lt maskradius) then starmask[i,j]=1
    endfor
  endfor
  starmask=starmask*sdss_flagval('MANGA_DRP3PIXMASK','FORESTAR')
  for q=0,cube_zsize-1 do cube_mask[*,*,q]=cube_mask[*,*,q] or starmask
endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Make some broadband images from the data cube
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Pass a version of ivar that has been zeroed everywhere that
; the mask is DONOTUSE
tempivar=cube_ivar
index=where((cube_mask and sdss_flagval('MANGA_DRP3PIXMASK','DONOTUSE')) ne 0,nindex)
if (nindex gt 0) then tempivar[index]=0.

g_img=ml_mangatosdssimg(cube_flux,tempivar,wave,cubehdr,'g',imghdr=g_hdr)
r_img=ml_mangatosdssimg(cube_flux,tempivar,wave,cubehdr,'r',imghdr=r_hdr)
i_img=ml_mangatosdssimg(cube_flux,tempivar,wave,cubehdr,'i',imghdr=i_hdr)
z_img=ml_mangatosdssimg(cube_flux,tempivar,wave,cubehdr,'z',imghdr=z_hdr)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Write out LOG data files
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Write out the RSS-format file
rssname=outdir+outnameroot+'LOGRSS.fits'
ml_mwrfits, dummyext, rssname, hdr=rss0hdr, /create
ml_mwrfits, flux_arr, rssname, hdr=rsshdr, extname='FLUX'
ml_mwrfits, ivar_arr, rssname, hdr=rss_ivar_hdr, extname='IVAR'
ml_mwrfits, mask_arr, rssname, hdr=rss_mask_hdr, extname='MASK'
ml_mwrfits, disp_arr, rssname, extname='DISP'
ml_mwrfits, wave, rssname, extname='WAVE'
ml_mwrfits, specres, rssname, extname='SPECRES'
ml_mwrfits, specres_stddev, rssname, extname='SPECRESD'
ml_mwrfits, obsinfo, rssname, extname='OBSINFO'
ml_mwrfits, x_arr, rssname, extname='XPOS'
ml_mwrfits, y_arr, rssname, extname='YPOS'

; Write out the cube
cubename=outdir+outnameroot+'LOGCUBE.fits'
ml_mwrfits, dummyext, cubename, hdr=cube0hdr, /create
ml_mwrfits, cube_flux, cubename, hdr=cubehdr, extname='FLUX'
ml_mwrfits, cube_ivar, cubename, hdr=cube_ivar_hdr, extname='IVAR'
ml_mwrfits, cube_mask, cubename, hdr=cube_mask_hdr, extname='MASK'
ml_mwrfits, wave, cubename, extname='WAVE'
ml_mwrfits, specres, cubename, extname='SPECRES'
ml_mwrfits, specres_stddev, cubename, extname='SPECRESD'
ml_mwrfits, obsinfo, cubename, extname='OBSINFO'
ml_mwrfits, g_img, cubename, hdr=g_hdr, extname='GIMG'
ml_mwrfits, r_img, cubename, hdr=r_hdr, extname='RIMG'
ml_mwrfits, i_img, cubename, hdr=i_hdr, extname='IIMG'
ml_mwrfits, z_img, cubename, hdr=z_hdr, extname='ZIMG'
ml_mwrfits, psfcube[*,*,0], cubename, hdr=g_hdr, extname='GPSF'
ml_mwrfits, psfcube[*,*,1], cubename, hdr=r_hdr, extname='RPSF'
ml_mwrfits, psfcube[*,*,2], cubename, hdr=i_hdr, extname='IPSF'
ml_mwrfits, psfcube[*,*,3], cubename, hdr=z_hdr, extname='ZPSF'

; gzip output files
spawn, ['gzip', '-f', rssname], /noshell
spawn, ['gzip', '-f', cubename], /noshell

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;; Do the LINEAR format calculations ;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; We need to actually re-read all of the input files from the LIN
; instead of LOG mgCFrame files, and redo all of the bad pixel
; flagging and cube creation.  However, we can skip all of the
; astrometry since the astrometric solution in LOG wavelength
; space can easily be interpolated to the LINEAR space with
; no loss of accuracy.  We can also simply carry over or
; interpolate things like OBSINFO, spectral resolution, etc.

splog,'Processing LINEAR format data'

; Save memory by zeroing out much of the old arrays
flux_arr=0
ivar_arr=0
mask_arr=0
disp_arr=0
cube_flux=0
cube_ivar=0
cube_mask=0

; Data has the MaNGA linear wavelength solution
linwave=ml_setwcalib(/linear)
nlinwave=n_elements(linwave)

; Set up arrays for the data
flux_arr=fltarr(nlinwave,nftotal)
ivar_arr=fltarr(nlinwave,nftotal)
mask_arr=lonarr(nlinwave,nftotal)
disp_arr=fltarr(nlinwave,nftotal)
x_linarr=fltarr(nlinwave,nftotal)
y_linarr=fltarr(nlinwave,nftotal)

cube_zsize=nlinwave
cube_flux=fltarr(cube_xsize,cube_ysize,cube_zsize)
cube_ivar=fltarr(cube_xsize,cube_ysize,cube_zsize)
cube_mask=lonarr(cube_xsize,cube_ysize,cube_zsize)

; Running index for which fiber in nftotal a frame starts on
fiberstart=0L

; Construct the vector of LINEAR filenames
; and do a safety test that all of these can actually be found
; on the filesystem
linframenames=finalframenames
missing=0
for i=0,nframes-1 do begin
  linframenames[i]=ml_strreplace(finalframenames[i],'-LOG.fits','-LIN.fits')
  missing += (~(file_test(linframenames[i])))
endfor

; Only process exposures if none are missing
if (missing ne 0) then begin
  splog,'Cant find all LINEAR format files expected, skipping!'
endif else begin

; Loop over individual exposures
for i=0,nframes-1 do begin
  splog,''
  splog,strcompress('Processing exposure number '+string(i+1))
  splog,'File = '+linframenames[i]

  ; Allow for multiple values of finalifuDesign
  if (ndesign eq 1) then thisdesign=finalifuDesign $
  else thisdesign=finalifuDesign[i]

  ; Read slitmap from file first
  ml_mgframeread,linframenames[i],slitmap=slitmap,slithdr=slithdr

  ; index is the lines from slitmap appropriate for this IFU
  index=where(slitmap.ifuDesign eq thisdesign)

  ; Read data from file for only these ifu lines only
  ml_mgframeread,linframenames[i],index,objflux=objflux,objivar=objivar, mask=objmask, $
    dispimg=dispimg, slitmap=partialslit, hdr=hdr

  ; If this is the first frame, save the FITS header
  if (i eq 0) then hdrsave=hdr

  ; Number of fibers read in from this exposure
  nfiber=n_elements(partialslit)

  ; Check that loglam (the wavelength solution from the data file)
  ; has nlinwave elements
  if ((size(objflux))[1] ne nlinwave) then $
    message,'ERROR- Incorrect number of wavelength samples.'

  ; Combine all previous 2d maskbits into the mask array
  newmask=objmask

  ; Decide what 2d maskbits are important enough that these pixels will
  ; be rejected when constructing the cube.  For easy reference, add
  ; the '3DREJECT' bit for each of these.
  ; Look for where the MANGA_DRP2PIXMASK bits 'FULLREJECT', 'COMBINEREJ',
  ; 'LOWFLAT', or 'NODATA' are set.  Set these pixels to zero inverse variance
  ; and activate the '3DREJECT' bit
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','COMBINEREJ')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','FULLREJECT')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','LOWFLAT')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif
  index=where((objmask and sdss_flagval('MANGA_DRP2PIXMASK','NODATA')) ne 0,nindex)
  if (nindex ne 0) then begin
    objivar[index]=0.
    newmask[index]=newmask[index] OR sdss_flagval('MANGA_DRP2PIXMASK','3DREJECT')
  endif

  ; Figure out where the broken fibers were using the slitmap.gbu flag
  ; and add this information to the MANGA_DRP2PIXMASK bitmask
  index=where(partialslit.gbu ne 0,nindex)
  if (nindex ne 0) then begin
    tempmask=replicate(0L,nfiber)
    tempmask[index]=sdss_flagval('MANGA_DRP2PIXMASK','DEADFIBER')
    newmask = newmask OR rebin(reform(tempmask,1,nfiber),nlinwave,nfiber)
  endif

  ; Look for where the inverse variance is zero, and zero out the
  ; flux at these locations too so that they don't confuse anyone
  index=where(objivar eq 0.,nindex)
  if (nindex gt 0) then objflux[index]=0.

  ; Stuff the flux, inverse variance, and new bitmask into the big arrays
  flux_arr[*,fiberstart:fiberstart+nfiber-1]=objflux[*,*]
  ivar_arr[*,fiberstart:fiberstart+nfiber-1]=objivar[*,*]
  mask_arr[*,fiberstart:fiberstart+nfiber-1]=newmask[*,*]
  disp_arr[*,fiberstart:fiberstart+nfiber-1]=dispimg[*,*]

  ; Increment the fiber running index
  fiberstart+=nfiber
endfor

; Derive x_linarr and y_linarr from the astrometry already derived
; for x_arr and y_arr by simple interpolation
for i=0,nftotal-1 do begin
  x_linarr[*,i]=interpol(x_arr[*,i],wave,linwave)
  y_linarr[*,i]=interpol(y_arr[*,i],wave,linwave)
endfor


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Create data cube in linear wavelength grid
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Loop over all wavelengths figuring out the cube at that wavelength
; (unless we're in bright-time mode, in which case don't do this loop;
; cubes will be just empty placeholders.  This saves ENORMOUS cpu time.)
if (srvymode ne 'APOGEE') then begin
  for i=0,cube_zsize-1 do begin
    ; Print a message to splog every 10% of the loop
    if (i mod fix(cube_zsize/10) eq 0) then $
      splog,'Constructing cube: '+strcompress(string(round(i*100./cube_zsize)),/remove_all)+'% complete'

    ; Create input vectors for this wavelength slice
    xvec[*]=x_linarr[i,*]/cubescale+cube_xsize/2.
    yvec[*]=y_linarr[i,*]/cubescale+cube_ysize/2.
    fvec[*]=flux_arr[i,*]
    ivec[*]=ivar_arr[i,*]
    mvec[*]=mask_arr[i,*]

    ; Scale correction factor is the ratio of area between a PI area fiber
    ; (in arcsec^2) and the output pixel size in arcsec^2
    ; The result means that the cube will be in calibrated units/pixel
    scale=cubescale*cubescale/!PI

    cube_flux[*,*,i]=ml_griddata(xvec,yvec,fvec,[cube_xsize,cube_ysize],1.6/cubescale,0.7/cubescale,ivar=ivec,invarimg=invarimg,maskvec=mvec,maskimg=maskimg,scale=scale)
    ; Populate inverse variance and mask arrays
    cube_ivar[*,*,i]=invarimg
    cube_mask[*,*,i]=maskimg
  endfor
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Work out the median spectral resolving power
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

specreslin=interpol(specres,wave,linwave)
specreslin_stddev=interpol(specres_stddev,wave,linwave)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Make the basic cube and RSS file headers
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

mkhdr,rsshdr,flux_arr,/image
mkhdr,cubehdr,cube_flux,/image
; Add MaNGA keycards
mdrp_drp3dhead,drp3qual=drp3qual,cubehdr=cubehdr,rsshdr=rsshdr,obsinfo=obsinfo,mgcf_hdr=hdrsave,wave=linwave,cubescale=cubescale,cubefwhm=cubefwhm,cube0hdr=cube0hdr,rss0hdr=rss0hdr

; Make headers for the ivar and mask extensions
mkhdr,rss_ivar_hdr,ivar_arr,/image
mkhdr,rss_mask_hdr,mask_arr,/image
mkhdr,cube_ivar_hdr,cube_ivar,/image
mkhdr,cube_mask_hdr,cube_mask,/image

; Add info to the global headers to indicate flux, error, mask extensions
mdrp_hduclass,cubehdr,errhdr=cube_ivar_hdr,maskhdr=cube_mask_hdr,type='CUBE'
mdrp_hduclass,rsshdr,errhdr=rss_ivar_hdr,maskhdr=rss_mask_hdr,type='IMAGE'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Quality check- look for known foreground stars in the data cube
; so that we can flag them
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

maskradius=2.5; arcsec
thismangaid=strcompress(fxpar(cubehdr,'MANGAID'),/remove_all)
knownstars=yanny_readone(concat_dir(ml_getenv('MANGACORE_DIR'),'platedesign/foregroundstars/foregroundstars.par'),'FORESTAR')
index=where(strcompress(knownstars.mangaid,/remove_all) eq thismangaid,nindex)
starmask=intarr(cube_xsize,cube_ysize)

for k=0,nindex-1 do begin
  maskra=knownstars[index[k]].raj2000
  maskdec=knownstars[index[k]].dej2000
  starmask[*,*]=0L
  for i=0,cube_xsize-1 do begin
    for j=0,cube_ysize-1 do begin
      xyad,cubehdr,i,j,thisra,thisdec
      deltara=3600.*(maskra-thisra)*cos(thisdec*!PI/180.)
      deltadec=3600.*(maskdec-thisdec)
      if (sqrt((deltara*deltara)+(deltadec*deltadec)) lt maskradius) then starmask[i,j]=1
    endfor
  endfor
  starmask=starmask*sdss_flagval('MANGA_DRP3PIXMASK','FORESTAR')
  for q=0,cube_zsize-1 do cube_mask[*,*,q]=cube_mask[*,*,q] or starmask
endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Write out LIN data files
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Write out the RSS-format file
rssname=outdir+outnameroot+'LINRSS.fits'
ml_mwrfits, dummyext, rssname, hdr=rss0hdr, /create
ml_mwrfits, flux_arr, rssname, hdr=rsshdr, extname='FLUX'
ml_mwrfits, ivar_arr, rssname, hdr=rss_ivar_hdr, extname='IVAR'
ml_mwrfits, mask_arr, rssname, hdr=rss_mask_hdr,extname='MASK'
ml_mwrfits, disp_arr, rssname, extname='DISP'
ml_mwrfits, linwave, rssname, extname='WAVE'
ml_mwrfits, specreslin, rssname, extname='SPECRES'
ml_mwrfits, specreslin_stddev, rssname, extname='SPECRESD'
ml_mwrfits, obsinfo, rssname, extname='OBSINFO'
ml_mwrfits, x_linarr, rssname, extname='XPOS'
ml_mwrfits, y_linarr, rssname, extname='YPOS'

; Write out the cube
cubename=outdir+outnameroot+'LINCUBE.fits'
ml_mwrfits, dummyext, cubename, hdr=cube0hdr, /create
ml_mwrfits, cube_flux, cubename, hdr=cubehdr, extname='FLUX'
ml_mwrfits, cube_ivar, cubename, hdr=cube_ivar_hdr, extname='IVAR'
ml_mwrfits, cube_mask, cubename, hdr=cube_mask_hdr, extname='MASK'
ml_mwrfits, linwave, cubename, extname='WAVE'
ml_mwrfits, specreslin, cubename, extname='SPECRES'
ml_mwrfits, specreslin_stddev, cubename, extname='SPECRESD'
ml_mwrfits, obsinfo, cubename, extname='OBSINFO'
ml_mwrfits, g_img, cubename, hdr=g_hdr, extname='GIMG'
ml_mwrfits, r_img, cubename, hdr=r_hdr, extname='RIMG'
ml_mwrfits, i_img, cubename, hdr=i_hdr, extname='IIMG'
ml_mwrfits, z_img, cubename, hdr=z_hdr, extname='ZIMG'
ml_mwrfits, psfcube[*,*,0], cubename, hdr=g_hdr, extname='GPSF'
ml_mwrfits, psfcube[*,*,1], cubename, hdr=r_hdr, extname='RPSF'
ml_mwrfits, psfcube[*,*,2], cubename, hdr=i_hdr, extname='IPSF'
ml_mwrfits, psfcube[*,*,3], cubename, hdr=z_hdr, extname='ZPSF'

; gzip output files
spawn, ['gzip', '-f', rssname], /noshell
spawn, ['gzip', '-f', cubename], /noshell

; End the ELSE criterion that required no missing LINEAR format files
endelse

return
end