;+
;NAME:
;     HELIOCENTRIC
;PURPOSE:
;     Computes the heliocentric angle given solar latitude and longitude
;CATEGORY:
;CALLING SEQUENCE:
;     result = heliocentric(lonlat,date)
;INPUTS:
;     lonlat = longitude,latitude in degrees, W and N are positive.
;              fltarr(2,N):  lonlat(0,*) = longitude
;                            lonlat(1,*) = latitude
;OPTIONAL INPUT PARAMETERS:
;     date = date in any format.  Defaults to now if date is not
;            specified.
;KEYWORD PARAMETERS
;     dir2center = direction from the heliographic coordinates to the center
;                  of the solar disk as projected on the plane of the sky.  
;                  Degrees CCW from solar north.
;     /terrestrial = dir2center returns degrees CCW from terrestrial north.
;OUTPUTS:
;     result = heliographic angle in degrees: fltarr(N)
;COMMON BLOCKS:
;SIDE EFFECTS:
;RESTRICTIONS:
;PROCEDURE:
;     Calls pb0r to compute the solar b angle
;MODIFICATION HISTORY:
;     T. Metcalf  June 17, 1993
;     20-feb-1997 - S.L.Freeland - use <get_rb0p> in place of "old" pb0r
;-

function heliocentric, lonlat, date, dir2center=dir2center, $
                       terrestrial=terrestrial

   if n_elements(date) LE 0 then date = !stime
   tim = anytim2ints(date)
   b0  =get_rb0p(tim,/b0angle)
   p   =get_rb0p(tim,/pangle)   
   result=reform( acos( sin(b0)*sin(lonlat(1,*)*!dtor) $
                      + cos(lonlat(1,*)*!dtor)*cos(b0)*cos(lonlat(0,*)*!dtor) $
                      )*!radeg $
                )

   dir2center = reform(conv_h2a(lonlat,date))
   dir2center = atan(dir2center(0,*),-dir2center(1,*))*!radeg  ; CCW from N
   if keyword_set(terrestrial) then $
      dir2center = temporary(dir2center) + p*!radeg
   if n_elements(dir2center) EQ 1 then dir2center=dir2center(0)

   if n_elements(result) EQ 1 then return,result(0) $
   else return,result

end