;+
;
; 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
;               
;-            
; sample:
;quickcme,64,'20071231_025220','20071231_022220',/c3
;
; use of /c3 means it will look for the c3 file closest to the given
; times.  Ditto for /c2

pro quickCME,n,cmedate,priordate,zmask=zmask,xyzmask=xyzmask,$
             checksn=checksn,dn=dn,excerpt=excerpt,c2=c2,c3=c3,cor1=cor1,$
             sidemask=sidemask,nozero=nozero,$
             mxit=mxit,occulttop=occulttop,occultbottom=occultbottom,$
             occultleft=occultleft,occultright=occultright,pseudo=pseudo,$
             maskratio=maskratio,nadjust=nadjust,setuponly=setuponly,$
             brightpixels=brightpixels,sncutoff=sncutoff,$
             standardplots=standardplots,flatnoise=flatnoise,best=best

  on_error,2
  pixon_env,'c',/local; ensure we are using local, latest version

  msbperdn=2.7E-12

  standardplots=setyesno(standardplots)
  flatnoise=setyesno(flatnoise)
  c2=setyesno(c2)
  c3=setyesno(c3)
  cor1=setyesno(cor1)
  
  checksn=setyesno(checksn)
  setuponly=setyesno(setuponly)
  dn=setyesno(dn)
  pseudo=setyesno(pseudo)
  brightpixels=setyesno(brightpixels)
  zmask=setyesno(zmask)  
  xyzmask=setyesno(xyzmask)  
  if (xyzmask) then zmask=0; redundant
  sidemask=setyesno(sidemask)
  nozero=setyesno(nozero)
  excerpt=setyesno(excerpt)
  if (n_elements(sncutoff) eq 0) then sncutoff=-1
  occultleft=setyesno(occultleft)
  occultright=setyesno(occultright)
  occulttop=setyesno(occulttop)
  occultbottom=setyesno(occultbottom)

  best=setyesno(best)
  if (best) then begin
      nozero=1
      xyzmask=1
  end

  while (n_elements(N) eq 0) do begin
      print,'Input a value for N of 64, 128, or 256 (256 is risky)'
      read,'N = ? ',Nval
      if (Nval eq 64 or Nval eq 128 or Nval eq 256) then N=floor(Nval)
  end
  
  if (n_elements(cmedate) eq 0) then begin
      cmedate=""
      print,'Enter the datestamp of the CME, e.g. 20071104_222220'
      read,cmedate
  end

;   Nd=2
;   if (lasco) then ++Nd
;   if (zmask) then ++Nd
;   if (xyzmask) then Nd=Nd+3
  
  ; cmedate='20071104_212220'
  ;cmedate='20071104_222220'
  ; cmedate='20071104_232220'

  cpieces=strsplit(cmedate,'_',/extract)

  if (n_elements(priordate) eq 0) then begin  
      priordate=""
      print,'To run as a difference image, provide the datestamp of the'
      print,' immediately prior frame (for a running difference, e.g.'
      print,' 20071104_212220) or of an earlier quiet frame (for a'
      print,' difference image, e.g. 20071104_022220).  If you wish to'
      print,' instead just subtract a standard background, type "background"'
      read,priordate
  
  end

  if (strmid(priordate,0,4) eq 'back') then priordate='back'

  Nfactor=2048/float(N)         ; e.g. 1/16 for N=128 so Nfactor=16
  
  print,'Setting up for ',cmedate,' - ',priordate

  ; ** SECCHI FILENAMES
  
  ;quietdate='20071104_022220'; using an early frame
  ;runningdate='20071104_212220'; a running difference
  
