;+
;
; 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 flyby_reorder,dset,lset,nper,skip=skip

  if (n_elements(nper) eq 0) then nper=6; default is 6 per row
  if (n_elements(skip) eq 0) then skip=1; default is use each element

  ; subselect if requested
  if (skip ne 1) then begin
      Nd=n_elements(dset)
      dset2=ptrarr(1+Nd/skip)
      lset2=fltarr(3,1+Nd/skip)
;      print,Nd,Nd/skip
      j=0
      for i=0,n_elements(dset)-1,skip do begin
          dset2[j]=dset[i]
          lset2[*,j]=lset[*,i]
          ++j
      end
      dset=dset2
      lset=lset2
  end

  ; do a 'best top-down re-sort'
  dtemp=ptrarr(nper)
  ltemp=fltarr(3,nper)

  nrows=n_elements(dset)/nper
  dset=dset[0:nrows*nper-1]; truncate to even rows

  itop=nrows-1
;  print,'rows: ',nrows,', per row: ',nper,', total ele: ',nrows*nper,$
;    ' vs actual data set: ',n_elements(dset),' top row id is ',itop

  for ir=0,(nrows/2)-1 do begin
      
      i0=ir*nper
      i1=i0+(nper-1)
      i2=itop*nper
      i3=i2+(nper-1)
;      print,i0,i1,' to ',i2,i3

      dtemp=dset[i0:i1]
      dset[i0:i1]=dset[i2:i3]
      dset[i2:i3]=dtemp

      ltemp=lset[*,i0:i1]
      lset[*,i0:i1]=lset[*,i2:i3]
      lset[*,i2:i3]=ltemp

      --itop

  end

END

FUNCTION flyby_default_views

  views=[  [0.0,  0.0, 1.0],$
           [20.0, 0.0, 1.0],$
           [40.0, 0.0, 1.0],$
           [60.0, 0.0, 1.0],$
           [80.0, 0.0, 1.0],$
           [90.0, 0.0, 1.0],$
           [0.0, 0.0, 1.0],$
           [0.0, 20.0, 1.0],$
           [0.0, 40.0, 1.0],$
           [0.0, 60.0, 1.0],$
           [0.0, 80.0, 1.0],$
           [0.0, 90.0, 1.0],$
           [10.0, 90.0, 1.0],$
           [20.0, 90.0, 1.0],$
           [40.0, 90.0, 1.0],$
           [60.0, 90.0, 1.0],$
           [80.0, 90.0, 1.0],$
           [90.0, 90.0, 0.9],$
           [90.0, 90.0, 0.8],$
           [90.0, 90.0, 0.7],$
           [90.0, 90.0, 0.6],$
           [90.0, 90.0, 0.5]]

  return,views

END

FUNCTION flyby_default_zooms

  views=[  [0.0,  0.0, 1.0],$
           [0.0,  0.0, 0.9],$
           [0.0,  0.0, 0.8],$
           [0.0,  0.0, 0.7],$
           [0.0,  0.0, 0.6],$
           [0.0,  0.0, 0.5]]

  return,views

END

FUNCTION flyby_default_cube,N

  imgcube=fltarr(N,N,N)
  imgcube[N/4:N/2,N/4:N/2,N/4:N/2]=10.0^6

  return,imgcube

END


; skip means to only use each 'skip'th element
PRO print_flyby,filename,delay=delay,local=local,inverse=inverse,$
                noscreen=noscreen,skip=skip,nper=nper

  noscreen=setyesno(noscreen)

  if (n_elements(nper) eq 0) then nper=6
  if (n_elements(skip) eq 0) then skip=1

  if (file_exist(filename)) then begin
      restore,filename; contains ptr arry 'dset' and angles 'views'

      flyby_reorder,dset,views,nper,skip=skip; put into order

      nele=n_elements(views[0,*])
      lset=strarr(nele)
      for i=0,n_elements(views[0,*])-1 do $
        lset[i]=strjoin(int2str(floor(views[*,i])),',')

