;-------------------------------------------------------------
;+
; NAME:
;       q_Tsrn
; PURPOSE:
;       To calculate Tisserand's parameter
; CALLING SEQUENCE:
;       T = q_Tsrn(q, e, i, a_p)
; INPUTS:
;       q = semimajor axis of the body
;       e = eccentricity of the body
;       i = inclination of the body
;       a_p = semimajor axis of a perturbing larger body
;
; OUTPUTS:
;       T = Tisserand's parameter
; NOTES:
;       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, 04 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, q, e, i, a_p, help=hlp, degrees=degrees
 
  if (n_params(0) lt 4) or keyword_set(hlp) then begin
    print,' To calculate Tisserand parameter of a body.'
    print,' T = q_Tsrn(q, e, i, a_p)'
    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,'  a_p = semimajor axis of a perturbing larger body in AU           in'
    print,'  T = Tisserand parameter                                         out'
    print,' Keywords:'
    print,'   /DEGREES means angle is in degrees, else radians.'
    return, -1
  endif
 

  cf = 1.D
  if keyword_set(degrees) then cf = 180.D/!dpi
 
  if (e NE 1.D) then begin
    a = q/(1-e)
    T = a_p/a + 2*sqrt(a/a_p*(1-e^2))*cos(i/cf)
  endif else begin
    T = 2*sqrt(q/a_p*(1+e))*cos(i/cf)
  endelse
  
  return, T
  end