;  tbdira='/net/pluto/p1/secchi/P0/a/cor2/'
;  tbdirb='/net/pluto/p1/secchi/P0/b/cor2/'
  tbdira=getenv('SECCHI_P0')+'/a/cor2/'
  tbdirb=getenv('SECCHI_P0')+'/b/cor2/'
  
  cmeA= tbdira + cpieces[0] + '/' + cmedate + '_0B4c2A.fts'
  cmeB= tbdirb + cpieces[0] + '/' + cmedate + '_0B4c2B.fts'

  if (file_exist(cmeA) eq 0 or file_exist(cmeB) eq 0) then begin
      print,'Missing input file CME, exiting'
      return
  end
  
  dcA=rebin(sccreadfits(cmeA,hdrA),N,N)
  dcB=rebin(sccreadfits(cmeB,hdrB),N,N)*2.8/2.7
  
  hdrA.exptime=1                ; required for scc_getbkgimg
  hdrB.exptime=1                ; required for scc_getbkgimg
  
  
  if (priordate eq 'back') then begin

      dpA=rebin(scc_getbkgimg(hdrA,/match,/totalb),N,N)
      dpB=rebin(scc_getbkgimg(hdrB,/match,/totalb),N,N)*2.8/2.7

  end else begin

      ppieces=strsplit(priordate,'_',/extract)
      priorA= tbdira + ppieces[0] + '/' + priordate + '_0B4c2A.fts'
      priorB= tbdirb + ppieces[0] + '/' + priordate + '_0B4c2B.fts'
  

      if (file_exist(priorA) eq 0 or file_exist(priorB) eq 0) then begin
          print,'Missing input prior file, exiting'
          return
      end

      dpA=rebin(sccreadfits(priorA),N,N)
      dpB=rebin(sccreadfits(priorB),N,N)*2.8/2.7
  
  end

  rawdata=[[[dcA]],[[dcB]]]
  priordata=[[[dpA]],[[dpB]]]
  
  diffA=dcA-dpA
  diffB=dcB-dpB

  dtemp=[[[diffA]],[[diffB]]]
  ie=where(dtemp gt 0)
  if (ie[0] eq -1) then dmin=0.0 else dmin=min(dtemp[ie])

  ; ** LASCO FILENAMES
  
  if (c2 or c3) then begin

      lasco1=find_c2c3(time1)
      lasco2=find_c2c3(time2)
;      lascodir='/home/antunes/RUNS/lasco/c3/'
      cmeL=lasco1
      priorL=lasco2
;      p_ingest_fits,cmeL
      print,cmeL,priorL
      lasco_prep,cmeL,hdrL,dcLfull
      lasco_prep,priorL,hdrLignore,dpLfull

      if (c2) then begin
          calfac=c2_calfactor(hdrL)
      end else begin
          calfac=c3_calfactor(hdrL)
      end
      propalter=calfac/msbperdn ; normalize DN for later bulk convert

      if (excerpt) then begin
      ; subselect item
      ; assume C3, aka FOV=2*Cor2
      ; so 0-1023 becomes 0-511
          No=1024
          Np=512
          dcLhalf=dblarr(Np,Np)
          dcLhalf[0:Np-1,0:Np-1]=dcLfull[No/4:(3*No/4)-1,No/4:(3*No/4)-1]
          dpLhalf=dblarr(Np,Np)
          dpLhalf[0:Np-1,0:Np-1]=dpLfull[No/4:(3*No/4)-1,No/4:(3*No/4)-1]

          dcL=rebin(dcLhalf,N,N) * propalter
          dpL=rebin(dpLhalf,N,N) * propalter

      end else begin
          ; full frame

          dcL=rebin(dcLfull,N,N) * propalter
          dpL=rebin(dpLfull,N,N) * propalter

      end

