pro coord_sphere,rsun,p,p0,l0,b0,grid
;+
; Project     : Solar Physics 
;
; Name        : COORD_SPHERE 
;
; Category    : Graphics 
;
; Explanation : This routine overplots a spherical coordinate grid onto contourplot of sun
;		where X-axis and Y-axis are given in degrees relative to the suncenter
;
; Syntax      : IDL> coord_sphere,rsun,p,p0,l0,b0,grid
;
; Inputs      : 
;		RSUN	=radius of optical sun (in units of used coordinate system)
;		P	=Position angle of solar North (from map North to East)
;		P0	=rotation angle of image (e.g. CROTA2=-POS if image previously
;	 		 has already been oriented to solar North, otherwise CROTA2=0)
;		B0	=heliographic latitude of solar disk center
;		L0      =heliographic longitude of solar disk center
;		GRID 	=spacing of spherical grid in degrees (e.g. 10 or 15 deg)
;			 The parameters POS,BLAT,BLONG can be found in "Nautical Almanach" and
;			 are given in degrees.
;
; Outputs     : 
;
; History     : March 1998, Version 1 written by Markus J. Aschwanden, UMd 
;
; Contact     : aschwanden@lmsal.com
;-

!noeras	=1
!p.linestyle=0
phi	=2.*!pi*findgen(361)/360.
oplot,rsun*cos(phi),rsun*sin(phi)

!p.linestyle=1
x0	=0.
y0	=0.
ilat	=fix(89./grid)
nlat	=2.*ilat+1
lat	=grid*(-ilat+findgen(nlat))
l	=findgen(361)
r	=(l*0+1.)*rsun
for i=0,nlat-1 do begin
 b	=l*0.+lat(i)
 helio_trans2,p,p0,l0,b0,l,b,r,x,y,z,rsun
 ind	=where(z gt 0,nind)
 if (nind ge 2) then begin
  igap	=where((ind(1:nind-1)-ind(0:nind-2)) ne 1,ngap)
  if (ngap ge 1) then ind=[ind(igap(0)+1:nind-1),ind(0:igap(0))]
  oplot,x(ind),y(ind)
 endif
endfor

nlon	=fix(360./grid)
lon	=grid*findgen(nlon)
b	=findgen(361)
for i=0,nlon-1 do begin
 l	=b*0.+lon(i)
 helio_trans2,p,p0,l0,b0,l,b,r,x,y,z,rsun
 ind	=where(z gt 0,nind)
 if (nind ge 2) then begin
  igap  =where((ind(1:nind-1)-ind(0:nind-2)) ne 1,ngap)
  if (ngap ge 1) then ind=[ind(igap(0)+1:nind-1),ind(0:igap(0))]
  oplot,x(ind),y(ind)
 endif
endfor
!p.linestyle=0
end