; docformat = 'rst'
;
; NAME: 
;    yarkovsky
; PURPOSE:
;    Compute the Yarkovsky drift of an asteroid
;+
; :Description:
;    Compute the Yarkovsky drift of an asteroid
;
; :Categories:
;    Asteroid
;
; :Returns: The Yarkovsky drift [in 1e-4 au/Myr]
;
; :Params:
;
; :Keywords:
;    period: in, required, type=float
;      The rotation period of the asteroid (h)
;    diameter: in, required, type=float
;      The diameter of the asteroid (km)
;    albedo: in, required, type=float
;      The geometrical visible albedo of the asteroid
;    inertia: in, required, type=float
;      The thermal inertia of the asteroid (SI)
;    meanMotion: in, required, type=float
;      The asteroid orbital mean motion
;    dSun: in, required, type=float
;      The average distance to the Sun (au)
;    obliquity: in, required, type=float
;      The obliquity of the asteroid (deg)
;    mass: in, required, type=float
;      The mass of the asteroid (kg)
;
;    G: in, optional, type=float, default=0.15
;      The (phase) slope
;    emissivity: in, optional, type=float, default=0.9
;      The emissivity of the surface
;    beaming: in, optional, type=float, default=1.0
;      The beaming parameter
;
; :Examples:
;
; :Uses:
;    
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Written in July 2018, B. Carry (OCA)
;-
function yarkovsky, period=period, diameter=diameter, albedo=albedo, inertia=inertia, $
                    meanMotion=meanMotion, dSun=dSun, obliquity=obliquity, mass=mass, $
                    G=G, beaming=beaming, emissivity=emissivity
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--I.1-- Test on Inputs  ----------------------------------------------------------------------;
  if not keyword_set(period) or $
     not keyword_set(diameter) or $
     not keyword_set(albedo) or $
     not keyword_set(inertia) or $
     not keyword_set(dSun) or $
     not keyword_set(obliquity) or $
     not keyword_set(meanMotion) or $
     not keyword_set(mass) then begin
    message, /ioError, 'Missing parameter!'
    return, 0
  endif

 ;--I.2-- Default Values  ----------------------------------------------------------------------;
  if not keyword_set(G)          then G=0.15
  if not keyword_set(emissivity) then emissivity=0.9
  if not keyword_set(beaming)    then beaming=1.0

 ;--I.3-- Physical Constant  -------------------------------------------------------------------;
  cSpeed  = 299792458.D     ;- m s-1
  au      = 149597870700.D  ;- m
  solCst  = 1360.8          ;- W m-2
  stefBolt= 5.67e-8         ;- W m-2 K-4




;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Computation of Yarkovsky                           -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Subsolar Temperature  ---------------------------------------------------------------;
  bond = albedo*(0.29+0.684*G)
  Tss = ( (1-bond)*solCst / ( dSun*dSun*beaming*stefBolt*emissivity) )^0.25

 ;--II.2-- Rotation Frequency  -----------------------------------------------------------------;
  rotFreq = 2*!PI/(period*3600.)

 ;--II.3-- Diurnal Contribution  ---------------------------------------------------------------;
  ThetaD = inertia* sqrt(rotFreq)    / ( emissivity * stefBolt* Tss^3.)
    WD = -0.5*ThetaD / (1 + ThetaD + ThetaD*ThetaD/2.)

 ;--II.4-- Annual Contribution  ----------------------------------------------------------------;
  ThetaA = inertia* sqrt(meanMotion) / ( emissivity * stefBolt * Tss^3.)
    WA = -0.5*ThetaA / (1 + ThetaA + ThetaA*ThetaA/2.)

 ;--II.4-- Compute Yarkovsky Drift  ------------------------------------------------------------;
  yarko =  ( solCst/ (dSun*dSun) ) * $
           ( (1-bond) / meanMotion ) * $
           !PI*diameter*diameter*(1e6)/ (4* mass * cSpeed) * $
           ( -8.*WD* cos(obliquity*!DTOR)    /9.  + $
              4.*WA*(sin(obliquity*!DTOR)^2.)/9 ) 

 ;--II.5-- Return Yarkovsky Drift  -------------------------------------------------------------;
  return, yarko * (365.2524*86400.) / (1e-10*au)

end