;      if (noscreen eq 0) then $
;        tv_multi,fixed=64,oom=2,/cap,dset[0:nele-2],/slide,winid=winid,$
;        /log,labels=lset,line=nper,local=local,inverse=inverse
      if (noscreen eq 0) then $
        tv_multi,fixed=64,oom=1,/cap,dset,/slide,winid=winid,$
        /log,labels=lset,line=nper,local=local,inverse=inverse
      tv_multi,fixed=96,oom=2,/cap,dset,/slide,winid=winid,$
        /log,labels=lset,line=nper,/ps,file=filename+'.ps',local=local,$
        inverse=inverse

      print,'wrote ',filename+'.ps'
  end

END

PRO show_flyby,fname1,fname2,fname3,fname4,fname5,fname6,fname7,fname8,$
               delay=delay,log=log,oom=oom,cap=cap,fixed=fixed,$
               movie=movie,flist=flist,mask=mask,maxval=maxval

  movie=setyesno(movie)

  if (n_elements(mask) eq 0) then mask=1.0

  fset=[fname1]
  if (n_elements(fname2) ne 0) then fset=[fset,fname2]
  if (n_elements(fname3) ne 0) then fset=[fset,fname3]
  if (n_elements(fname4) ne 0) then fset=[fset,fname4]
  if (n_elements(fname5) ne 0) then fset=[fset,fname5]
  if (n_elements(fname6) ne 0) then fset=[fset,fname6]
  if (n_elements(fname7) ne 0) then fset=[fset,fname7]
  if (n_elements(fname8) ne 0) then fset=[fset,fname8]

  ndatasets=n_elements(fset)
  dptrset=ptrarr(ndatasets)
  dptr=ptrarr(ndatasets)

  dmax=0.0

  if (file_exist(fname1)) then begin

      for iset=0,ndatasets-1 do begin
          restore,fset[iset] ; contains ptr arry 'dset' and angles 'views'
          if (n_elements(nframes) eq 0) then nframes=n_elements(dset)

          for i=0,nframes-1 do *(dset[i])=*(dset[i])*mask
          for i=0,nframes-1 do dmax=max([dmax,max(*(dset[i]))])
          dptrset[iset]=ptr_new(dset)
      end
      print,'dmax is ',dmax

      if (n_elements(maxval) ne 0) then dmax=maxval

      for i=0,nframes-1 do begin
          head=int2str(floor(views[0,i]))+','+int2str(floor(views[1,i]))

          for iset=0,ndatasets-1 do begin
              dtemp=*(dptrset[iset])
              dptr[iset]=dtemp[i]
          end

          tv_multi,fixed=fixed,oom=oom,cap=cap,dptr,$
            /slide,winid=winid,/local,log=log,head=head,maxval=dmax
          if (movie) then begin
              ftemp=fname1 + '.' + int2str(i,3) + '.gif'
              write_gif,ftemp,tvrd()
              if (i eq 0) then flist=[ftemp] else $
                flist=[flist,ftemp]
          end
          if (n_elements(delay) ne 0) then wait,delay
      end

  end

END

; best default is /grandtour, which does a 180 ecliptic then a 180 polar

PRO flyby,imgcube,rtwl=rtwl,both=both,views=views,N=N,zoom=zoom,mask=mask,$
          timing=timing,mass=mass,delay=delay,savefile=savefile,$
          eclipticloop=eclipticloop,zloop=zloop,addaxis=addaxis,cap=cap,$
          grandtour=grandtour,test=test,angleper=angleper,noisy=noisy

  start=systime(/seconds)

  rtwl=setyesno(rtwl)
  both=setyesno(both)
  zoom=setyesno(zoom)
  timing=setyesno(timing)
  mass=setyesno(mass)
  eclipticloop=setyesno(eclipticloop)
  zloop=setyesno(zloop)
  grandtour=setyesno(grandtour)
  mask=setyesno(mask)
  addaxis=setyesno(addaxis)
  cap=setyesno(cap)
  test=setyesno(test)
  noisy=setyesno(noisy)
  if (n_elements(angleper) eq 0) then angleper=10; spacing for autoloops

  if (n_elements(savefile) eq 0) then savefile='flyby.out'

  nele=floor(360/angleper) + 1 ; default we may need for autoloops
  if (test) then views=make_array(nele,/float,/index)*angleper

  if (grandtour) then begin
      ; grand tour is ecl 0-180, then polar 0-180, then both to return
      nele=3*floor(180/angleper) + 4
      views=fltarr(3,nele)
      i=0
      for angle = 0,180,angleper do begin
          views[0,i]=angle
          views[1,i]=0.0
          ++i
      end
      for angle=0,180,angleper do begin
          views[0,i]=180
          views[1,i]=angle
          ++i
      end
      for angle = 180,360,angleper do begin
          views[0,i]=angle
          views[1,i]=angle
          ++i
      end
      views[0,i]=0.0
      views[1,i]=0.0
      ++i

