;+
;
; Project   : s3drs reconstruction software
;                   
; Name      : (various test routines)
;               
; Purpose   : part of the s3drs test kit
;               
; Explanation: (archived here for completeness, not intended for users use)
;               
; Use       : (self-documented, see notes in code)
;    
; Written     : Sandy Antunes, NRC, 2005-2009
;               
;-            
; now allows optional insert of FM fluxrope, e.g.
; e.g. restore,/v,'18_imgp.sav'
;      supercone,N=64,fm=sol_image,/show,prob_cutoff=0
;

PRO supercone,N=N,fm=fm,cone=cone,show=show,prob_cutoff=prob_cutoff,$
              occulternix=occulternix,pure=pure
  
  ; ** iterative method.

  cone=setyesno(cone)
  if (n_elements(N) eq 0) then N=64
  show=setyesno(show)
  occulternix=setyesno(occulternix); removes xyzmask false occulter region
  pure=setyesno(pure); just reruns solver
  
  pixon_env,'c'

  infiles=deriveseq('20071231')

  for ifile=0,n_elements(infiles)-1 do $
    infiles[ifile]=infiles[ifile]+'_20071231_002220_'+int2str(N)+'_zmask_solved.sav'

  
  for ifile=0,n_elements(infiles)-1 do begin

    ; 1) Start with the current Pixon best solution.
    
;    fname='../dec31seq64/' + infiles[ifile]
      fname= infiles[ifile]
      print,'checking ',fname
      if (file_exist(fname)) then begin

          restore,/v,fname
          
         ; 2) make sure you find the z-element, usually # 2
          print,pxwrapper.sr.rho,pxwrapper.sr.phi
          zele=2
          
         ; 3) set iterations to something reasonable
          pxwrapper.sp.mxit=100
          opixon->set,mxit=pxwrapper.sp.mxit
          
         ; 4) generate a probability map for the z-axis view
          
          data=stereo_fetch(pxwrapper,opixon,/data)

         ; note choice of 'prob space' or an input FM to determine Zmask
          if (n_elements(fm) ne 0) then sol=fm else $
            sol=simple3drv(data[*,*,0:1],pxwrapper=pxwrapper,/usemax)

          fakez=pr_render(pxwrapper,sol,zele)
          fakez=fakez/max(fakez) ; normalize as 0-1

          if (show) then tv_multi,fakez,res=2,/log,title='fakez'

          if (n_elements(prob_cutoff) eq 0) then prob_cutoff=0.0
;            if (n_elements(fm) eq 0) then prob_cutoff=0 else prob_cutoff=0.1
          
         ; 5) extract the rest of our original input items, for later reinsertion
          
          sigma=stereo_fetch(pxwrapper,opixon,/sigma)
          masko=p_getmask(sr=pxwrapper.sr)
          if (occulternix) then mask=p_getmask(sr=pxwrapper.sr,rocc=0) else $
            mask=masko

         ; test-- just insert as-is
         ;p_insert_data,pxwrapper,data,opixon,sigma=sigma,masks=masks
         ;opixon->full
          
          
         ; now make z-mask of '1' outside of region of interest with
         ; corresponding data of '0' everywhere
;    zdat=data[*,*,zele]; should be filled with uniform 'dmin' value
          
          zmask=~mask ; currently 1 in occulter and edge, i.e. enforces zeroness

          if (pure eq 0 and cone eq 0) then begin
              ie=where (fakez le prob_cutoff)
              ;also set to 1 in low-prob region to enforce zeroness
              if (ie[0] ne -1) then zmask[ie]=1
          end else if (cone) then begin
              conemask=mask * 0.0 + 1.0
              conemask[N/2:*,0:N/2]=0
              conemask=rot(conemask,35)
              ic=where(conemask eq 1)
              zmask[ic]=1
          end
          
          
          
          datz=data[*,*,zele]
          sigz=datz*10.0        ; wildly noise
          sigma[*,*,zele]=sigz
          
 
; hack for 025220
;maskA=mask         
;maskA[0:8,11:25]=0
          
          if (pure) then zmask=~mask; ignores above, just for rerunning/testing
          
         ; only element we need to alter is the z-mask
          masks=[[[masko]],[[masko]],[[zmask]],[[~mask]],[[~mask]]]
          
         if (show) then tv_multi,masks,title='Input Masks'

         if (show) then tv_multi,masks*data,title='Masked Input Data',/log,oom=2

          p_insert_data,pxwrapper,data,opixon,sigma=sigma,masks=masks
         ;   7c) Run!
          opixon->set,mxit=400
          opixon->full
         ;p_standardplots,pxwrapper,opixon,/best,/quick
          
          
          print,'saving ',fname+'.sav'
          save,opixon,pxwrapper,file=fname+'.sav'
          
         ; different discarded approaches here
          
         ; 7) now update z-data to enforce region, iteratively, by constantly
         ;    weighting the z-solution by the normalized region mask in order
         ;    to nudge z to spec.
          
         ;   7a) make fake Z data, and then generate an approximate noise
;    render=stereo_fetch(pxwrapper,opixon,/render)
;    datz=render[*,*,zele]
;    sigz=2.7E-12 * sqrt( ratio(datz,2.7E-12) )
         ; now jam down the signal in the parts with poor s/n to constraint
         ;ie=where (ratio(datz,sigz) le 10)
         ;   datz[ie]=0.0  zmask=mask, zmask[ie]=0
         ;if (ie[0] ne -1) then sigz[ie]=datz[ie]
          
          
;    zmask=mask
;    ie=where (fakez lt 0.1)
;    if (ie[0] ne -1) then zmask[ie]=0.0
;    if (ie[0] ne -1) then sigz[ie]=datz[ie]*0.1
;    if (ie[0] ne -1) then datz[ie]=0.0
         ;zmask=zmask*fakez; do not take Z as seriously as the others
          
;    masks=[[[mask]],[[mask]],[[zmask]],[[~mask]],[[~mask]]]
          
;    data[*,*,zele]=datz*5.0; artifically jam up z signal
;    sigz=2.7E-12 * sqrt( ratio(datz,2.7E-12) ); recompute sigma
;    sigma[*,*,zele]=sigz
          
         ; tv_multi,fixed=196,ratio(data,sigma) > 1 < 10
         ;   7b) insert
;    p_insert_data,pxwrapper,data,opixon,sigma=sigma,masks=masks
         ;   7c) Run!
;    opixon->full
         ;p_standardplots,pxwrapper,opixon,/best,/quick
          
      end
  end
  
END