;+
;
; 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
;               
;-            
pro break_rtwl

  mycase=1

  if (mycase eq 0) then begin
     ; fastest seg fault yet!
      N=64
      timage=pixonmodels(N,'saito')
      rsun=30.0/float(N)        ; e.g. for N=64, 0.46875
      rsunrad=0.0046524         ; from Allen
      fov=rsunrad * rsun        ; e.g. for N=64, 0.0021808
      modparam=[N,N,N,N/2,N/2,N/2,rsun,$
                reform(timage,n_elements(timage))]
      raytracewl,sbt,spt,sne,$
        imsize=[64,64],losnbp=long(64),losrange=[-15,15],$
        modelid=25,modparam=modparam,neang=[0.0,0,0]*!dtor,$
        fovpix=fov,obspos=[0.0,0.0,-214.94299]

;,limbdark=0.805

  end else if (mycase eq 1) then begin
     ; fails

      restore,/v,'debug1a.sav'

      rtwl=rtwl_wrapper(pxwrapper,imgK,neonly=eonly)
      rtwlS=rtwl_wrapper(pxwrapper,timage,neonly=eonly); ,/interp)
     
      tv_multi,fixed=128,/log,rtwlS

  end else if (mycase eq 2) then begin
     ; same as 1 but not using rtwl wrapper, works
      restore,/v,'debug1a.sav'

      N=64
      reorder=[3,2,1]       ; converts RWTL standard to Pixon standard
      rsun=0.4687500
      modparam=[N,N,N,N/2,N/2,N/2,rsun,$
                reform(rearrange(imgK,reorder),n_elements(timage))]
      raytracewl,sbt,spt,sne,$
        imsize=[64,64],losnbp=long(64),losrange=[-15,15],$
        modelid=25,modparam=modparam,neang=[0.0,0,0]*!dtor,$
        fovpix=0.0021808,obspos=[0.0,0.0,-214.94299],limbdark=0.805
      
      modparam=[N,N,N,N/2,N/2,N/2,rsun,$
                reform(rearrange(timage,reorder),n_elements(timage))]
      raytracewl,sbt,spt,sne,$
        imsize=[64,64],losnbp=long(64),losrange=[-15,15],$
        modelid=25,modparam=modparam,neang=[0.0,0,0]*!dtor,$
        fovpix=0.0021808,obspos=[0.0,0.0,-214.94299],limbdark=0.805
      
      rtwlS=sbt.im
      
      tv_multi,fixed=128,/log,rtwlS
      
  end else begin
                                ;works
      restore,/v,'debug1b.sav'
      rtwlS=rtwl_wrapper(pxwrapper,timage,neonly=eonly) ;,/interp)
      
      tv_multi,fixed=128,/log,rtwlS
  end
  
end

pro rtwltest

  N=64

  prep_3drv,N=N,pxwrapper,/cor2

  dialme,'saito'
  restore,/v,'saito+saito_noise0saitoscl=0.00000.img'
  timage=congrid(timage,N,N,N)

  reorder=[3,2,1]; converts RWTL standard to Pixon standard
  print,pxwrapper.sr.l0tau
  rsun=0.4687500
  modparam=[N,N,N,N/2,N/2,N/2,rsun,$
            reform(rearrange(timage,reorder),n_elements(timage))]

  print,pxwrapper.sr.l0
  distance=214.94299
  obspos=[0.0,0.0,0.0-distance]

  print,pxwrapper.sr.gamma
  gamma=[0.0]

  print,pxwrapper.sr.limbdark
  limbdark=0.805

  it=25; interp

  fovcor2=rsun
  rsunrad=0.0046524             ; from Allen

  raytracewl,sbt,spt,sne,$
    imsize=[64,64],losnbp=long(64),losrange=[-15,15],$
    modelid=25,modparam=modparam,neang=[0.0,0,0]*!dtor,$
    fovpix=0.0021808,obspos=[0.0,0.0,-214.94299],limbdark=0.805

  im1=sbt.im

  print,N,it,modparam[0:6],gamma[0],fovcor2,obspos,limbdark

  raytracewl,sbt,spt,sne,$
    imsize=[N,N],losnbp=long(N),losrange=[-15,15],$
    modelid=it,modparam=modparam,neang=[gamma[0],0,0]*!dtor,$
    fovpix=fovcor2*rsunrad,obspos=obspos,limbdark=limbdark

  im2=sbt.im

  im3=rtwl_wrapper(ignore,modparam=modparam,gamma=gamma,distance=distance)

  im4=rtwl_wrapper(pxwrapper,timage)

  tv_multi,fixed=128,[[[im1]],[[im2]],[[im3]],[[im4]]],/local,/log,$
    labels=['rtwl hard','rtwl var','wrap modparam','wrap timage']

end


