;+
; NAME:
; oplot_dipole
;
;
; PURPOSE:
; overplot dipole field lines on a plot
;
;
; CATEGORY:
; plotting
;
;
; INPUTS:
; L: array of L values to plot
;
;
;
; KEYWORD PARAMETERS:
; NPTS: number of points to plot on each field line
; RIGHT_ONLY: only plot the right side
; LEFT_ONLY: only plot the left side
; _EXTRA: extra keywords to oplot
;
;
; OUTPUTS:
; oplot on the current graphics window
;
;
; RESTRICTIONS:
; only plots in equitorial plane
;
;
; EXAMPLE:
; plot, findgen(20)-10, findgen(20)-10, /nodata, /iso
; oplot_dipole, [2,3,4,5,6,7,8]
;
;
; MODIFICATION HISTORY:
;
;       Mon Apr 28 11:54:45 2008, Brian Larsen
;
;		added _extra to plot_earth call
;
;       Mon Jul 23 15:05:58 2007, Brian Larsen
;
;		written and tested
;
;-
PRO oplot_dipole, L, $
                  NPTS=npts, $
                  RIGHT_ONLY=ro, $
                  LEFT_ONLY=lo, $
                  _EXTRA=_extra

if (N_ELEMENTS(npts) eq 0) then begin
 npts = 400.
endif

er = make_array(npts, /float, value=1)

thr = (!dpi)*(findgen(npts+1)/npts-0.5)
thl = (!dpi)*(findgen(npts+1)/npts)+!dpi/2.


th = findgen(npts+1)*2.*!dpi/npts-!dpi/2.


plot_earth, _STRICT_EXTRA=_extra ; fix this so there is no line...

if (n_elements(L) eq 0 ) then L= 6.6

for i=0l, n_elements(L)-1 do begin
    if (not keyword_set(lo)) then begin
        ind = where(L[i]*sin(thr+!dpi/2.)^2 gt 1, ind_num)
        if ind_num gt 0 then oplot, L[i]*sin(thr[ind]+!dpi/2.)^2, thr[ind], /polar, _STRICT_EXTRA=_extra
    ENDIF

    if (not KEYWORD_SET(ro)) then BEGIN
        ind = where(L[i]*sin(thl+!dpi/2.)^2 gt 1, ind_num)
        if ind_num gt 0 then oplot, L[i]*sin(thl[ind]+!dpi/2.)^2, thl[ind], /polar, _STRICT_EXTRA=_extra
    ENDIF
ENDFOR



end