;-------------------------------------------------------------
;+
; NAME:
;       q_Tsrn_solpep
; PURPOSE:
;       To calculate Tisserand's parameter. 
; CALLING SEQUENCE:
;       T = q_Tsrn_solpep(q, e, i, j, y, m, d, UT)
; INPUTS:
;       q = semimajor axis of the body
;       e = eccentricity of the body
;       i = inclination of the body
;       j = the j-th planet in the solar system as the main perturber (0<j<9)
;       y = year of the epoch
;       m = month of the epoch
;       d = date of the epoch
;       UT = time of the epoch in decimals
;
; OUTPUTS:
;       T = Tisserand's parameter
;       
; ROUTINES USED: 
;       JDCNV - converts Gregorian dates to Julian days
;
; NOTES:
;       [1] It is a improved version of q_Tsrn.pro with constraints of planets in the
;           solar system. Users do not have to insert the value of semimajor axis of the
;           needed planet as what users have to do with q_Tsrn, but only type the order
;           of the perturber.  
;       [2] Please ALWAYS remember that non-negative j>8 is invalid. 
;       [3] The Tisserand's parameter is derived from one of so called Delaunay standard
;           variables, used to study the perturbed Hamiltonian in 3-body system. Ignoring
;           higher order perturbation terms, the following value is conserved:
;           \sqrt{a (1-e^2)} \cos i
;           Consequently, perturbations may lead to the resonance between the orbit
;           inclination and eccentricity, known as Kozai resonance. Near circular, highly
;           inclined orbits can thus become very eccentric (in exchange for lower 
;           inclination). As example, such mechanism can produce Sun-grazing comets.
;           * TJ, Tisserand’s parameter with respect to Jupiter as perturbing body, 
;             is frequently used to distinguish asteroids (typically TJ > 3) from 
;             Jupiter-family comets (typically 2 < TJ < 3).
;           * The roughly constant value of the parameter before and after the
;             interaction (encounter) is used to determine whether or not an observed
;             orbiting body is the same as a previously observed in Tisserand's Criterion.
;           * The quasi-conservation of Tisserand's parameter constrains the orbits
;             attainable using gravity assist for outer Solar system exploration.
;           * TN, Tisserand's parameter with respect to Neptune, has been suggested to
;             distinguish Near Scattered Objects (believed to be affected by Neptune)
;             from Extended Scattered trans-Neptunian objects (e.g. 90377 Sedna).
;; 
; MODIFICATION HISTORY:
;       M.-T. Hui, 12 Oct, 2011
;
; Copyright (C) 2011, M.-T. Hui
; This software may be used, copied, or redistributed as long as it is not sold and this 
; copyright notice is reproduced on each copy made. This routine is provided as is without
; any express or implied warranties whatsoever. Other limitations apply as described in the
; file disclaimer.txt.
; Converted to IDL V8.0   M.-T. Hui   October 2011
;-
;-------------------------------------------------------------
 
  function q_Tsrn_solpep, q, e, i, j, y, m, d, UT, help=hlp, degrees=degrees
 
  if (n_params(0) lt 8) or keyword_set(hlp) then begin
    print,' To calculate Tisserand parameter of a body.'
    print,' T = q_Tsrn_solp(q, e, i, j, y, m, d, UT)'
    print,'  q = perihelion distance of the body in AU                              in'
    print,'  e = eccentricity of the body                                           in'
    print,'  i = inclination of the body in degrees                                 in'
    print,'  j = the j-th planet in the solar system as the main perturber (0<j<9)  in'
    print,'  y = year of the epoch                                                  in'
    print,'  m = month of the epoch                                                 in'
    print,'  d = date of the epoch                                                  in'
    print,'  UT = time of the epoch in decimal                                      in'
    print,'  T = Tisserand parameter                                               out'
    print,' Keywords:'
    print,'   /DEGREES means angle is in degrees, else radians.'
    return, -1
  endif
 
  cf = 1d
  if keyword_set(degrees) then cf = 180d/!dpi

  JD0 = 2451545.0d
; convert to Julian Date
  JDCNV, y, m, d, UT, JD
  
; days in centuries since 2000.0
  cy = (JD - JD0)/36525d
  
; mean semimajor axis in AU, with respect to the mean ecliptic and equinox of J2000
  a_Me = 0.38709927d
  a_Ve = 0.72333566d
  a_EM = 1.00000261d
  a_Ma = 1.52371034d
  a_Ju = 5.20288700d
  a_Sa = 9.53667594d
  a_Ur = 19.18916464d
  a_Ne = 30.06992276d

  ap0 = [a_Me, a_Ve, a_EM, a_Ma, a_Ju, a_Sa, a_Ur, a_Ne]
  
; the time rate of change of the semimajor axis (AU/century), valid for the time-
; interval 1800 AD - 2050 AD.

  da_Me = 0.00000037d
  da_Ve = 0.00000390d
  da_EM = 0.00000562d
  da_Ma = 0.00001847d
  da_Ju = -0.00011607d
  da_Sa = -0.00125060d
  da_Ur = -0.00196176d
  da_Ne = 0.00026291d
  
  dap = [da_Me, da_Ve, da_EM, da_Ma, da_Ju, da_Sa, da_Ur, da_Ne]
 
  lst = j - 1
  
  if (e NE 1.D) then begin
    a = q/(1 - e)
    T = (ap0[lst] + cy*dap[lst])/a + 2*sqrt(a/(ap0[lst] + cy*dap[lst])*(1 - e^2))*cos(i/cf)
  endif else begin
    T = 2*sqrt(q/(ap0[lst] + cy*dap[lst])*(1 + e))*cos(i/cf)
  endelse
  
  return, T
  end