pro renderoffcomps

  rsun=[2.0,4.0,6.0,7.9]
  pixonA=[1.342,0.353,0.161,0.092]
  pixonB=[4244.37,1116.96,509.69,291.1]
  rtwlA=[0.0073,0.00193,0.000884,0.000505]
  rtwlB=[23.187,6.1226,2.796,1.598]
  analytic=[0.000771,0.000625,0.00051665,0.000434]


end


pro renderoff,N=N,rsun=rsun,plot=plot,saito=saito,factor=factor,$
              choice=choice,eonly=eonly,disablelimb=disablelimb,shell=shell,$
              ratioreturn=ratioreturn,mask=mask,skippixeltest=skippixeltest,$
              debug=debug,ps=ps

  on_error,2
  
  if (n_elements(N) eq 0) then N=64
  if (n_elements(factor) eq 0) then factor=6.5 ; or 6.31
  if (n_elements(rsun) eq 0) then rsun=2.0
  plot=setyesno(plot)
  eonly=setyesno(eonly)
  disablelimb=setyesno(disablelimb)
  debug=setyesno(debug)
  skippixeltest=setyesno(skippixeltest)
  shell=setyesno(shell)
  mask=setyesno(mask)
  saito=setyesno(saito)
  if (saito or shell) then plot=1
  ps=setyesno(ps)
  
  if (shell) then lab='Halfshell' else lab='Saito'
  if (mask) then m='masked_' else m=''

; choose a comp
; choice=1 is Pixon with and without thomson
; choice=2 is Pixon no-thom vs RTWL
; choice=3 is Pixon vs RTWL (default)
  if (n_elements(choice) eq 0) then choice=3
  
  angle=0.0
  sangle='0'
  
  ergsper=2.457e+06
  
  prep_3drv,pxwrapper,N=N,gamma=[angle],/cor2
  
  pxwrapper.sr.imodel=10
;print,pxwrapper.sr.imodel,pxwrapper.sr.polar
;pxwrapper.sr.lla0=long(N/2)
  
  if (disablelimb) then pxwrapper.sr.chl=0 ; turns off limb darkening
  
  pxwrapper2=pxwrapper
  pxwrapper2.sr.polar='n'
;pxwrapper.sr.tau=0
  
;pxwrapper.sr.imodel=3
  
  PixelsPerRsun=1.0/pxwrapper.sr.l0tau[0]
  
  center=pxwrapper.sr.lla0[0]
  
  r2A=floor(center + PixelsPerRsun*rsun)
  r2B=floor(center + PixelsPerRsun*rsun)+1
  ractualA=(r2A-center)/PixelsPerRsun
  ractualB=(r2B-center)/PixelsPerRsun
  if ( abs(ractualA-rsun) le abs(ractualB-rsun))then r2=r2a else r2=r2b
  
  ractual=(r2-center)/PixelsPerRsun
  
  if (skippixeltest eq 0) then begin
      
      imgK=fltarr(N,N,N)
      
;;if (skippixeltest eq 0) then begin
;factor=6.31
;factor=6.5
      
      imgK[N/2,r2,N/2]=10^factor
;threeview,fixed=128,imgK,/surp,imgset
      
      pr_thom=pr_render(pxwrapper,imgK)
      pr_nothom=pr_render(pxwrapper2,imgK)
      rtwl=rtwl_wrapper(pxwrapper,imgK,neonly=eonly)
      pr_thom=rtwl
      pr_nothom=rtwl
      
; pick which pair to compare here
      ren1=pr_thom
      ren2=pr_nothom
      ren3=rtwl
      
      ir1=where(ren1[*,*,0] ne 0)
      ir2=where(ren2[*,*,0] ne 0)
      ir3=where(ren3[*,*,0] ne 0)
      
;print,'lit pixels for Pixon: ',ir1
;print,'lit pixels for Pixon: ',ir2
;print,'lit pixels for RTWL: ',ir3
      
      rows1=ir1 MOD N
      cols1=ir1/N
      rows2=ir2 MOD N
      cols2=ir2/N
      rows3=ir3 MOD N
      cols3=ir3/N
      
      ne1=n_elements(ir1)
      ne2=n_elements(ir2)
      ne3=n_elements(ir3)
      
; find overlapping ROI
      rows=[min([rows1[0],rows2[0],rows3[0]]),$
            max([rows1[ne1-1],rows2[ne2-1],rows3[ne3-1]])]
      cols=[min([cols1[0],cols2[0],cols3[0]]),$
            max([cols1[ne1-1],cols2[ne2-1],cols3[ne3-1]])]
      
      cosy = cos(asin(1.0/r2))
      analytic = (4.0/3.0) - ( cosy + (1.0/3.0)*(cosy^3.0) )
      
