;+
; NAME:
;       DIPOLE_LAMBERT
;
;
; PURPOSE:
;       To plot a given set of celestial coordinates over a
;       CMB dipole simulated map, using lambert projection.
;       A full emisphere is plotted, either in galactic or celestial coordinates
;
;
; CALLING SEQUENCE:
;       DIPOLE_LAMBERT, ra, dec, ngsp, [/GAL_COORDS]
;
; 
; INPUTS:
;       ra:   Right Ascension (hours)
;       dec:  declination (degrees)
;       ngsp: +1 if northern emisphere, -1 if southern 
;
;
; KEYWORD PARAMETERS:
;       GAL_COORD: if set, use galactic coordinate system, else celestial.
;
;
; EXAMPLE:
;       Plot the set of coordinates defined in the vectors ra, dec,
;       over a map of the Southern Galactic Emisphere:
;
;       IDL> DIPOLE_LAMBERT, ra, dec, -1, /gal
;
;       Same thing for the Northern Celestial Emisphere:
;       IDL> DIPOLE_LAMBERT, ra, dec, +1
;
;
; PROCEDURES USED:
;       DIPOLE: Simulates CMB dipole signal
;       DISPLAY: Image display routine
;       EULER: Converts between coordinates (eg. galactic-celestial)
;       LAMBERT_XY_TO_LB: Gives lat-long corresponding to pixel
;          positions in Lambert projection
;       LAMBERT_LB_TO_XY: Gives pixel positions corresponding to
;          lat-long in Lambert projection
;       PRECESS: Precesses equinoxes
;
;
; MODIFICATION HISTORY:
;       Created, Amedeo Balbi, November 1998
;
;-
pro DIPOLE_LAMBERT, rascan0, decscan0, nsgp, GAL_COORD=gal_coord

; defining pixels number:
npix=1024

; converting to degrees:
rascan=rascan0*15.
decscan=decscan0

; defining labels for the coordinates system:
lab1=intarr(12)
lab2=intarr(6)
lab1(*)=1
lab2(*)=1

; creating pixels grid:
xcoord=(lindgen(npix,npix) mod npix - npix/2.)
ycoord=(lindgen(npix,npix) / long(npix) - npix/2.)
xs=lindgen(npix)-npix/2.
ys=lindgen(npix)-npix/2.

if keyword_set(gal_coord) then begin

; galactic coordinates:
;  reading lat-long corresponding to pixels grid:
    print,' ** creating coordinates grid... **'
    LAMBERT_XY_TO_LB,xcoord,ycoord,nsgp,512,l,b
    print,' ** ...done **'
; converting from galactic to celestial coord - then, precess equinox:
    print,' ** converting coordinates... **'
    l1=reform(l,long(npix)*long(npix))
    b1=reform(b,long(npix)*long(npix))
    EULER,l1,b1,ra1,dec1,2
    PRECESS,ra1,dec1,1950,2000
    ra=reform(ra1,npix,npix)
    dec=reform(dec1,npix,npix)
    print,' ** ...done **'    
; reading dipole values and plotting:
    print,' ** calculating dipole signal... **'
    DIPOLE,ra/15.,dec,dipole
    print,' ** ...done **'    
    DISPLAY,dipole,xs,ys
;  plotting labels:
    contour,l,xs,ys,/over,level=360./12.*findgen(12), c_lab=lab1, c_charsize=1.5
    contour,b,xs,ys,/over,level=(180./6.*findgen(6)-90), c_lab=lab2, c_charsize=1.5 
;  converting from celestial to galactic coord - first, precess equinox:
    print,' ** converting coordinates... **'
    PRECESS,rascan,decscan,2000,1950
    EULER,rascan,decscan,lscan,bscan,1
    print,' ** ...done **'
; reading pixel positions corresponding to the set of input coord and overplotting:
    LAMBERT_LB_TO_XY,lscan,bscan,nsgp,512,xscan,yscan
    oplot,xscan,yscan,psy=3

endif else begin

;celestial coordinates:
    print,' ** creating coordinates grid... **'
    LAMBERT_XY_TO_LB,xcoord,ycoord,nsgp,512,ra,dec
    print,' ** ...done **'
    print,' ** calculating dipole signal... **'
    DIPOLE,ra/15.,dec,dipole
    print,' ** ...done **'
    DISPLAY,dipole,xs,ys
    contour,ra,xs,ys,/over,level=360./12.*findgen(12), c_lab=lab1, c_charsize=1.5
    contour,dec,xs,ys,/over,level=(180./6.*findgen(6)-90), c_lab=lab2, $
      c_charsize=1.5
    LAMBERT_LB_TO_XY,rascan,decscan,nsgp,512,xscan,yscan
    oplot,xscan,yscan,psy=3

endelse

return
end