;+
; NAME:
;   GEOGREAD
;
; 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:
;   Read gravity model from disk file
;
; MAJOR TOPICS:
;   Physics, Gravity, Geodesy, Spacecraft Navigation
;
; CALLING SEQUENCE:
;   GEOGREAD, ROOTFILE, MODEL [, STATUS=STATUS, ERRMSG=ERRMSG,
;                               /COEFF_ERR]
;
; DESCRIPTION:
;
;   GEOGREAD reads a gravity model from a disk file.  The gravity
;   model must have already been prepared.  There are a number of
;   freely available models.
;
;   Each model must have a "description" file which describes, in IDL
;   syntax, the name, content and format of the model file.  The
;   ROOTFILE parameter is the name of this description file.  The
;   description must provide the name of the model file (the FILENAME
;   field), which must reside in the same directory.
;
; FILE FORMAT:
;
;   The format of the description file (and hence also the format of
;   the structure returned in GEOGMOD), is as follows, an example
;   modified from egm96.desc:
;
; { $
;   name: 'EGM96', $              ;; Title of the model
;   type: 'GRAVITY', $            ;; Type of model 'GRAVITY' or 'BFIELD'
;   filename: 'EGM96.GEO', $      ;; Model coefficient file name (same dir)
;   reference: 'Lemoine, ...'     ;; Complete literature reference
;   url: 'ftp://ftp.csr.utexas.edu/pub/grav/EGM96.GEO.Z', $   ;; Source URL
;   nmax: 360L, $                 ;; Maximum order (inclusive)
;   mmax: 360L, $                 ;; Maximum degree (inclusive)
;   normalized: 1, $              ;; Coefficients are normalized (1=yes, 0=no)
;   mu: 398600.44150D+09, $       ;; GM for central body [m^3/s^2]
;   a: 6378136.30d, $             ;; Mean equatorial radius [m]
;   tide: 'ZERO', $               ;; Tide system (ZERO, FREE, or MEAN)
;   epoch: 1986.0d, $             ;; Epoch of model coefficients (Julian year)
;   C21: -.1869876359548955D-09,$ ;; C21 coefficient (if not in Cnm table)
;   S21:  .1195280120306540D-08,$ ;; S21 coefficient (if not in Cnm table)
;   C20_dot: 1.16275534D-11,$     ;; C20 rate (unitless; yr^-1)
;   C21_dot: -0.32d-11, $         ;; C21 rate (unitless; yr^-1)
;   S21_dot: +1.62d-11, $         ;; S21 rate (unitless; yr^-1)
;   rowstart: 4L, $               ;; Coefficient starting row (first row = 0)
;   nrows: 65338L, $              ;; Number of coefficient rows in file
;   ncolrange: [6,8], $           ;; Column range for degree (first col = 0)
;   mcolrange: [9,11], $          ;;    "     "    "  order 
;   Ccolrange: [12,30], $         ;;    "     "    "  C coefficients
;   Scolrange: [31,49], $         ;;    "     "    "  S coefficients
;   dCcolrange: [50,62], $        ;;    "     "    "  C std deviation
;   dScolrange: [63,75] $         ;;    "     "    "  S std deviation
; }
;
;   The xCOLRANGE fields describe which character columns in the model
;   file, inclusive, contain the quantity of interest.  You can use a
;   text editor which reports the column number to find these values.
;   Exclude any character columns that contain field delimiters such
;   as commas.
;
;   Since the C21 and S21 coefficients are commonly not included in
;   the table itself, their values are allowed to be specified in the
;   description file.  If the coefficients *are* in the table, then
;   they must be set to zero in the description file to avoid double
;   computations.  The coefficient rates can be used to extrapolate to
;   different epochs from the reference epoch (specified by EPOCH).
; 
;
; INPUTS:
;
;   ROOTFILE - scalar string, the name of the model description file.
;
;   GEOGMOD - upon return, an IDL structure containing the model
;             information.  In addition to the fields listed above,
;             other fields are appended which contain (pointers to)
;             the coefficient data, etc.
;
; KEYWORD PARAMETERS:
;
;   STATUS - upon return, a status indicator.  A value of 1 is OK, 0
;            indicates an error condition.
;
;   ERRMSG - upon return, an error message, if any.  If no error
;            occurred, then ERRMSG is set to ''.
;
;   CEOFF_ERR - if set, then coefficient standard deviations are also
;               read in.q
;
;
; EXAMPLE:
;   GEOGREAD, 'egm96', egm96
;   GEOGRAV, egm96, r, phi, a
;
;   Read the gravity model "EGM96" and evaluate it at position "R" in
;   body coordinates.  The potential and acceleration are returned in
;   PHI and A.
;
; MODIFICATION HISTORY:
;   Documentation additions, CM, 26 Sep 2004
;
;  $Id: geogread.pro,v 1.3 2004/09/26 14:58:19 craigm Exp $
;
;-
; Copyright (C) 2004, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
pro geogread, rootfile, geogmod, status=status, errmsg=errmsg, $
              coeff_err=cerr

  status = 0
  errmsg = ''

  if n_params() EQ 0 then begin
      message, 'USAGE: GEOGREAD, ROOTFILE, GEOGMOD [, STATUS=, ERRMSG=]', /info
      return
  endif

  openr, unit, rootfile, /get_lun, error=err
  if err NE 0 then begin
      errmsg = 'ERROR: could not open '+rootfile
      return
  endif

  rootstr = ''
  on_ioerror, DONE_ROOTFILE
  while NOT eof(unit) do begin
      s = ''
      readf, unit, s
      s = strtrim(s,2)
      if s NE '' then begin
          p = strpos(s, '$')
          if p GE 0 then s = strtrim(strmid(s,0,p),2)
          p = strpos(s, ';')
          if p GE 0 then s = strtrim(strmid(s,0,p),2)

          if s NE '' then rootstr = rootstr + s
      endif
  endwhile
  DONE_ROOTFILE:
  free_lun, unit

  if rootstr EQ '' then begin
      errmsg = 'ERROR: no data found in '+rootfile
      return
  endif

  temp = strlowcase(rootstr)
  if strpos(temp,'execute') GE 0 OR strpos(temp,'call_') GE 0 $
    OR strpos(temp,'spawn') GE 0 OR strpos(temp,'openw') GE 0 then begin
      errmsg = 'ERROR: '+rootfile+' contains insecure commands'
      return
  endif
  
  cmdstr = 'geogmod = '+rootstr
  dummy = execute(cmdstr)
  if dummy NE 1 then begin
      errmsg = 'ERROR: could not parse structure in '+rootfile
      return
  endif
  sz = size(geogmod) 
  if sz(sz(0)+1) NE 8 then begin
      errmsg = 'ERROR: structure not found in '+rootfile
      return
  endif

  tn = strlowcase(tag_names(geogmod))
  if max(tn EQ 'filename') EQ 0 OR max(tn EQ 'type') EQ 0 OR $
    max(tn EQ 'a') EQ 0 OR max(tn EQ 'mu') EQ 0 OR $
    max(tn EQ 'nmax') EQ 0 OR max(tn EQ 'mmax') EQ 0 OR $
    max(tn EQ 'nrows') EQ 0 OR max(tn EQ 'rowstart') EQ 0 OR $
    max(tn EQ 'ncolrange') EQ 0 OR max(tn EQ 'mcolrange') EQ 0 OR $
    max(tn EQ 'ccolrange') EQ 0 OR max(tn EQ 'scolrange') EQ 0 then begin
      errmsg = 'ERROR: structure did not contain required fields in '+rootfile
      return
  endif

  ps = path_sep()

  p = rstrpos(rootfile, ps)
  if p GE 0 then path = strmid(rootfile,0,p+1) $
  else           path = ''

  file = path + geogmod.filename

  openr, unit, file, /get_lun, error=err
  if err NE 0 then begin
      errmsg = 'ERROR: could not open '+file
      return
  endif

  nrows = geogmod.nrows

  reading = 1
  on_ioerror, DONE_READING
  str = strarr(geogmod.rowstart)
  readf, unit, str

  str = strarr(nrows)
  readf, unit, str
  reading = 0
  
  DONE_READING:
  free_lun, unit
  if reading then begin
      errmsg = 'ERROR: could not read coefficient rows from '+file
      return
  endif

  rng = geogmod.ncolrange & n = long(strmid(str,rng(0),rng(1)-rng(0)+1))
  rng = geogmod.mcolrange & m = long(strmid(str,rng(0),rng(1)-rng(0)+1))

  rng = geogmod.Ccolrange & C = double(strmid(str,rng(0), rng(1)-rng(0)+1))
  rng = geogmod.Scolrange & S = double(strmid(str,rng(0), rng(1)-rng(0)+1))

  if keyword_set(cerr) then begin
      rng = geogmod.dCcolrange & dC = double(strmid(str,rng(0), rng(1)-rng(0)+1))
      rng = geogmod.dScolrange & dS = double(strmid(str,rng(0), rng(1)-rng(0)+1))
  endif
  str = 0
  
  Cnm = dblarr(geogmod.nmax+1,geogmod.mmax+1)
  Snm = Cnm

  if keyword_set(cerr) then begin
      dCnm = Cnm
      dSnm = Cnm
      dCnm(n,m) = temporary(dC) & dSnm(n,m) = temporary(dS)
  endif

  Cnm(0,0) = 1
  Cnm(n,m) = temporary(C)  &  Snm(n,m) = temporary(S)

  geogmod = create_struct(geogmod, $
                          'Cnm', ptr_new(Cnm, /no_copy), $
                          'Snm', ptr_new(Snm, /no_copy))
  if keyword_set(cerr) then begin
      geogmod = create_struct(geogmod, $
                              'dCnm', ptr_new(dCnm, /no_copy), $
                              'dSnm', ptr_new(dSnm, /no_copy))
  endif

  return
end