;      print,i,nele
      views[2,*]=1.0; all at d=1

  end

  if (n_elements(N) eq 0) then N=64
  if (n_elements(views) eq 0) then $
    if (eclipticloop or zloop or grandtour) then $
    views=make_array(nele,/float,/index)*angleper else $
    if (zoom) then views=flyby_default_zooms() else views=flyby_default_views()
  if (n_elements(imgcube) eq 0) then imgcube=flyby_default_cube(N)
  if (n_elements(delay) eq 0) then delay=0

  N=n_elements(imgcube[0,0,*])

  ; if given just a 1-D array, it's an spin and needs rho/z/d
  if (size(views,/n_dim) eq 1) then begin
      tviews=fltarr(3,n_elements(views))
      rho=views
      z=0.0
      if (zloop) then z=views
      if (eclipticloop) then rho=views
      tviews[2,*]=1.0
      tviews[1,*]=z
      tviews[0,*]=rho
      views=tviews
  end

  iviews=n_elements(views[0,*])

  ; if using rtwl, made our modparam reorder once now for speed
  if (rtwl or both) then begin
      reorder=[3,2,1]; converts RWTL standard to Pixon standard
      modparam=[N,N,N,N/2,N/2,N/2,0.0,$
                reform(rearrange(imgcube,reorder),n_elements(imgcube))]
  end

  if (rtwl) then lset=['RTWL'] else if both then lset=['RTWL','PR'] else $
    lset=['PR']

  if (mask) then maskd=p_getmask('COR2',N) else maskd=fltarr(N,N)+1

  dset=ptrarr(iviews)

  if (addaxis) then begin
      ; use a 64x64x64 axis cube
      Na=64
      axisimage=fltarr(64,64,64)
      axisimage[32,32:40,32]=1
      axisimage[32:40,32,32]=1
      axisimage[32,32,32:40]=1
  end

  for ifr=0,iviews-1 do begin
      prep_3drv,pxfly,N=N,/cor2,$
        gamma=views[0,ifr],theta=views[1,ifr],LL=views[2,ifr]

      if (rtwl or both) then begin
          modparam[6]=pxfly.sr.l0tau[0]
          datr=maskd * rtwl_wrapper(pxfly,modparam=modparam,neonly=mass)
      end
      if (rtwl eq 0 or both) then begin
          pxfly.sr.tau[0]=0.0; switch to infinite geometry for zooming
          if (mass) then pxfly.sr.polar[0]='n'
          datp=maskd * pr_render(pxfly,imgcube,0)
      end

      if (n_elements(dmax) eq 0) then dmax=max(datp); initialize maximum
      if (cap) then datp=datp < dmax
      if (cap and (rtwl+both)) then datr=datr < dmax

      ; horribly slow way to add a ring
      if (addaxis) then begin
          gamma=pxfly.sr.rho[0]
          theta=pxfly.sr.theta[0]
          timage=axisimage
          cuberot3d,timage,[0.0,theta,0.0-gamma],/reverse
          d=ijkproj(timage,'jk'); jk best, ij,ji,ki,kj wrong, ik backwards
;          d[where (d gt 0)] = 1.0; make it uniform, with no depth
          d = d * max(datp)
          datp += rebin(d,N,N)
      end

      if (rtwl) then dat=datr else $
        if (both) then dat=[[[datr]],[[datp]]] else $
        dat=datp

      dset[ifr]=ptr_new(datp)

      head=int2str(floor(views[0,ifr]))+','+int2str(floor(views[1,ifr]))

      if (noisy) then print,'generated ',ifr+1,' of ',iviews

;      if (noscreen eq 0) then $
;        tv_multi,fixed=256,oom=2,/cap,dat,/slide,$
;        winid=winid,/local,labels=lset,/log,head=head
      if (delay ne 0) then wait,delay
  end

  save,file=savefile,views,dset
  print,"show_flyby,'"+savefile+"' or print_flyby,'"+savefile+"'"

  if (timing) then print,(systime(/seconds)-start)/60.0,' minutes'

END