;      dcL=rebin(sccreadfits(cmeL,hdrL),N,N)
;      dpL=rebin(sccreadfits(priorL),N,N)

      ; handle quicklook artifacts
      missingL=where(dcL eq 0)
      missingLb=where(dpL eq 0)
      if (missingLb[0] ne -1) then begin
          if (missingL[0] eq -1) then missingL=missingLb else $
            missingL=[missingL,missingLb]
      end

      diffL=dcL-dpL

  end                           ; lasco

  ; ** DATA IN
  
  ie= where (dcA lt 0)
  if (ie[0] ne -1) then dcA[ie]=0.0
  ie= where (dcB lt 0)
  if (ie[0] ne -1) then dcB[ie]=0.0
  ie= where (dpA lt 0)
  if (ie[0] ne -1) then dpA[ie]=0.0
  ie= where (dpB lt 0)
  if (ie[0] ne -1) then dpB[ie]=0.0
  
  ; ** FITS INGEST
  
  gamma=[0.0,0.0]           ; assume 2 initial views
  theta=[0.0,0.0]          ; assume 2 initial views
  if (zmask or xyzmask) then begin
      gamma=[gamma,0.0]; add Z
      theta=[theta,90.0]; add Z
      if (xyzmask) then gamma=[gamma,90.0,0.0]; also add XY
      if (xyzmask) then theta=[theta,0.0,0.0]; also add XY
  end
  if (c2) then theta=[0.0,theta] ; front-add space for LASCO, too
  if (c3) then theta=[0.0,theta] ; front-add space for LASCO, too
  if (cor1) then theta=[0.0,0.0,theta] ; front-add space for Cor1, too
  if (c2) then gamma=[0.0,gamma] ; front-add space for LASCO, too
  if (c3) then gamma=[0.0,gamma] ; front-add space for LASCO, too
  if (cor1) then gamma=[0.0,0.0,gamma] ; front-add space for Cor1, too
  prep_3drv,pxwrapper,N=N,theta=theta,gamma=gamma,/cor2
  
  if (dn) then pxwrapper.sr.imodel=5; for input in DN

  ; store the specifics of this setup, for record-keeping
  if (xyzmask) then pxwrapper.sc.xyzmask=1
  if (zmask) then pxwrapper.sc.zmask=1
  if (sidemask) then pxwrapper.sc.sidemask=1
  if (nozero) then pxwrapper.sc.nozero=1
  if (excerpt) then pxwrapper.sc.excerpt=1
  if (occultleft or occultright or occulttop or occultbottom) then $
    pxwrapper.sc.occult=1
  if (brightpixels) then pxwrapper.sc.brightpixels=1
  if (n_elements(nadjust) ne 0) then pxwrapper.sc.nadjust=nadjust

  pxwrapper.sc.gname=cmedate+'_'+int2str(N)+'_'+priordate
  if (c2) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_c2'
  if (c3) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_c3'
  if (cor1) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_cor1'
  if (zmask) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_zmask'
  if (xyzmask) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_xyzmask'
  if (sidemask) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_sidemask'
  if (nozero) then pxwrapper.sc.gname = pxwrapper.sc.gname + '_nozero'
  
  p_ingest_fits,cmeA,ignore2,maskA,N=N,pxwrapper=pxwrapper,iele=0
  p_ingest_fits,cmeB,ignore2,maskB,N=N,pxwrapper=pxwrapper,iele=1
  if (c2) then $
    p_ingest_fits,cmeL,ignore2,maskc2,N=N,pxwrapper=pxwrapper,iele=2
  if (c2 and c3) then $
    p_ingest_fits,cmeL,ignore2,maskL,N=N,pxwrapper=pxwrapper,iele=3

  if (n_elements(maskratio) ne 0) then begin
      maskA=p_getmask(sr=pxwrapper.sr,iele=0,ratio=maskratio)
      maskB=p_getmask(sr=pxwrapper.sr,iele=1,ratio=maskratio)
      if (lasco) then maskL=p_getmask(sr=pxwrapper.sr,iele=2,ratio=maskratio)
  end

  if (excerpt) then begin
      ; subsection alteration
      pxwrapper.sr.l0tau[2]=pxwrapper.sr.l0tau[2]*0.5
      pxwrapper.sr.tau[2]=pxwrapper.sr.tau[2]*0.5
      ; dvoxfrac, flibm and istd will be recomputed in p_osetup later
  end

  ; compensate for SECCHI frames that are actually 2x3=6sec
  if (long(pxwrapper.sr.exptime[0]) eq 2) then pxwrapper.sr.exptime[0]=6.0
  if (long(pxwrapper.sr.exptime[1]) eq 2) then pxwrapper.sr.exptime[1]=6.0



  ; ** NOISE
  ;  noise(DN_sec)=noise(DN_tot)/6sec=noise(p*15)/6sec
  ;  noise(p) = sqrt(p) = sqrt(DN_tot/15) = sqrt(DN_sec*6sec/15)
  ; so noise(DN_sec)=sqrt(DN_sec*6/15)/6sec = sqrt(6/15)/6.0 * sqrt(DN_sec)
  ; so noise(DN_sec)=0.105409 * sqrt(DN_sec)
  ; Also noise total is sqrt of sum of 2 frames, roughly 1.4 each frame:
  ;   noise(DN_sec_diff) ~= 0.15 sqrt(DN_sec_raw)
  ;
  ; Since DN_diff = Q * DN_raw and noise(DN_sec_cme)=0.15 sqrt(DN_sec_raw)
  ; then  S/N = DN_sec_diff/noise(DN_sec_raw)=Q*DN_sec_raw/0.15 sqrt(DN_sec_raw)
  ;    so S/N = (Q/0.15) sqrt(DN_sec_raw),
  ; so S/N > 1 if sqrt(DN_sec_raw > (0.15/Q)
  
;;  const=1.0 / (sqrt(90.0) * Nfactor)
;;  
;;  ncA=const * sqrt(dcA)
;;  ncB=const * sqrt(dcB)
;;  if (lasco) then ncL=const * sqrt(dcL)
;;  
;;  npA=const * sqrt(dpA)
;;  npB=const * sqrt(dpB)
;;  if (lasco) then npL=const * sqrt(dpL)
;;  
;;  ndiffA=sqrt(ncA^2 + npA^2)
;;  ndiffB=sqrt(ncB^2 + npB^2)
;;  if (lasco) then ndiffL=sqrt(ncL^2 + npL^2)
  
;  nadjust=15.0
  if (n_elements(nadjust) eq 0) then nadjust=1.0

  if (flatnoise) then begin
      ; units are currently in DN, so noise is
      ; sqrt(DN/15.0) * 15.0 * (N/2048.) * 1.0/sqrt(exp) for 2 frames
      ; = sqrt(DN) * N * (15/ [sqrt(15) * 2048 * sqr(6.0)] ) * sqrt(2)
      ; = 0.00109183 * N * sqrt(DN)
      ndiffA=sqrt(dcA) * 0.00109183 * N/nadjust
      ndiffB=sqrt(dcB) * 0.00109183 * N/nadjust
      if (lasco) then ndiffL=sqrt(dcL) * 0.00109183 * N/nadjust
  end else begin
      ; better noise treatment
      ndiffA=p_quicknoise(dcA,dpA,pxwrapper=pxwrapper,/ingestdn,iele=0)/nadjust
      ndiffB=p_quicknoise(dcB,dpB,pxwrapper=pxwrapper,/ingestdn,iele=1)/nadjust
      if (lasco) then $
        ndiffL=p_quicknoise(dcL,dpL,pxwrapper=pxwrapper,/ingestdn,iele=2)$
        /nadjust
  end

; SN CUTOFF/FEATURES ONLY mod
; mod to extract only features and zero out data where SN<cutoff
  if (sncutoff gt 0) then begin
      isn=where (ratio(diffA,ndiffA) le sncutoff)
      if (isn[0] ne -1) then diffA[isn]=0.0
      isn=where (ratio(diffB,ndiffB) le sncutoff)
      if (isn[0] ne -1) then diffB[isn]=0.0
      if (lasco) then begin
          isn=where (ratio(diffL,ndiffL) le sncutoff)
          if (isn[0] ne -1) then diffL[isn]=0.0
      end
  end
  
  

; *** MASKING

  ; for zmasking, the data is the mask, mask=inv data
  if (zmask or xyzmask) then $
    blank=p_getmask(sr=pxwrapper.sr,iele=2,rocc=-1)
  if (zmask or xyzmask) then datumZ=blank*0.0 + dmin*0.1; low everywhere
  if (zmask or xyzmask) then maskZ= ~blank
  if (zmask or xyzmask) then sigmaZ= datumZ * 10; wildly noisy
;  if (zmask or xyzmask) then print,'zmask min and max are: ',minmax(datumZ)


  ; extremely obscure option, cuts off right side of both A & B with mask

  if (occultright) then begin
      maskA[N/2:*,*]=0.0
      maskB[N/2:*,*]=0.0
      if (lasco) then maskL[N/2:*,*]=0.0
  end
  if (occultleft) then begin
      maskA[0:N/2,*]=0.0
      maskB[0:N/2,*]=0.0
      if (lasco) then maskL[0:N/2,*]=0.0
  end
  if (occulttop) then begin
      maskA[*,N/2:*]=0.0
      maskB[*,N/2:*]=0.0
      if (lasco) then maskL[*,N/2:*]=0.0
  end
  if (occultbottom) then begin
      maskA[*,0:N/2]=0.0
      maskB[*,0:N/2]=0.0
      if (lasco) then maskL[*,0:N/2]=0.0
  end

  ; primitive hot pixel remover, adds to mask and doesn't change data
  if (brightpixels) then begin

      dcapA=percentiles(diffA,value=[0.99])
      ibrightA=where(diffA gt dcapA)
      help,ibrightA
      for i=0,n_elements(ibrightA)-1 do begin
          myi=ibrightA[i]
          if (myi ne 0 and myi ne N-1) then $
            if (diffA[myi] gt diffA[myi-1] + diffA[myi+1]) then maskA[myi]=0
      end

      dcapB=percentiles(diffB,value=[0.99])
      ibrightB=where(diffB gt dcapB)
      for i=0,n_elements(ibrightB)-1 do begin
          myi=ibrightB[i]
          if (myi ne 0 and myi ne N-1) then $
            if (diffB[myi] gt diffB[myi-1] + diffB[myi+1]) then maskB[i]=0
      end

  end


  ; to add a sigma and value to side masks is tricky
  ; mask is 1 for good data, 0 for bad
  ; so we want to:
  ;   1) where mask=0, set data=dmin*0.1
  ;   2) where mask=0, set sigma=dmin*10
  ;   3) set mask=1 everywhere
  if (sidemask) then begin
      ieA=where(maskA eq 0)
      if (ieA[0] ne -1) then diffA[ieA]=dmin*0.1
      if (ieA[0] ne -1) then ndiffA[ieA]=dmin*10.0
      maskA=maskA*0.0+1.0

      ieB=where(maskB eq 0)
      if (ieB[0] ne -1) then diffB[ieB]=dmin*0.1
      if (ieB[0] ne -1) then ndiffB[ieB]=dmin*10.0
      maskB=maskB*0.0+1.0

      if (lasco) then ieL=where(maskL eq 0) else ieL=[-1]
      if (ieL[0] ne -1) then diffL[ieL]=dmin*0.1
      if (ieL[0] ne -1) then ndiffL[ieL]=dmin*10.0
      if (lasco) then maskL=maskL*0.0+1.0
      if (lasco) then if (missingL[0] ne -1) then maskL[missingL]=0.0

  end else begin

      if (lasco) then if (missingL[0] ne -1) then maskL[missingL]=0.0

  end