; PRINT RESULTS
      sfactor=strtrim(string(factor),2)
      
      
      if (eonly) then print,'Ne only'
      if (disablelimb) then print,'(no limb darkening for Pixon)'
      
      tp=total(ren1[rows1[0]:rows1[ne1-1],cols1[0]:cols1[ne1-1]])*ergsper
      tr=total(ren3[rows3[0]:rows3[ne3-1],cols3[0]:cols3[ne3-1]])*ergsper
      
      print,'At equatorial R='+strtrim(string(rsun),2)+$
        ' for log('+sfactor+')e-, B/B0 @ 7000A = 1.7E-08',$
        'at an actual image cube distance of ',ractual,'                  '
      print,$
        'Pixon non-zero values: ',ren1[rows1[0]:rows1[ne1-1],cols1[0]:cols1[ne1-1]],$
        'totalling',tp
;  ' vs                                                            ',$
;'Pixon non-thompson non-zero values: ',$
;  ren2[rows2[0]:rows2[ne2-1],cols2[0]:cols2[ne2-1]],$
; 'totalling',total(ren2[rows2[0]:rows2[ne2-1],cols2[0]:cols2[ne2-1]])*ergsper,
      print,' vs                                                            ',$
        'RTWL non-zero values: ',ren3[rows3[0]:rows3[ne3-1],cols3[0]:cols3[ne3-1]],$
        'totalling',tr
      
      print,'vs analytic factor ',analytic
      
      print,'*** ratio Pixon/RTWL: ',tp/tr,' ,Pixon/analytic: ',tp/analytic,$
        ' RTWL/analytic: ',tr/analytic
      
      if (plot eq 0) then return
;pause
      
                                ; now PLOT
      
;lset=['rtwl '+sangle,'pixon no-thom','pixon '+sangle]
      lset=['rtwl '+sangle,'pixon '+sangle,'diff']
;dset=[[[ren3]],[[ren2]],[[ren1]]]
;tv_multi,fixed=256,line=1,dset,labels=lset,/unique,/show,/raw
      
      dsub1=ren1[rows[0]-1:rows[1]+1,cols[0]-1:cols[1]+1]
      dsub2=ren2[rows[0]-1:rows[1]+1,cols[0]-1:cols[1]+1]
      dsub3=ren3[rows[0]-1:rows[1]+1,cols[0]-1:cols[1]+1]
      
      dsubd=dsub1-dsub3
      
;dsub=[[[dsub3]],[[dsub2]],[[dsub1]]]
      dsub=[[[dsub3]],[[dsub1]],[[dsubd]]]
      
      tv_multi,fixed=128,dsub,labels=lset,/show,/border,/raw,$
        header='Lit Area',ps=ps,/land,file=m+'renderoff_lit.ps'
      
      
;window,/free,retain=2
;!P.MULTI=[0,0,3]
;plot,ren1[rows[0]:rows[1],cols[0]:cols[1]],title='Pixon',font=7
;plot,ren2[rows[0]:rows[1],cols[0]:cols[1]],title='Pixon non-thom',font=7
;plot,ren3[rows[0]:rows[1],cols[0]:cols[1]],title='RTWL',font=7
;!P.MULTI=0
      
  end
  
  if (saito eq 0 and shell eq 0) then return
  
; now compare same regions using a saito or an offset half shell
  if (shell) then timage=pixonmodels(N,'hshell2') else $
    timage=pixonmodels(N,'saito')
  
  print,'equatorial log(Ne) at r=2 is 6.31'
  print,'log(Ne) density at rx=',ractual,' is ',alog10(timage[r2,N/2,N/2])
  print,'log(Ne) density at ry=',ractual,' is ',alog10(timage[N/2,r2,N/2])
  
; shell is offset so choose a new line
  if (shell) then r2s=(N/4+N/8) else r2s=r2
  if (shell) then slab='slice across N/4' else slab='slice across r=2Rsun'
  
  pr_thomS=pr_render(pxwrapper,timage)
  pr_nothomS=pr_render(pxwrapper2,timage)
  if (debug) then save,pxwrapper,timage,imgK,eonly,file='debug1a.sav'
  if (debug) then save,pxwrapper,timage,eonly,file='debug1b.sav'
  rtwlS=rtwl_wrapper(pxwrapper,timage,neonly=eonly) ;,/interp)
  if (debug) then save,pxwrapper,timage,eonly,rtwlS,file='debug2.sav'
  
  if (debug) then tv_multi,rtwlS,fixed=128,/log,title='debug rtwl'
  if (debug) then tv_multi,pr_thomS,fixed=128,/log,title='debug pr'
  
