; docformat = 'rst'
;
; NAME: 
;    shapePhaseFunction
; PURPOSE:
;    Compute the phase function for lightcurve/spin in KOALA format
;
;+
; :Description:
;    Compute the phase function for lightcurve/spin in KOALA format
;
; :Categories:
;    Shape, Disk I/O, Lightcurve
;
; :Returns: An array of structure containing each lightcurve and its
;    phase function. Fields are::
;      .JD: The original JD of observation
;      .flux: The original flux 
;      .phase: The phase angle for each observation
;      .value: The associated value of the phase function
;
; :Params:
;    lc: in, required, type=structure/string
;      A LC structure, or the path to a LC file (DAMIT format).
;    spin: in, required, type=structure/string
;      A spin structure, or the path to a spin file (DAMIT format).
;
; :Examples:
;    Write a the first lightcurve of a previously larger data set::
;      IDL> lc = koalaReadLC( 'file.lc' )
;      IDL> koalaWriteLC, 'newfile.lc', lc[0]
;
; Uses:
;    readSpin, koalaReadLC
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in February 2017, B. Carry (OCA)
;-
function shapePhaseFunction, lc, spin
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--I.1-- Check Input Syntax  ------------------------------------------------------------------
  if n_params() lt 2 then begin
    message, /ioError, 'Syntax: phaseStructure = shapePhaseFunction( lc, spin )'
    return, -1
  endif

 ;--I.2-- Check Input Lightcurves  -------------------------------------------------------------
  if size(lc,/type) ne 8 then begin
    if test_file(lc,/read) then begin
      lc = koalaReadLC(lc)
    endif else begin
      message, /ioError, 'Cannot read input lightcurve file: '+strTrim(lc,2)
      return, -2
    endelse
  endif
  
 ;--I.3-- Check Input Spin  --------------------------------------------------------------------
  if size(spin,/type) ne 8 then begin
    if test_file(spin,/read) then begin
      spin = readSpin(spin)
    endif else begin
      message, /ioError, 'Cannot read input lightcurve file: '+strTrim(spin,2)
      return, -3
    endelse
  endif

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Build the Phase Function                           -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--I.1-- Number of Lightcurves & Output Variables  --------------------------------------------
  nbLC = n_elements( lc ) 
  empty = { jd:ptr_new(), flux:ptr_new(), phase:ptr_new(), value:ptr_new() }
  res = replicate( empty, nbLC )

 ;--I.2-- For Each Epoch  ----------------------------------------------------------------------
  for kLC=0L, nbLC-1 do begin


   ;--I.3-- Phase and Distances  ---------------------------------------------------------------
    phase = fltarr( lc[kLC].N )
    dObs  = sqrt( total( (*lc[kLC].obs)^2., 2 ) )
    dSun  = sqrt( total( (*lc[kLC].sun)^2., 2 ) )

    for kEp=0, lc[kLC].N-1 do begin
       vSun = reform((*lc[kLC].sun)[kEp,*] / dSun[kEp])
       vObs = reform((*lc[kLC].obs)[kEp,*] / dObs[kEp])
       phase[kEp] = acos(  transpose(vSun) # vObs ) / !DTOR
    endfor

   ;--I.4-- Phase Function  --------------------------------------------------------------------
    f = spin.phase[0] * exp( -(phase*!DTOR)/(spin.phase[1]) ) + spin.phase[2]*phase*!DTOR + 1.
    f /= mean(F)

   ;--I.5-- Fill Output Structure  -------------------------------------------------------------
    res[kLC].phase = ptr_new(/Allocate_Heap)
    res[kLC].value = ptr_new(/Allocate_Heap)
    res[kLC].JD    = ptr_new(/Allocate_Heap)
    res[kLC].flux  = ptr_new(/Allocate_Heap)

    *res[kLC].phase = phase
    *res[kLC].value = f

    *res[kLC].JD    = *lc[kLC].JD
    *res[kLC].flux  = *lc[kLC].flux / mean( *lc[kLC].flux )
    
  endfor

  return, res
end