;+
;DESCRIPTION
;
;  Overplots solar grid lines on an existing plot.  Make a plot first (using
;    /nodata switch in plot command, if no data is to be plotted), then call
;    stonyhurst.
;
;USAGE
;
;  stonyhurst,b_angle,p_angle,R_sun,deg=d_between_lines,dsize=d_between_points,cm=cm
;
;EXAMPLES
;
;1) To plot a solar disk without any other data:
;
;     plot,[-1000,1000],[-1000,1000],xstyle=1,ystyle=1,/nodata
;     stonyhurst,-5.3,19.23,947.0
;
;2) To overplot onto a contour plot with 30 degree grid lines:
;
;     contour,image,x,y
;     stonyhurst,-5.3,19.23,947.0,deg=30
;
;INPUTS
;
;  b_angle        Angle of Sun's rotation axis from plane of the sky in degrees
;
;  p_angle        Angle of Sun's north pole in degrees east of Earth North.
;
;  R_sun          Radius of Sun, in units of plot, as seen from Earth.  If the
;                   plot axes are in arcsec, R_sun should be in arcsec.
;
;KEYWORDS
;
;  deg            Degrees of solar longitude and latitude between grid lines
;                   (default is 10 degrees).
;
;  dsize          Number of degrees between points drawn for each grid line
;                   (default is 5 degrees).  If smooth grid lines are desired
;                   on high-resolution plot devices, but at the cost of speed,
;                   dsize can be made smaller.
;  CM             Angle of central meridian (useful for the Moon).  If omitted,
;                   zero is used.
;
;MODIFICATION HISTORY
;
;  Created by D.E. Gary  3-Oct-1993
;
;  Fixed a problem that occurred when b_angle was exactly 30 or 60.  In this
;  case, iph0 was undefined.  D.E. Gary  6-Oct-1993
;
;  11-Jan-2006  DG
;    Added CM argument, useful for the Moon, which sets the Central Meridian
;    angle.
;-

pro stonyhurst,ba,pa,rsun,deg=deg,dsize=dsize,center=center,linesty=linesty,CM=cm

if (not keyword_set(CM)) then CM = 0
;
; Negative b angles are handled by turning the whole plot upside down (adding
; 180 degrees to pangle).
;
if(ba lt 0) then begin
   b = -ba
   p = -pa+180
endif else begin
   b = ba
   p = -pa
endelse
;
; If b angle is zero, we have problems.  Avoid this by making b small but
; non-zero.
;
if (b eq 0) then b = 0.001
;
; Plot the outer circle (limb) of the Sun
;
oplot,rsun*cos(findgen(361)*!dtor),rsun*sin(findgen(361)*!dtor)
;
; Set keyword defaults
;
if not keyword_set(dsize) then dsize = 5
if not keyword_set(deg) then deg = 10

size1 = 180/dsize
size2 = 360/dsize
b = b*!dtor
;
; Latitude lines
;
x = fltarr(size1+1,size2+1)
y = fltarr(size1+1,size2+1)
maxind = intarr(size2+1)
rlat = (findgen(size1+1)*dsize-90)*!dtor
arg = tan(rlat)*tan(b)
for i=0,size1 do begin
   if(arg(i) ge 1) then begin
      iph0 = 180
   endif else if(arg(i) gt -1) then begin
      iph0 = fix(acos(-arg(i))/!dtor)
   endif
   if(arg(i) gt -1) then begin
      maxind(i) = (iph0*2)/dsize + 1
      index = findgen(maxind(i))
      iphi = index*dsize - iph0
      phi = iphi*!dtor
      x(i,0:maxind(i)-1) = sin(phi)*cos(rlat(i))
      y(i,0:maxind(i)-1) = sin(rlat(i))*cos(b) - cos(phi)*cos(rlat(i))*sin(b)
   endif
endfor
; Rotate x,y for P angle
prad = p*!dtor
xp = x*cos(prad) + y*sin(prad)
yp = -x*sin(prad) + y*cos(prad)
; Plot latitude lines
for i = 0, size1, deg/dsize do begin
   if(maxind(i) gt 0) then begin
      if (keyword_set(center)) then plots,/dev,xp(i,0:maxind(i)-1)*rsun+center[0],$
                                               yp(i,0:maxind(i)-1)*rsun+center[1] $
                            else oplot,xp(i,0:maxind(i)-1)*rsun,yp(i,0:maxind(i)-1)*rsun
   endif
endfor
;
; Longitude lines
;
cos2b = cos(b)*cos(b)
sin2b = sin(b)*sin(b)
rlam = (findgen(size1+1)*dsize-90+CM)*!dtor
th0 = 0.5*atan(2*cos(rlam)*sin(b)*cos(b), $
        cos2b-sin(rlam)^2-(cos(rlam)^2)*sin2b)
th0(where(th0 lt 0))=th0(where(th0 lt 0)) + !pi
for i = 0, size1 do begin
   phi = (findgen(size1+1)*dsize - 90+CM)*!dtor + th0(i)
   x(i,0:size1) = cos(phi)*sin(rlam(i))
   y(i,0:size1) = sin(phi)*cos(b) - cos(phi)*cos(rlam(i))*sin(b)
endfor
; Rotate x,y for P angle
xp = x*cos(prad) + y*sin(prad)
yp = -x*sin(prad) + y*cos(prad)
; Plot longitude lines
for i = 0, size1, deg/dsize do begin
   if (keyword_set(center)) then plots,/dev,xp(i,0:size1)*rsun+center[0],$
                                            yp(i,0:size1)*rsun+center[1] $
                         else oplot,xp(i,0:size1)*rsun,yp(i,0:size1)*rsun
endfor
return
end