;  tv_multi,res=3,[[[maskA]],[[diffA]],[[ndiffA]]],title='one',/local,/log

  if (c2 or c3) then begin
      ; change solar north to ecliptic north
      DEFSYSV, "!DRADEG", 180.0d/!DPI
      DEFSYSV, "!DDTOR", !DPI/180.0d
      xyz=*(pxwrapper.sr.hae_coords[2])
      angle = festival_ecl_solar_north(xyz[0], xyz[1], xyz[2])
      xy=floor(p_getk0(pxwrapper.sr,2))
      diffL=rot(diffL,angle,1.0,xy[0],xy[1])
      ndiffL=rot(ndiffL,angle,1.0,xy[0],xy[1])
      maskL=rot(maskL,angle,1.0,xy[0],xy[1])
  end

  
  ; *** false zero mod, to remove actual data zeroes, now moved to later
;  if (nozero) then begin
;      ieA=where(diffA le 0)
;      if (ieA[0] ne -1) then diffA[ieA]=dmin*0.1
;      if (ieA[0] ne -1) then ndiffA[ieA]=dmin*10.0

;      ieB=where(diffB le 0)
;      if (ieB[0] ne -1) then diffB[ieB]=dmin*0.1
;      if (ieB[0] ne -1) then ndiffB[ieB]=dmin*10.0

