;+
; NAME:
;   JPLEPHMAKE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Make a new ephemeris suitable for interpolation by JPLEPHINTERP
;
; MAJOR TOPICS:
;   Planetary Orbits, Interpolation
;
; CALLING SEQUENCE:
;   JPLEPHREAD, INFO, RAW, OBJ, T, CX, CY, CZ, $
;               [ POSUNITS=, AUTHOR=, DATE=, OBJECTNAME=, ]
;               [ KEYWORDS=, KEYVALUES=, /RESET ]
;
; DESCRIPTION:
;
;   JPLEPHMAKE is a utility routine which forms an ephemeris table
;   suitable for interpolation with JPLEPHINTERP.  This is a way for
;   users to make or augment an ephemeris of solar system bodies not
;   already present in the JPL planetary ephemeris.  This routine only
;   creates new ephemerides in memory.  No facility is provided to
;   write to disk.
;
;   The user must have already estimated the Chebyshev polynomial
;   coefficients for the body of interest.  One way to do this is with
;   CHEBGRID from the Markwardt library.  
;
;   The two options are either to create a new ephemeris or to augment
;   an existing one.  Augmentation merely means that new columns are
;   added to an existing ephemeris table.  The JPL ephemeris itself
;   can be augmented.
;
;   Even when creating a new ephemeris from scratch, passing an
;   existing INFO structure based on another epehemeris is strongly
;   recommended, because the structure usually contains planetary
;   masses, physical constants, etc. which are relevant.
;
;
;
; PARAMETERS: 
;
;
;   INFO - upon input, an existing INFO structure based on a known
;          ephemeris.  Upon output, a modified INFO structure.
;
;          If INFO is undefined upon input, or the RESET keyword is
;          set, then the returned INFO is set to a generic header.
;
;   RAW - upon input, an existing set of Chebyshev coefficients.  Upon
;         output, the new or augmented set of coefficients.
;
;         If RAW is undefined upon input, or if the RESET keyword is
;         set, then the returned RAW variable is initialized to a new
;         set of keywords.
;
;   OBJ - scalar string, name of the object.
;
;   T - array of times, in Julian Days (TDB), which refer to the
;       *start* epoch of each granule.  [ In the terminology of the
;       JPL ephemeris and CHEBGRID, a "granule" is a single
;       subinterval over which a Chebyshev polynomial is fitted. ] If
;       an existing ephemeris is to be augmented, then T must overlap
;       exactly.
;
;   CX, CY, CZ - arrays of Chebyshev polynomial coefficients.
;
;
;
; KEYWORD PARAMETERS:
;
;   POSUNITS - a scalar string, the units of position as fitted by CX,
;              CY, and CZ.  Allowed values:
;                 'KM' - kilometers  (default)
;                 'CM' - centimeters
;                 'AU' - astronomical units
;                 'LT-S' - light seconds
;
;   NSUBINTERVALS - Number of granules per time sample.
;                   Default: 1
;
;   RESET - if set, then a new ephemeris table is created.  Any
;           Chebyshev coefficients in RAW are overwritten.
;
;   AUTHOR - a scalar string, an identifier giving the author of the
;            new ephemeris.
;            Default: ''
;
;   DATE - a scalar string, the creation date of the ephemeris.
;          Default: SYSTIME(0)
;
;   KEYWORDS - an optional string array, giving any header keywords to
;              be added to the ephemeris (in conjunction with
;              KEYVALUES).
;              Default: (none)
;
;   KEYVALUES - an optional double array, giving any header values for
;               the keywords specified by KEYWORDS.
;
;               Default: (none)
;
;
; EXAMPLE:
;
;
; REFERENCES:
;
;   JPL Export Ephmeris FTP Site
;      ftp://navigator.jpl.nasa.gov/pub/ephem/export/
;      (ephemeris files are available here, however, they must be
;      converted to FITS format using the "bin2eph" utility found in
;      AXBARY)
;
;
; SEE ALSO
;   JPLEPHREAD, JPLEPHINTERP, JPLEPHTEST
;   
; MODIFICATION HISTORY:
;   Written and Documented, CM, Mar 2002
;   Corrected way that ephemerides are merged, also
;     way that AUTHOR field is filled, 29 May 2002, CM
;
;  $Id: jplephmake.pro,v 1.4 2002/05/29 20:07:41 craigm Exp $
;
;-
; Copyright (C) 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy and distribute unmodified copies for
; non-commercial purposes, and to modify and use for personal or
; internal use, is granted.  All other rights are reserved.
;-

