;+
; NAME:
;       DUST_GNOM
;
;
; PURPOSE:
;       To plot a given set of celestial coordinates over the
;       IRAS/DIRBE 100 micron dust maps, using gnomic projection.
;       The area of the map is fixed to 600x600 pixels, each
;       pixel side being 10 arcmin, so that the region is 10x10 square
;       degrees wide.
;
;
; CALLING SEQUENCE:
;       DUST_GNOM, racenter, decenter, side, ra, dec
;
; 
; INPUTS:
;       racenter:  Right Ascension of the center of the map (degrees)
;       deccenter: declination of the center of the map (degrees)
;       side:      angular size of the side of the map (degrees)
;       ra:        Right Ascension (hours)
;       dec:       declination (degrees)
;
;
; EXAMPLE:
;       Plot the set of coordinates defined in the vectors ra, dec,
;       over a 100x100 degs map, taking the North Celestial Pole
;       (90,90) as center:
;
;       IDL> DUST_GNOM, 90, 90, 100, ra, dec
;
;
; PROCEDURES USED:
;       DISPLAY: Image display routine
;       EULER: Converts between coordinates (eg. galactic-celestial)
;       ISSA_PROJ_GNOM: Gives pixel positions corresponding to
;          lat-long in gnomic projection
;       ISSA_PROJ_IGNOM: Gives lat-long corresponding to pixel
;          positions in gnomic projection
;       PRECESS: Precesses equinoxes
;
;
; MODIFICATION HISTORY:
;       Created, Amedeo Balbi, November 1998
;
;-
pro DUST_GNOM, racenter0, deccenter0, side, rascan0, decscan0

; defining image size:
pxsize=10              ; pixel side in arcminutes
scale=3437.7468/pxsize ; scale (needed in issa routines)

xsize=fix(side*60./pxsize)  ; map dimensions
ysize=fix(side*60./pxsize)

; issa routines work in radians, converting everything:
rascan  = rascan0*15.*!dtor
decscan = decscan0*!dtor
racenter  = racenter0*!dtor
deccenter = deccenter0*!dtor

; creating pixels grid:
xbox = lindgen(xsize, ysize) MOD long(xsize) - xsize/2.
ybox = lindgen(xsize, ysize) / long(xsize) - ysize/2.
xs = lindgen(xsize) - xsize/2.
ys = lindgen(ysize) - ysize/2.

; reading lat-long corresponding to pixels grid:
print, ' ** building coordinates grid... ** '
ISSA_PROJ_IGNOM, xbox, ybox, racenter, deccenter, scale, ra, dec
print, ' ** ...done **'

; back to degrees:
ra  = ra*!radeg
dec = dec*!radeg

; allocating grid in vectors (needed in order to use PRECESS and EULER):
ra1  = reform(ra, long(xsize)*long(ysize))
dec1 = reform(dec, long(xsize)*long(ysize))

; converting from celestial to galactic coord - first, precess
;                                               equinox:
print, ' ** converting coordinates... **'
PRECESS, ra1, dec1, 2000., 1950.
EULER, ra1, dec1, l1, b1, 1
print, ' ** ...done **'

; recreating grid of galactic coordinates:
l = reform(l1, xsize, ysize)
b = reform(b1, xsize, ysize)

; reading dust values:
print,' ** reading dust intensity values... **'
image = DUST_GETVAL(l, b, map='I100', /noloop, /verbose)
print,' ** ...done **'

; displaying map:
DISPLAY,image,xs,ys,min=-3,max=80

; plotting coordinates labels:
lab1=intarr(12)
lab2=intarr(6)
lab1(*)=1
lab2(*)=1
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

; reading pixel positions corresponding to the set of input coord and overplotting:
ISSA_PROJ_GNOM,rascan,decscan,racenter,deccenter,scale,xrascan,ydecscan
oplot,xrascan,ydecscan
print, 'ra'
print, xrascan[0:5]
print, 'dec'
print, ydecscan[0:5]

return
end