;+
;
; 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
;               
;-            
; besides different N and choices, you can set /blowup to preserve
; the same pixel scale across different N, rather than keeping the
; same FOV (and thus reducing pixel scale for larger N)

PRO ctrtests,N=N,choice=choice,blowup=blowup,gamma=gamma

  blowup=setyesno(blowup)
  if (n_elements(N) eq 0) then N=4; or N=16
  if (n_elements(choice) eq 0) then choice=1; 1=lit, 2=asym, 3=saito
  if (choice le 0 or choice gt 3) then stop,'choice must be 1, 2 or 3'
  if (n_elements(gamma) eq 0) then gamma=0.0

  if (choice eq 1) then print,'lit cube' else $
    if (choice eq 2) then print,'asym cube' else print,'saito cube'

  n1=(N-1)/2
  n2=N/2
  f1=(N-1)/2.0
  f2=N/2.0
  print,'pixel (N-1)/2:N/2 is ',n1,n2
  print,'coord (N-1)/2:N/2 is ',f1,f2

  if (choice eq 1) then lit=fltarr(N,N,N)
  if (choice eq 1) then lit[n1:n2,n1:n2,n1:n2]=10^6.5
  if (choice eq 2) then asym=fltarr(N,N,N)
  if (choice eq 2) then asym[n1,n1,n1]=10^6.5
  if (choice eq 3) then saito=pixonmodels(N,'saito')

  if (choice eq 1) then timage=lit else $
    if (choice eq 2) then timage=asym else timage=saito

  print,'cube ctr: ',timage[n1:n2,n1:n2,n1:n2]

  rsun=30.0/float(N)            ; e.g. for N=64, 0.46875
  if (blowup) then rsun=3.75; fixed ratio, based on N=8
  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=[N,N],losnbp=long(N),losrange=[-15,15],$
    modelid=25,modparam=modparam,neang=[gamma,0,0]*!dtor,$
    fovpix=fov,obspos=[0.0,0.0,-214.94299],/quiet

  rtwl=sbt.im

  prep_3drv,pxwrapper,N=N,gamma=[gamma]
  pxwrapper.sr.tau=0.0; set infinite geometry
  if (blowup) then pxwrapper.sr.l0tau=rsun
  if (blowup) then pxwrapper.sr.tau=pxwrapper.sr.l0tau/(2.0*pxwrapper.sr.l0)
  pixon=pr_render(pxwrapper,timage)

  ;'r' is centered on (N-1)/2.0, so pixel(n1)=pixel(n2) and 0 is between them
  r=(make_array(N,/float,/index)-((N-1)/2.0))*rsun

  ir1=max(where(r lt 2.0))+1
  ir2=min(where(r gt -2.0))-1

; now print results

  print,'RTWL center coords: ',N/2,N/2
  print,'Pixon render cubectr coords: ',pxwrapper.sr.lla0
  print,'Pixon render optctr coords: ',p_getk0(pxwrapper)
  print,'Pixon render sunctr coords: ',p_getk0(pxwrapper,/sun)

  tv_multi,/raw,fixed=256,[[[sbt.im]],[[pixon]]],title='Renders',$
    labels=['RTWL','Pixon']

  print,'RTWL render ctr : ',$
    rtwl[n1,n1],rtwl[n1,n2],rtwl[n2,n1],rtwl[n2,n2]
  print,'Pixon render ctr: ',$
    pixon[n1,n1],pixon[n1,n2],pixon[n2,n1],pixon[n2,n2]

  if (choice ne 3) then return

  ; extra items for Saito check
  print,' '
  print,'RTWL far left edge : ',rtwl[0,n1],rtwl[0,n2]
  print,'RTWL far right edge: ',rtwl[N-1,n1],rtwl[N-1,n2]
  print,'Pixon far left edge : ',pixon[0,n1],pixon[0,n2]
  print,'Pixon far right edge: ',pixon[N-1,n1],pixon[N-1,n2]
  print,' '
  print,'RTWL r=',r[ir1],':',rtwl[ir1,n1],rtwl[ir1,n2]
  print,'RTWL r=',r[ir2],':',rtwl[ir2,n1],rtwl[ir2,n2]
  print,'Pixon r=',r[ir1],':',pixon[ir1,n1],pixon[ir1,n2]
  print,'Pixon r=',r[ir2],':',pixon[ir2,n1],pixon[ir2,n2]
  window,/free,retain=2,xsize=400,ysize=400

  !P.MULTI=[0,0,2]
  plot,r,rtwl[*,n1],title='RTWL slice'
  oplot,r,rtwl[*,n2]
  oplot,[-1,-1,1,1],[0,1,1,0]
  plot,r,pixon[*,n1],title='Pixon slice'
  oplot,r,pixon[*,n2]
  oplot,[-1,-1,1,1],[0,1,1,0]
  !P.MULTI=0

END