; choose a comp
; choice=1 is Pixon with and without thomson
; choice=2 is Pixon no-thom vs RTWL
; choice=3 is Pixon vs RTWL
  if (n_elements(choice) eq 0) then choice=3
  if (choice eq 2) then data1=pr_nothomS else data1=pr_thomS
  if (choice eq 1) then data2=pr_nothomS else data2=rtwlS
  if (choice eq 1) then ll=['Pixon','Pxn no-thom'] else $
    if (choice eq 2) then ll=['Pxn no-thom','RTWL'] else $
    ll=['Pixon','RTWL']
  
  if (debug) then tv_multi,fixed=128,[[[data1]],[[data2]]],title='debug',/local
  
  ratiop=ratio(data1,data2)
  
  ratioreturn=ratiop            ; optional return value
  
  lrset=['log(ratio)',ll[1]+' '+sangle,ll[0]+' '+sangle]
  if (mask) then maskdat=p_getmask(sr=pxwrapper.sr) else maskdat=fltarr(N,N)+1.0
  if (mask) then for i=0,n_elements(lrset)-1 do lrset[i]=lrset[i]+' masked'
  
  dsetobj=[[[alog10(ratiop)>0]],[[data2]],[[data1]]]
  tv_multi,fix=128,dsetobj*[[[maskdat]],[[maskdat]],[[maskdat]]],$
    title=lab,labels=lrset,/log,line=1,/show,/raw,header=lab,$
    ps=ps,/land,file=m+'renderoff_full.ps' ;,/reverse
  
  if (skippixeltest eq 0) then begin 
      
      dsub1obj=data1[rows[0]-1:rows[1]+1,cols[0]-1:cols[1]+1]
      dsub2obj=data2[rows[0]-1:rows[1]+1,cols[0]-1:cols[1]+1]
      dsubratio=ratio(dsub1obj,dsub2obj)
      dsubobj=[[[dsubratio]],[[dsub1obj]],[[dsub2obj]]]
      dlset=['log(ratio)',ll[1],ll[0]]
      tv_multi,fixed=128,dsubobj,$
        title=lab,label=dlset,/local,/show,/border,/raw,/log,$
        header=lab+' Subregion '+ll[0]+'/'+ll[1],ps=ps,/land,$
        file=m+'renderoff_pixeltest.ps'
      
;window,/free,retain=2
;!P.MULTI=[0,0,2]
;plot,data1[rows[0]:rows[1],cols[0]:cols[1]],title=ll[0],font=7
;plot,data2[rows[0]:rows[1],cols[0]:cols[1]],title=ll[1],font=7
;!P.MULTI=0
      
      
  end
  
  current_device=!d.name
  
  if (ps) then begin
      set_plot,'PS'
      device,filename=m+'renderoff_slices_graph.ps',/color,bits_per_pixel=8,$
        xsize=10.5,ysize=8.0,/inches,xoffset=0.5,yoffset=0.5,/landscape,$
        font_size=9
  end else begin
      window,/free,retain=2,xsize=600,ysize=400
      flip_colors
  end
  !P.MULTI=[0,0,3]
  plot,maskdat[*,r2s]*data2[*,r2s],title=slab,font=7,yrange=[1.5e-10,2.5e-10]
  oplot,maskdat[*,r2s]*data1[*,r2s]
  plot,alog(maskdat[*,r2s]*data2[*,r2s]),title='log '+slab,font=7
  oplot,alog(maskdat[*,r2s]*data1[*,r2s])
  plot,maskdat[*,r2s]*ratiop[*,r2s],title='Ratio '+ll[0]+'/'+ll[1],font=7
  !P.MULTI=0
  
  if (ps) then begin
      device,/CLOSE_FILE
      set_plot,current_device
  end else begin
      flip_colors
  end
  
; no values in middle for shell so return early, for saito add one more plot
  if (shell) then return
  
  c2=N/2
  
  if (ps) then begin
      set_plot,'PS'
      device,filename=m+'renderoff_midslice_graph.ps',/color,bits_per_pixel=8,$
        xsize=10.5,ysize=8.0,/inches,xoffset=0.5,yoffset=0.5,/landscape,$
        font_size=9
  end else begin
      window,/free,retain=2,xsize=600,ysize=400
      flip_colors
  end
  !P.MULTI=[0,0,3]
  plot,maskdat[*,c2]*data1[*,c2],title='slice across middle',font=7
  oplot,maskdat[*,c2]*data2[*,c2]
  plot,alog(maskdat[*,c2]*data1[*,c2]),title='log slice across middle',font=7
  oplot,alog(maskdat[*,c2]*data2[*,c2])
  plot,maskdat[*,c2]*ratiop[*,c2],title='Ratio '+ll[0]+'/'+ll[1],font=7
  !P.MULTI=0
  
  if (ps) then begin
      device,/CLOSE_FILE
      set_plot,current_device
  end else begin
      flip_colors
  end
  
  
end