;+
; NAME:
;	SOLARGRID
; Purpose:
;	The SOLARGRID procedure draws the graticule of parallels and meridians;
;	 
; CALLING SEQUENCE:
;	solargrid,b0,dx=dx, dy=dy, p0=p0, rsun=rsun,$
;		labels=[long,lat], /nolong, /nolat,/nsew,$
;		color=color, range=range, linestyle=linestyle, cx=cx,cy=cy
; INPUTS:
;	b0 = solar B0 or date string (for ex. '1999-05-25')
;
; OPTIONAL INPUT KEYWORDS:
;	p0 - solar position angle
;	rsun - solar radius, default=960 arcsec
;	dx - longitude increment, default=10 degree
;	dy - latitude increment, default=10 degree
;	cx,cy - solar center, default 0,0
;	range - grid range [x0,y0,x1,y1]
;	labels - draw grid labels
;	/nsew - print labels with N,S,E,W marks
;	/nolong - draw only latitude lines
;	/nolat - draw only longitude lines
;
;
; MODIFICATION HISTORY:
; 	Written by:	Vladimir Garaimov, May 2002
;-

pro sunadtoxy,a,d,b0,p0,x,y,z=z

  a1=a*!dtor
  d1=d*!dtor

  y=sin(d1)
  x=sin(a1)*cos(d1)
  z=cos(d1)*cos(a1)

  if b0 ne 0 then begin
	b=-b0*!dtor
	z1=z & y1=y
	z=z1*cos(b)-y1*sin(b)
	y=z1*sin(b)+y1*cos(b)
	end

  i=where(z gt 0)
  if i(0) ne (-1) then begin
	x=x(i) & y=y(i)
  endif else begin
	x=-2 & y=-2 &return
  end

  if p0 ne 0 then begin
	x1=x &y1=y
	p=-p0*!dtor
	x=x1*cos(p)-y1*sin(p)
	y=x1*sin(p)+y1*cos(p)
	end
end

pro solargrid,_b0,dx=dx,dy=dy,labels=labels,color=color,$
	p0=p0,rsun=rsun,linestyle=linestyle,range=range,nsew=nsew,$
	cx=cx,cy=cy, nolat=nolat, nolong=nolong

if n_params(0) eq 0 then begin
	doc_library,'solargrid'
	return
end
b0=_b0
sx=size(b0)
sx=sx(n_elements(sx)-2)
if sx eq 7 then begin
 b0=(pb0r(b0(0)))(1)
 print,'b0=',b0
end
if n_elements(rsun) eq 0 then _rsun=960. else _rsun=float(rsun(0))
if n_elements(dx) eq 0 then _xstep=10. else _xstep=float(dx(0))
if n_elements(dy) eq 0 then _ystep=10. else _ystep=float(dy(0))
if n_elements(p0) ne 0 then _p0=float(p0(0)) else _p0=0.
if n_elements(cx) eq 0 then sx=0. else sx=float(cx(0))
if n_elements(cy) eq 0 then sy=0. else sy=float(cy(0))

if n_elements(color) eq 0 then col=!p.color else col=color(0)
if n_elements(linestyle) eq 0 then ls=1 else ls=linestyle(0)

if n_elements(range) ne 4 then begin
	x0=-180. & x1=180. & y0=-90. & y1=90. &ks=1.
endif else begin
	ks =0.
	range=float(range)
	x0=range(0)>(-180) & x1=range(2)<180
	y0=range(1)>(-90) & y1=range(3)<90
	if x0 ge x1 then begin &x0=-180.&x1=180.& end
	if y0 ge y1 then begin &y0=-90.& y1=90.& end
endelse

xl=500 & yl=500
if keyword_set(labels) then begin
  if n_elements(labels) eq 2 then begin
	if abs(labels(0)) lt 90 then xl=labels(0)
	if abs(labels(1)) lt 90 then yl=labels(1)
  endif else begin
	xl=(x1+x0)/2
	yl=(y1+y0)/2
  endelse
end 

nn=fix((x1-x0)/_xstep)+1
mm=fix((y1-y0)/_ystep)+1

if keyword_set(nolat) then goto, xlat

yy=fltarr(nn)
for i=y0+ks*_ystep,y1-_ystep*ks,_ystep do begin
 xx=findgen(nn)*_xstep+x0
 yy(*)=float(i)
 sunadtoxy,xx,yy,b0,_p0,zx,zy 
 if zx(0) ne -2 then begin
  zx=zx*_rsun &zy=zy*_rsun
  plots,zx+sx,zy+sy,color=col,lines=ls,noclip=0
 end
 if xl ne 500 then begin
   sunadtoxy,xl+_xstep/3.,i,b0,_p0,zx,zy 
   if zx(0) ne -2 then begin
     zx=zx*_rsun &zy=zy*_rsun
     if keyword_set(nsew) then begin
	j=strtrim(string(fix(abs(i))),2)
	case 1 of
	 i lt 0: j='S'+j
	 i gt 0: j='N'+j
	 else:
	endcase
     endif else j=strtrim(string(fix(i)),2)
     if (zx(0) le !x.crange(0)) or (zx(0) ge !x.crange(1)) then goto, xl_end
     if (zy(0) le !y.crange(0)) or (zy(0) ge !y.crange(1)) then goto, xl_end
     xyouts,zx,zy,j,col=col,noclip=0
   end 
   xl_end: 
 end
end

xlat: if keyword_set(nolong) then goto,xend

xx=fltarr(mm)
for i=x0+ks*_xstep,x1-_xstep*ks,_xstep do begin
 yy=findgen(mm)*_ystep+y0
 xx(*)=float(i)
 sunadtoxy,xx,yy,b0,_p0,zx,zy 
 if zx(0) ne -2 then begin
  zx=zx*_rsun &zy=zy*_rsun
  plots,zx+sx,zy+sy,color=col,lines=ls,noclip=0
 end
 if yl ne 500 then begin
   sunadtoxy,i,yl+_ystep/4.,b0,_p0,zx,zy 
   if zx(0) ne -2 then begin
     zx=zx*_rsun &zy=zy*_rsun
     if keyword_set(nsew) then begin
	j=strtrim(string(fix(abs(i))),2)
	case 1 of
	 i lt 0: j='E'+j
	 i gt 0: j='W'+j
	 else:
	endcase
     endif else j=strtrim(string(fix(i)),2)
     if (zx(0) le !x.crange(0)) or (zx(0) ge !x.crange(1)) then goto, yl_end
     if (zy(0) le !y.crange(0)) or (zy(0) ge !y.crange(1)) then goto, yl_end
     xyouts,zx,zy,j,col=col,noclip=0
   end 
   yl_end: 
 end

end

xend: empty
end