pro jplephmake, info, raw, obj, t, cx, cy, cz, reset=reset, $
                posunits=posunits, author=author, date=date, $
                keywords=newkeywords, keyvalues=newkeyvalues, $
                nsubintervals=ns0

  if n_elements(info) GT 0 then begin
      c = info.c
      au = info.au
      keywords = info.keywords
      keyvalues = info.keyvalues
  endif

  if n_elements(author) EQ 0 then author = ''
  if n_elements(date) EQ 0 then date = systime()

  fresh = n_elements(info) EQ 0 OR keyword_set(reset)

  if fresh then begin
      nrows = n_elements(t)-1
      tstart = min(t) & tstop = max(t)
      timedel = t(1)-t(0)
      format= 'JPLEPHMAKE'
      denum = 0L
      if n_elements(c) EQ 0 then $
        c = 299792458d                ;; Speed of light in km/s
      if n_elements(au) EQ 0 then $
        au = 499.00478380613566287138 ;; AU in lt-sec
      jdlimits = [tstart, tstop]
      jdrows = nrows
  endif else begin
      nrows = info.nrows
      tstart = info.tstart & tstop = info.tstop
      timedel = info.timedel
      format = info.format
      denum = info.denum
      jdlimits = info.jdlimits
      jdrows = info.jdrows
      
      objname = info.objname
      ptr = info.ptr
      ncoeff = info.ncoeff
      nsub = info.nsub
      if author EQ '' then author = info.author
  endelse

  ;; Figure the number of subintervals.  If it's specified exactly,
  ;; well okay.  If this is an augmented ephemeris, then it may be
  ;; possible to assign the number of subintervals.
  if n_elements(ns0) EQ 0 then begin
      if fresh then begin
          ns = 1L
      endif else begin
          ns = round(timedel/(t(1)-t(0)))
      endelse
  endif else begin
      ns = round(ns0(0))
  endelse

  sz = size(cx)
  nc = sz(1)
  nr = sz(2) / ns

  xraw = transpose([[[cx]],[[cy]],[[cz]]], [0,2,1])
  xraw = reform(xraw,ns*nc*3L,nr, /overwrite) 

  c_km  = c/1000d
  au_km = au*c_km
  if n_elements(posunits) GT 0 then begin
      case strupcase(strtrim(posunits(0),2)) of
          'KM':   km = 1 ;; Dummy statement
          'CM':   xraw = xraw / 1d5   ;; Convert from CM to KM
          'AU':   xraw = xraw * au_km ;; Convert from AU to KM
          'LT-S': xraw = xraw * c_km  ;; Convert from LT-S to KM
          ELSE: message, 'ERROR: Unrecognized position units'
      endcase
  endif

  if fresh then begin
      objname = [strupcase(strtrim(obj(0),2))]
      ptr = [1L]
      ncoeff = [nc]
      nsub = [ns]

      raw = temporary(xraw)
  endif else begin

      if abs((t(1)-t(0))*ns - timedel) GT 1d-6 then begin
          message, 'ERROR: time sampling does not match'
      endif
      
      if abs(jdlimits(0) - min(t)) GT 1d-6 $
        OR abs(jdlimits(1) - max(t)) GT 1d-6 then begin
          message, 'ERROR: start and stop times do not match'
      endif

      if jdrows NE nr then begin
          message, 'ERROR: number of rows does not match'
      endif

      iprev = n_elements(ncoeff)-1
      objname = [objname, strupcase(strtrim(obj(0),2))]
      ptr = [ptr, max(ptr)+ncoeff(iprev)*nsub(iprev)*3L]
      ncoeff = [ncoeff, nc]
      nsub = [nsub, ns]
      
      szr = size(raw)
      newraw = make_array(szr(1)+nc*3, nr)
      newraw(0,0) = raw
      newraw(szr(1),0) = xraw
      raw = temporary(newraw)
  endelse

  if n_elements(newkeywords) EQ 0 then begin
      newkeywords = ['']
      newkeyvalues = ['']
  endif 
  if n_elements(keywords) EQ 0 then begin
      keywords = newkeywords
      keyvalues = newkeyvalues
  endif else if newkeywords(0) NE '' then begin
      keywords = [newkeywords, keywords]
      keyvalues = [newkeyvalues, keyvalues]
  endif

  auth = author
  info = {filename: '', edited: 1L, $
          creation_date: date, author: auth, $
          nrows: nrows, tstart: tstart, tstop: tstop, $
          timedel: timedel, format: format, denum: denum, $
          c: c, au: au, jdlimits: jdlimits, jdrows: jdrows, $
          objname: objname, ptr: ptr, ncoeff: ncoeff, $
          nsub: nsub, keywords: keywords, keyvalues: keyvalues}

  return
end