; docformat = 'rst'
;
; NAME: 
;    randomLC
; PURPOSE:
;    Randomly creates artificial lightcurves based on Fourier series
;
;+
; :Description:
;    Randomly creates artificial lightcurves based on Fourier series
;
; :Categories:
;    Observations, Simulation
;    
; :Params:
;    t: in, required, type=float
;      The time variable for lightcurve generation
;    seed: in, optional, type=int, default=1
;      The seed for random number generation
;
; :Keywords:
;    pMin: in, optional, type=float, default=2.5
;      Minimum rotation period (same unit as t)
;    pMax: in, optional, type=float, default=10
;      Maximum rotation period (same unit as t)
;    nbP: in, optional, type=int, default=10
;      Number of random period to draw
;    aMin: in, optional, type=float, default=2.5
;      Minimum amplitude (mag)
;    aMax: in, optional, type=float, default=10
;      Maximum amplitude (mag)
;    nbA: in, optional, type=int, default=10
;      Number of random amplitude to draw
;    order: in, optional, type=int, default=2
;      Maximum order of Fourier series
;
;    oPeriod: out, optional, type=flat
;      The periods used to generate each lightcurve
;    oAmplitude: out, optional, type=flat
;      The amplitudes used to generate each lightcurve
;
; :Returns: An array of nbA*nbP lightcurves, evaluated at epochs in t.
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in July 2017, B. Carry (OCA)
;      2017-Aug. - B. Carry (OCA) - Corrected bug on amplitude scaling
;-
function randomLC, t, seed, pMin=pMin, pMax=pMax, nbP=nbP, aMin=aMin, aMax=aMax, nbA=nbA, $
                   order=order, oPeriod=oPeriod, oAmplitude=oAmplitude, pLog=pLog
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--I.1-- Check Syntax  ------------------------------------------------------------------------;
  if not keyword_set(t) then begin
    message, 'syntax: lc = randomLC(t [, seed, pMin, pMax, nbP, aMin, aMax, nbA] )'
    return, -1
  endif


 ;--I.2-- Default Parameters  ------------------------------------------------------------------;
 ;--I.2.1-- Randomness Seed
  if not keyword_set(seed) then seed=1

 ;--I.2.2-- Rotation period
  if not keyword_set(pMin) then pMin = 2.5
  if not keyword_set(pMax) then pMax = 10.
  if not keyword_set(nbP)  then nbP  = 10

 ;--I.2.3-- Lightcurve Amplitude
  if not keyword_set(aMin) then aMin = 0
  if not keyword_set(aMax) then aMax = 1.6
  if not keyword_set(nbA)  then nbA  = 10
  nbLC = nbA * nbP

 ;--I.2.4-- Orders of the Fourier serie
  if not keyword_set(order) then order=2
  coefA = randomu(seed+1, order, nbLC)
  coefB = randomu(seed+2, order, nbLC)
  

 ;--I.3-- Output Variable  ---------------------------------------------------------------------;
  nbT = n_elements(t)
  lcs = fltarr(nbT,nbLC)




;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Create the random Lightcurves                      -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Random Period, Amplitude, Dephasing Arrays  -----------------------------------------;
 ;--II.1.1-- Period
  pArr = pMin + randomu(seed+3,nbLC)*(pMax-pMin)
  if keyword_set(pLog) then begin
    pLogInt = alog10(pMin) + randomu(seed+3,nbLC)*( alog10(pMax)-alog10(pMin))
    pArr = 10.^( pLogInt )
  endif

 ;--II.1.2-- Amplitude
  aArr = aMin + randomu(seed+4,nbLC)*(aMax-aMin)

 ;--II.1.3-- Phase
  phase = randomu(seed+5,nbLC)

 ;--II.1.4-- Export Period and Amplitude
  oPeriod    = pArr
  oAmplitude = aArr


 ;--II.2-- Create Lightcurves  -----------------------------------------------------------------;
  for kLC=0, nbLC-1 do begin

   ;--II.2.1-- Create a Lightcurve from Fourier Expansion
    locPhase = phase[kLC]*pArr[kLC]
    locLC = fourierSerie(t, pArr[kLC], coefA[*,kLC], coefB[*,kLC], x0=locPhase )

   ;--II.2.2-- Create a Lightcurve Convering the Entire Period
    fullT  = pArr[kLC]* findgen(360)/360.
    fullLC = fourierSerie(fullT, pArr[kLC], coefA[*,kLC], coefB[*,kLC], x0=locPhase )

   ;--II.2.3-- Adjust Peak-to-Peak Amplitude 
    amp = max(fullLC) - min(fullLC)
    locLC  *= ( aArr[kLC]/amp )

   ;--II.2.3-- Keep Scaled Lightcurve
    lcs[*,kLC] = locLC

  endfor

 ;--II.3-- Return the Lightcurves  -------------------------------------------------------------;
  return, lcs

end