;      if (lasco) then ieL=where(diffL le 0) else ieL=[-1]
;      if (ieL[0] ne -1) then diffL[ieL]=dmin*0.1
;      if (ieL[0] ne -1) then ndiffL[ieL]=dmin*10.0
;  end


  ; *** Now sum everything together

  diff=[[[diffA]],[[diffB]]]
  if (c2) then diff=[[[diff]],[[diffc2]]]
  if (c3) then diff=[[[diff]],[[diffc3]]]
  if (cor1) then diff=[[[diff]],[[diffcor1]]]

  ndiff=[[[ndiffA]],[[ndiffB]]]
  if (c2) then ndiff=[[[ndiff]],[[ndiffc2]]]
  if (c3) then ndiff=[[[ndiff]],[[ndiffc3]]]
  if (cor1) then ndiff=[[[ndiff]],[[ndiffcor1]]]

  ; *** false zero mod, to remove actual data zeroes
  if (nozero) then p_nozero,diff,ndiff

  if (lasco) then mset=[[[maskA]],[[maskB]],[[maskL]]] else $
    mset=[[[maskA]],[[maskB]]]

  if (zmask) then mset=[[[mset]],[[maskZ]]]
  if (xyzmask) then mset=[[[mset]],[[maskZ]],[[maskZ]],[[maskZ]]]
  if (zmask) then diff=[[[diff]],[[datumZ]]]
  if (xyzmask) then diff=[[[diff]],[[datumZ]],[[datumZ]],[[datumZ]]]
  if (zmask) then ndiff=[[[ndiff]],[[sigmaZ]]]
  if (xyzmask) then ndiff=[[[ndiff]],[[sigmaZ]],[[sigmaZ]],[[sigmaZ]]]
  
