; this routine takes the existing forward model and uses it to
; tightly constrain the area of the inverse reconstruction.
;
; Assumes file 'fmfile' contains prior solution as a 3D datacube.
;   (It can have any name, but it must be the only 3D item in
;    the file for us to automatically find it).
;

PRO envelope_reconstruction,N=N,zmask=zmask,fmfile=fmfile,$
                            cmedate=cmedate,priordate=priordate

  if (n_elements(N) eq 0) then N=128

  zmask=setyesno(zmask)         ; default is no, use xyzmask
  xyzmask=~zmask                ; either one or the other

  if (n_elements(fmfile) eq 0) then fmfile='best_fluxrope_img.sav'
  if (n_elements(cmedate) eq 0) then cmedate='20071231_022220'
  if (n_elements(priordate) eq 0) then priordate='20071231_012220'
  
; bring in the forward envelope by reading any 3D element from file
  fmobj=obj_new('IDL_Savefile',fmfile)
  names=fmobj->names()
  idata=0
  for i=0,n_elements(names)-1 do $
    if ( (fmobj->size(names[i]))[0] eq 3) then idata=i
  fmobj->restore,names[idata]
  returncode=execute('sol_image='+names[idata])
  obj_destroy,fmobj

  restore,/v,fmfile ; contains sol_image
  sol_image=rebin(sol_image,N,N,N)
  
; make the surrounding structure, including a zmask for later replacement
  quickcme,N,cmedate,priordate,/cor2,zmask=zmask,xyzmask=xyzmask,mxit=2
  str='zmask'
  if (xyzmask) then str='xyzmask'
  fname=cmedate+'_'+priordate+'_'+int2str(N)+'_cor2_'+str+'_solved.sav'
  restore,/v,fname
  
  zele=2
  prob_cutoff=0.0
  
;extract z mask element
  masks=stereo_fetch(pxwrapper,oPixon,/mask)
  zmask=masks[*,*,zele:*]
  masks=p_getmask(sr=pxwrapper.sr)
; create new z-element data and mask
  
  faked=pr_render(pxwrapper,sol_image)
  fakez=faked[*,*,zele:*]
  fakez=fakez/max(fakez)        ; normalize as 0-1
  
  zmask=fakez*0.0
  ie=where (fakez le prob_cutoff)
  if (ie[0] ne -1) then zmask[ie]=1 ; enforce zeroness
  
  masks[*,*,zele:*]=zmask
  
  data=stereo_fetch(pxwrapper,oPixon,/data)
  sigma=stereo_fetch(pxwrapper,oPixon,/sigma)
  sigz=sigma[*,*,zele:*]
  sigz=data[*,*,zele:*]*10.0
  sigma[*,*,zele:*]=sigz
  
; replace
  p_insert_data,pxwrapper,data,opixon,sigma=sigma,masks=masks
  
; check
  tv_multi,/log,/local,/show,data,header='data'
  tv_multi,/log,/local,/show,sigma,header='sigma'
  tv_multi,/log,/local,/show,masks,header='masks'
  tv_multi,/log,/local,/show,stereo_fetch(pxwrapper,oPixon,/wt),header='wt'
  
; run!
  pxwrapper.sp.mxit=200
  opixon->set,mxit=pxwrapper.sp.mxit
  oPixon->full,/restart
  
  save,pxwrapper,opixon,file='envelope_solved_'+str+'_'+int2str(N)+'.sav'
  
  if (xyzmask) then line=5 else line=3
  im=stereo_fetch(pxwrapper,opixon,/im)
  
; cut out non-fluxrope entirely from cube, to remove artifacts
  cie=where(sol_image le 0)
  if (cie[0] ne -1) then im[cie]=0
  
  render=pr_render(pxwrapper,im)
  tv_multi,/log,/local,[[[[render]]],[[data]]],header='data,render',line=line
  tv_multi,/log,/local,[[[[render*masks]]],[[data*masks]]],line=line,$
    header='data,render (masked)'
  threeview,/log,/pixon,im

END