;+
;
; 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
;               
;-            
; test to check our Saito model versus our hand-calculated values at r=2Rsun

PRO saitoplot_plot,r,saito,pr,r2,N,Nstr

      window,/free,retain=2
      !P.MULTI=[0,0,2]
      plot,r,alog10(saito[r2,N/2,*]),title='log(Ne), N='+Nstr,font=7,$
        yrange=[1,8]
      plot,r,pr[r2,*],title='Render, N='+Nstr,font=7,yrange=[0,2e-8]
      !P.MULTI=0

END

PRO saitoplot,N=N,plot=plot

  on_error,2

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

  fov=30.0/float(N)
  r=(make_array(N,/index)-N/2.0)*fov
  r2=floor((N/2.0)-2.0/fov)
  r2actual=((N/2.0)-r2)*fov
  Nstr=strtrim(string(N),2)+' r='+strtrim(string(r2actual),2)

  saito=pixonmodels(N,'saito')
  flat=fltarr(N,N,N)+10^6.5

  prep_3drv,pxwrapper,N=N,/cor2

  prs=pr_render(pxwrapper,saito)
  prf=pr_render(pxwrapper,flat)

  if (plot) then saitoplot_plot,r,saito,prs,r2,N,Nstr

  print,'N=',N,', want r=2, R(actual) is ',r2actual
  print,'  For Saito:'
  print,'    log(Ne) at r=2: ',alog10(saito[r2,N/2,N/2]),' vs theory index 6.5'
  print,'    observed B/B0 at (r=2Rsun,N/2): ',prs[r2,N/2],' vs theory 1.7e-8'
  print,'  For flat 10^6.5:'
  print,'    log(Ne) at r=2: ',alog10(flat[r2,N/2,N/2]),' vs theory index 6.5'
  print,'    observed B/B0 at (r=2Rsun,N/2): ',prf[r2,N/2],' vs theory 1.7e-8'

end