;  tv_multi,diff[*,*,2],title='mask data'
;  tv_multi,mset[*,*,2],title='mask mask'
;  tv_multi,ndiff[*,*,2],title='mask noise'
  
  ; switch to MSB unless DN is given
  if (dn eq 0) then begin
      diff=diff*msbperdn
      ndiff=ndiff*msbperdn
  end

  if (n_elements(mxit) ne 0) then pxwrapper.sp.mxit=mxit

  save,pxwrapper,diff,ndiff,mset,rawdata,priordata,$
    file=cmedate+'_'+int2str(N)+'_inputdata.sav'
  
  ; ** CHECK S/N

  if (checksn) then begin
;      if (lasco) then $
;        sn=ratio(diff[*,*,0:2],ndiff[*,*,0:2])*mset[*,*,0:2] else $
;        sn=ratio(diff[*,*,0:1],ndiff[*,*,0:1])*mset[*,*,0:1]
;      tv_multi,fixed=256,sn > 1,/log,header='S/N ratio > 1',/showrange,$
;        title=pxwrapper.sc.gname,/raw

      sn=ratio(diff,ndiff)*mset > 1.0

      tsetd=make_array(/string,pxwrapper.sr.Nd,value='Data')
      tsetn=make_array(/string,pxwrapper.sr.Nd,value='Noise')
      tsetm=make_array(/string,pxwrapper.sr.Nd,value='Masks')
      tsets=make_array(/string,pxwrapper.sr.Nd,value='S/N')

      tset=[tsets,tsetm,tsetn,tsetd]

      tv_multi,title='data/noise/masks/sn',labels=tset,$
        [[[sn]],[[mset]],[[ndiff]],[[diff]]],$
        fixed=128,/raw,/log,/local,line=pxwrapper.sr.Nd,/showrange,oom=2

;      print,'models (and Istds) are: ',$
;        pxwrapper.sr.imodel,'(',pxwrapper.sr.istd,')'
  end

;  tv_multi,res=3,diff,title='one',/local,/log

  if (setuponly) then return

  p_generate_pixon,'classic',pxwrapper,oPixon
;  help,diff,ndiff,mset,pxwrapper.sr.Nd
  p_insert_data,pxwrapper,diff,oPixon,sigma=ndiff,masks=mset

  print,'Now starting Pixon reconstruction.'
  
  if (pseudo) then oPixon->pseudo,/restart else oPixon->full,/restart
 
  solution_image=oPixon->get(/im,/nop)
  solution_renders=stereo_project(pxwrapper,oPixon)
  input_data=oPixon->get(/data,/nop)

  savefile=cmedate+'_'+int2str(N)+'_solved.sav'
  
  p_update_runstats,pxwrapper,oPixon

  save,pxwrapper,oPixon,solution_image,solution_renders,input_data,$
    file=savefile
  
  print,savefile		 

  if (standardplots) then p_standardplots,/best,filein=savefile,/noscreen

end