; docformat = 'rst'
;+
;This function takes in a Julian Day time and translates it into its
;corresponding Greenwich Sidereal Day time.
;
;NOTES: The constants that are used in this function were taken from "The
;New Definition of Universal Time", McCarthy et al., 1992. 
;
; :Examples:
;   gs_time = gstime(jul_time)
;-
;+
; :Params:
;   jdut1 : in, required, type=double or array of double
;     Julian Date time that is to be converted to Greenwich Sidereal time.
;     May be a scalar or an array.
;     
; :Returns:
;   Returns the Greenwich Sidereal time that is associated with the Julian
;   Date time. 
;-
function gstime, jdut1
  
  if ~real_numeric_type(jdut1) then begin
    doc_library, 'gstime'
    return, -1
  endif

  twopi = 2.0d * !DPI;
  ;deg2rad = !DPI / 180.0d;
  deg2rad = !const.dtor

  tut1 = (jdut1 - 2451545.0d) / 36525.0d;
  a3=-6.2d-6* tut1 * tut1 * tut1;
  a2=+ 0.093104d * tut1 * tut1;
  a1=(876600.0d*3600 + 8640184.812866d) * tut1;
  a0=67310.54841d
  temp=a3+a2+a1+a0;
  temp = -6.2d-6* tut1 * tut1 * tut1 + 0.093104d * tut1 * tut1 + $
           (876600.0d*3600 + 8640184.812866d) * tut1 + 67310.54841d;  // sec
  temp = (temp * deg2rad / 240.0d) mod (twopi); //360/86400 = 1/240, to deg, to rad

  ; ------------------------ check quadrants ---------------------
  w=where(temp lt 0.0d,count)
  if count gt 0 then temp[w] += twopi;
  return,temp
end