pro terminator, slat, slon, tlat, tlon, npoints, plot=pflag
;
; Given the latitude SLAT and longitude SLON (degrees) of the
; sub-solar point, this procedure returns the latitudes TLAT and
; longitudes TLON (degrees) of npoints evenly-spaced points
; along the day-night terminator (assuming spherical Earth and
; Sun at infinity).
;
; B. Knapp, 2001-10-16
;
; We need a rotation matrix to transform coordinates from the
; "sub-solar point" reference frame to the Earth's x-y-z reference
; frame (z = north, x = 0 degrees longitude).  The "sub-solar
; point" reference frame is the same as the "south-east-zenith"
; reference frame (see Bate, Mueller & White, Fundamentals of
; Astrodynamics, section 6.4).  Our matrix consists of a
; negative rotation of 90-SLAT degrees about the SEZ y-axis,
; followed by a negative rotation of SLON degrees about the
; z-axis.
;
  D2R = !dpi/180.d0
  cl = cos(slon*D2R)
  sl = sin(slon*D2R)
  cb = cos(slat*D2R)
  sb = sin(slat*D2R)
  r = [[sb*cl, sb*sl, -cb], [-sl, cl, 0.d0], [cb*cl, cb*sl, sb]]
;
; How many points?
  if n_elements(npoints) gt 0 then $
     n = long(npoints) $
  else $
     n = 360L
  tlat = dblarr(n)
  tlon = dblarr(n)
;
  dlon =360.d0/n
  txyz = dblarr(6)
  for j=0,n-1 do begin
     jlon = j*dlon*D2R
     tsez = [cos(jlon), sin(jlon), 0.d0]
     txyz[0] = r # tsez
     sphxyz, -1, tsph, txyz
     tlon[j] = tsph[0]
     tlat[j] = tsph[1]
  endfor
;
; Plot?
  if keyword_set(pflag) then begin
     s = sort(tlon)
     get_dev,dev
     plot, tlon(s)/D2R, tlat(s)/D2r, $
       xrange=[0,360], xstyle=1, xticks=12, xminor=6, $
       yrange=[-90,90], ystyle=1, yticks=6, yminor=6
     oplot,[slon,slon],[slat,slat],psym=4
     send_plots,dev
  endif
  return
  end