;+
;
; Project   : s3drs and Pixon reconstruction software
;                   
; Name      : rtwl_wrapper
;               
; Purpose   : part of the s3drs GUI
;               
; Explanation: (part of the testing and legacy code, not documented herein)
;               
; Use       : 
;    
; Inputs    : 
;               
; Outputs   : 

; Keywords  :
;               
; Common    : 
;               
; Restrictions: none
;               
; Side effects: 
;               
; Category    : 
;               
; Prev. Hist. : (see CVS logs)
;
; Written     : Sandy Antunes, NRC, 2007-2009
;               
;-            
; simple little routine to render using RTWL and return a set
; of 'datum' images, mimicking Cor2
;
; for now, simple, later I will add other parameters
;

FUNCTION rtwl_wrapper,pxwrapper,timage,modparam=modparam,interp=interp,$
                      polonly=polonly,neonly=neonly,noisy=noise,$
                      rsunfactor=rsunfactor,gamma=gamma,theta=theta,$
                      distance=distance


;rtinitenv,forcelibfile='/home/antunes/dev/secchi/lib/linux/x86/libraytrace.so'
;;;rtinitenv,forcelibfile='/home/antunes/libraytrace.so'


  ; interp 25 = no interpolation, interp 26 = trilinear interpolation
  interp=KEYWORD_SET(interp)
  if (interp) then it=26 else it=25

  neonly=KEYWORD_SET(neonly)
  polonly=KEYWORD_SET(polonly)
  noisy=KEYWORD_SET(noisy)
  if (n_elements(rsunfactor) eq 0) then rsunfactor=1.0

  if (n_elements(gamma) eq 0) then gamma=pxwrapper.sr.gamma
  if (n_elements(theta) eq 0) then theta=pxwrapper.sr.theta

  if (n_elements(timage) ne 0) then begin
      ; convert timage to our modparam structure
      reorder=[3,2,1]; converts RWTL standard to Pixon standard
      N=n_elements(timage[0,0,*])
      rsun=pxwrapper.sr.l0tau[0]*rsunfactor
      ; Arnaud says use center N/2, not (N-1)/2
      modparam=[N,N,N,N/2,N/2,N/2,rsun,$
                reform(rearrange(timage,reorder),n_elements(timage))]

  end else if (n_elements(modparam) ne 0) then begin
      N=long(modparam[0])       ; get actual data size
  end else begin
      print,'No input image or modparam, exiting without rendering'
      return,-1
  end
;  help,timage,modparam[7]

  rsunrad=0.0046524             ; from Allen

;  pxwrapper.sr.Ncube=N; update structure just to be sure we are in sync

  datum=fltarr(N,N,n_elements(gamma))

  for i=0,n_elements(gamma)-1 do begin

      if (n_elements(distance) eq 0) then $
        obspos=[0.0,0.0,0.0-(pxwrapper.sr.l0*pxwrapper.sr.ll[i])] else $
        obspos=[0.0,0.0,0.0-distance]

      ; note that pxwrapper.sr.l0tau=pixel in Rsun
      ;   so pxwrapper.sr.l0tau*N = 30 for Cor2
      if (n_elements(pxwrapper) ne 0) then $
        fovcor2=rsunrad * pxwrapper.sr.l0tau[i] else $
        fovcor2=rsunrad * 30.0 / float(N)

      if (n_elements(pxwrapper) ne 0) then $
        limbdark=pxwrapper.sr.limbdark[i] else $
        limbdark=0.804

      ; using angles gamma,0,-90 puts RTWL into a LH, CCW frame
      ; where the rotation gamma is around the 'j' axis
      ; rotate by 0,0,-90 to put into LH, CCW frame, then add a
      ; -90 degree rotation so RTWL and Pixon match up

;      if (neonly) then print,'neonly *********' else print,'no ********'

      ; cludgy use of quiet because raytracewl does not accept quiet=0 as valid


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

      xang=gamma[i]*!dtor
      zang=0.0-theta[i]*!dtor

      if (noisy) then $
        raytracewl,sbt,spt,sne,$
        imsize=[N,N],losnbp=long(N),losrange=[-15,15],$
        modelid=it,modparam=modparam,neang=[xang,zang,0],$
        fovpix=fovcor2,obspos=obspos,limbdark=limbdark else $
        raytracewl,sbt,spt,sne,$
        imsize=[N,N],losnbp=long(N),losrange=[-15,15],$
        modelid=it,modparam=modparam,neang=[xang,zang,0],$
        fovpix=fovcor2,obspos=obspos,limbdark=limbdark,/quiet

;fovpix=0.00224954 ;fovpix=(1.0/N)*!dtor; 0.00816814; fovpix=7.02980e-5
; 0.00224954/!dtor=0.128889, *64=8.2489
; at fov=30 rsun and the sun 0.25 deg * !dtor= is 0.004363 so
;        30*0.25*!dtor = fov is 0.1309 radians across,
;        0.1309/64pixels = 0.00204531 
; Allen gives semi-diamater of 0.0046524 radians
;   so 30.0*0.0046524/64= 0.00218081

; using his example of (1.0/N)*!dtor yields 0.00816814,
; giving his main.cpp value is 7.02980e-5, yet /cor yields 0.00224954?


      ; "native" RTWL schema
;      raytracewl,sbt,imsize=[N,N],losnbp=long(N),losrange=[-20,20], $
;        modelid=25,modparam=modparam,neang=[gamma[i],0,0]*!dtor,/cor2

      
      if (polonly) then $
        datum(*,*,i)=spt.im else if (neonly) then $
        datum(*,*,i)=sne.im else datum(*,*,i)=sbt.im

  end

  return,datum

END