;+
; NAME:
;	smei_sky_read
; PURPOSE:
;	Reads orbital skymaps
; CATEGORY:
;	camera/idl/sky
; CALLING SEQUENCE:
	FUNCTION smei_sky_read, wanted_map,	$
		source		= source	, $
		nbin		= nbin		, $
		silent		= silent	, $

		xgrid		= xgrid 	, $
		ygrid		= ygrid 	, $
		xedge		= xedge 	, $
		yedge		= yedge 	, $

		degrees 	= degrees	, $

		hdrs		= hdrs		, $
		fixgain 	= fixgain	, $
		rmbkgnd		= rmbkgnd 	, $
		rmzld		= rmzld 	, $
		exten_no	= exten_no	, $

		use_flat	= use_flat	, $

		num			= num		, $
		_extra		= _extra
; INPUTS:
;	wanted_map 		scalar, array[1]; type: numerical, time structure, or string
;						indicates the orbit(s) to be processed
;						(passed to smei_getfile together with _extra).
; OPTIONAL INPUT PARAMETERS:
;	source=source	scalar; type: string
;						source directory for skymaps
;	(the following keywords are processed by href=smei_sky_field=)
;
;	/equatraw
;	/northraw
;	/southraw
;	/equatflat
;	/northflat
;	/southflat
;	/equatcnt
;	/northcnt
;	/southcnt
;	/equatstars
;	/northstars
;	/southstars
;	/dirtysky
;	/orbtime
;	/orbsecs
;	/psfn
;	/psfe
;	/badsky
;	/badhtm
;	/fovx
;	/thetax
;	/thetay
;	/hit
;	/badpixels
;					if none of these map types is selected then the hires maps
;					(equatorial,south and north pole) are put together into a
;					single equatorial map covering [-90,+90] degrees declination.
;					Note that this option also handles skymaps that contains
;					a single equatorial map covering the full latitude range.
;
;	nbin			scalar; type: integer; default: 5
;						rebinning factor applied to high-resolution maps.
;						A rebinning with factor 5 results in maps
;						with 5 deg bins in RA and Dec.
;
;	/fixgain		if set then sky brightnesses are corrected for changing
;					camera gain (see href=smei_camera_gain=)
;
;	/use_flat		if set for 'equ' and 'ecl' files then the
;					maps with subtracted star replaced by a flat background
;					are used
; OUTPUTS:
;	Result			array; type: float
;						skymap (-1 if no skymap was read)
; OPTIONAL OUTPUT PARAMETERS:
;	orbit			scalar; type: integer
;						actual orbit number used
;	xgrid			horizontal grid (bin centers)
;	ygrid			vertical grid (bin centers)
;	xedge			horizontal grid (bin edges)
;	yedge			vertical grid (bin edges)
;
;	num=num 		array; type: integer
;						array of same size as returned skymap
;						Contains 1,2 or 3 identifying camera
;						that observed correspondig skybin
;						(contains 0 where skymap is empty).
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsTime, IsType, destroyvar, boost
;	gridgen, readfits, headfits, fxpar, TimeSet, hide_env
;	ToRadians, ToDegrees, smei_getfile, smei_filename, BadValue
;	smei_camera_gain
; EXTERNAL:
;	smei_sky_hdr__define
; PROCEDURE:
;	Way too complicated.
; MODIFICATION HISTORY:
;	DEC-2004, Paul Hick (UCSD/CASS)
;	JUL-2006, Paul Hick (UCSD/CASS)
;		Added /fovx to extract Fits extension containing direction cosine
;		angle in long dimension. If the extension is not present then an
;		array with value of 15 degree everywhere is returned.
;	DEC-2007, Paul Hick (UCSD/CASS)
;		Added keyword /fixgain
;	APR-2008, Paul Hick (UCSD/CASS)
;		Added keyword /rmbkgnd and /rmzld to be able to remove
;		sidereal background and zodiacal light distribution.
;	MAY-2008, Paul Hick (UCSD/CASS)
;		Bugfixes in ZLD removal (missing degrees=degrees).
;	JUN-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added "num" keyword
;-
InitVar, nbin		, 5
InitVar, silent 	, 0

InitVar, fixgain	, 0
InitVar, rmbkgnd	, /key
InitVar, rmzld		, /key

InitVar, use_flat	, /key

rpm = ToRadians(degrees=degrees)
dpm = ToDegrees(degrees=degrees)

given_map = smei_getfile(wanted_map, $
	source	= source	, $
	count	= count 	, $
	_extra	= _extra	)

IF count EQ 0 THEN	$
	RETURN, -1

skytime = smei_filename(given_map[0], mode=mode, ccd=ccd)

IF ccd THEN BEGIN
	message, /info, given_map[0]
	message, 'you seem to be processing a CCD frame !!!'
ENDIF

; Find the extension to be read

InitVar, exten_no, smei_sky_field(given_map[0],_extra=_extra)

hires_map = smei_sky_field(given_map[0], exten_no=exten_no, /hires_map )
lores_map = smei_sky_field(given_map[0], exten_no=exten_no, /lores_map )
eclsun_map= smei_sky_field(given_map[0], exten_no=exten_no, /eclsun_map)
frame_map = smei_sky_field(given_map[0], exten_no=exten_no, /frame_map )
adus_map  = smei_sky_field(given_map[0], exten_no=exten_no, /adus_map  )
angle_map = smei_sky_field(given_map[0], exten_no=exten_no, /angle_map )

IF rmbkgnd THEN BEGIN
	tmp = filepath(root=smei_filepath(mode='sid'),'c0sid*.fts.gz')
	bkgnd_map = (file_search(tmp))[0]
	rmbkgnd = bkgnd_map NE ''
	IF NOT rmbkgnd THEN BEGIN
		message, /info, 'no sidereal background map, '+hide_env(tmp)
		message, /info, 'background subtraction skipped'
	ENDIF
ENDIF

CASE 1 OF

hires_map:

lores_map:

frame_map:

eclsun_map:

ELSE: BEGIN 						; Construct single map covering RA=[0,360], dec=[-90,+90]

	nra_bins  = nbin				; # input bins per output bin
	ndec_bins = nbin

	; Get information about the input map from the primary header in the
	; first file to be read

	; The first extension should refer to an equatorial map covering
	; dec=[-60,+60] (usually a 3600x1200 array, or it already is a complete map
	; covering [-90,+90] (usually a 3600x1800 array)

	skyfile = given_map[0]
									; Read primary header
	hdr = headfits(skyfile, silent=silent+1)

	; We assume that the equatorial input map has 10 bins per degree

	d_eq = fxpar(hdr,'CDELT1')
	IF !err GE 0 THEN d_eq = round(1.0d0/d_eq) ELSE d_eq = fxpar(hdr,'D_EQ')
									; Bins per degrees in input map
	dra_in  = 1.0d0/d_eq 			; Input bin width in degrees for RA
	ddec_in = dra_in				; Input bin width in degrees for dec

	nra_in  = round(360*d_eq)		; # Input RA bins needed to cover 360 degrees
	ndec_in = round(180*d_eq)		; # Input dec bins needed to cover 180 degrees

	ndim = fxpar(hdr, 'NAXIS*') 	; # Input dec bins covered by equatorial map in input files
	mdec_in = ndim[1] > 1200		; Will be 1800 (same as ndec_in) or 1200

	nra_out  = nra_in /nra_bins 	; # Output RA bins needed to cover 360 degrees
	ndec_out = ndec_in/ndec_bins	; # Output dec bins needed to cover 180 degrees
	mdec_out = mdec_in/ndec_bins	; # Output dec bins for dec range covered by input file

	; Sanity checks (probably redundant)

	IF nra_in/nra_out*nra_out NE nra_in THEN BEGIN
		message, /info, 'illegal nra; must be a factor of'+strcompress(nra_in)
		RETURN, -1
	ENDIF

	IF ndec_in/ndec_out*ndec_out NE ndec_in OR	$
	   mdec_in/mdec_out*mdec_out NE mdec_in THEN BEGIN
		message, /info, 'illegal ndec; must be a factor of '+	$
			strjoin(strcompress([ndec_in,mdec_in],/rem),' and ')
		RETURN, -1
	ENDIF

	dra_out  = 360.0d0/nra_out 		; Output resolutions (degrees per bin)
	ddec_out = 180.0d0/ndec_out

	cra_out  = -0.5d0				; Index for RA  = 0 in output equatorial map
	cdec_out =  0.5d0*(ndec_out-1)	; Index for Dec = 0 in output equatorial map

	IF silent LE 0 THEN $
		message, /info, 'rebin ' +strjoin(strcompress([nra_in  ,ndec_in  ],/rem),'x')+ $
			' equatorial map by '+strjoin(strcompress([nra_bins,ndec_bins],/rem),'x')+ $
			' to '				 +strjoin(strcompress([nra_out ,ndec_out ],/rem),'x')

END

ENDCASE

destroyvar, map, num

empty_sky = {smei_sky_hdr}

destroyvar, hdrs

FOR icam=0,count-1 DO BEGIN			; Merge cameras in order specified

	skyfile = given_map[icam]
	skytime = smei_filename(skyfile,camera=camera,mode=mode)

	skyhdr = empty_sky

	CASE 1 OF

	hires_map: BEGIN							; Get high-resolution map

		IF silent LE 2 THEN message, /info, hide_env(skyfile)

		CASE (file_search(skyfile))[0] EQ '' OF

		0: BEGIN

			sky = readfits(skyfile,hdr,exten_no=exten_no,silent=silent+1)

			skyhdr.file  = skyfile
			skyhdr.camera= camera
			skyhdr.mode  = mode

			main_hdr = headfits(skyfile, silent=silent)

			skyhdr.stime = TimeSet(fxpar(main_hdr,'STIME'))
			skyhdr.etime = TimeSet(fxpar(main_hdr,'ETIME'))

			skyhdr.oneorbit = fxpar(main_hdr,'ONEORBIT')
			skyhdr.orbit	= fxpar(main_hdr,'ORBIT'   )

			skyhdr.map		= fxpar(hdr,'MAP')

			skyhdr.cx		= fxpar(hdr, 'CRPIX1')
			IF !err GE 0 THEN BEGIN
				skyhdr.cx  -= 1
				skyhdr.cy	= fxpar(hdr, 'CRPIX2')-1
				skyhdr.bd	= 1.0/fxpar(hdr, 'CDELT2')
			ENDIF ELSE BEGIN
				skyhdr.cx	= fxpar(hdr, 'CXEQ')
				IF !err GE 0 THEN BEGIN
					skyhdr.cx -= 1						; CXEQ is for a 1-based array
					skyhdr.cy  = fxpar(hdr, 'CYEQ')-1	; CYEQ is for a 1-based array
					skyhdr.bd  = fxpar(hdr, 'D_EQ')		; Bins per degree
				ENDIF ELSE BEGIN
					skyhdr.cx  = fxpar(hdr, 'CXPL')-1	; CXPL is for a 1-based array
					skyhdr.cy  = fxpar(hdr, 'CYPL')-1	; CYPL is for a 1-based array
					skyhdr.bd  = fxpar(hdr, 'D_PL')		; Bins per degree
				ENDELSE
			ENDELSE

			CASE IsType(sky,/array) OF

			0: sky = float(sky) 					; Scalar; -1=read error

			1: BEGIN

				IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

				mtime = TimeOp(skyhdr.stime,skyhdr.etime,/mean)
				mgain = smei_camera_gain(mtime,camera=skyhdr.camera)

				CASE 0 OF
				adus_map: gain = 1.0
				fixgain : gain = 1.0
				ELSE: BEGIN
					gain = mgain
					IF silent LE 0 THEN message, /info, $
						'camera'+strcompress(skyhdr.camera) + $
						' gain correction of'+strcompress(gain)
				END
				ENDCASE

				IF adus_map THEN BEGIN

					IF rmbkgnd THEN BEGIN
						message, /info, 'subtract sidereal background'
						CASE smei_sky_hdr2mode(hdr) OF
						'eqmap': bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=0, /silent)
						'npmap': bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=1, /silent)
						'spmap': bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=2, /silent)
						ENDCASE
						n = where(sky NE 0)
						IF n[0] NE -1 THEN sky[n] -= bkgnd[n]*mgain
					ENDIF

					IF rmzld THEN BEGIN
						message, /info, 'subtract zodial light distribution'
						n = where(sky NE 0)
						IF n[0] NE -1 THEN BEGIN
							zld = ArrayLocation(n,dim=size(sky,/dim))
							zld = cvsmei(from_map=zld,mode=smei_sky_hdr2mode(hdr),/to_equatorial)
							sky[n] -= smei_zld_model(mtime,zld,/equatorial)*mgain
						ENDIF
					ENDIF

				ENDIF

				IF IsType(map,/undefined) THEN BEGIN
					n = size(sky,/dim)
					map = make_array(type=IsType(sky   ), dim=n)
					num = make_array(type=IsType(camera), dim=n)

					cx  = skyhdr.cx
					cy  = skyhdr.cy
					bd  = skyhdr.bd
				ENDIF

				n = where(sky NE 0)
				IF n[0] NE -1 THEN BEGIN
					map[n] = sky[n]/gain
					num[n] = camera
				ENDIF

				bad_data = fxpar(main_hdr,'BAD_DATA')

			END

			ENDCASE

		END

		1: BEGIN

			IF silent LE 1 THEN message, /info, 'no '+hide_env(skyfile)
			destroyvar, hdr
			sky = -1.0

		END

		ENDCASE

	END

	lores_map: BEGIN		; Extract one of the low-resolution maps

		IF silent LE 2 THEN message, /info, hide_env(skyfile)

		CASE (file_search(skyfile))[0] EQ '' OF
		0: sky = readfits(skyfile, hdr, silent=silent+1, exten_no=exten_no)
		1: BEGIN
			IF silent LE 1 THEN message, /info, 'no '+hide_env(skyfile)
			sky = -1.0
		END
		ENDCASE

		CASE IsType(sky,/array) OF

		0: sky = float(sky)

		1: BEGIN

			IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

			IF angle_map THEN sky /= dpm

			skyhdr.file  = skyfile
			skyhdr.camera= camera
			skyhdr.mode  = mode

			main_hdr = headfits(skyfile, silent=silent)

			skyhdr.stime = TimeSet(fxpar(main_hdr,'STIME'))
			skyhdr.etime = TimeSet(fxpar(main_hdr,'ETIME'))

			skyhdr.oneorbit = fxpar(main_hdr,'ONEORBIT')
			skyhdr.orbit	= fxpar(main_hdr,'ORBIT')
			skyhdr.map		= fxpar(hdr, 'MAP')

			skyhdr.cx		= fxpar(hdr, 'CRPIX1')
			IF !err GE 0 THEN BEGIN
				skyhdr.cx  -= 1
				skyhdr.cy	= fxpar(hdr, 'CRPIX2')-1
				skyhdr.bd	= 1.0/fxpar(hdr, 'CDELT2')
			ENDIF ELSE BEGIN
				skyhdr.cx	= fxpar(hdr, 'CXLO')-1	; CXLO is for a 1-based array
				skyhdr.cy	= fxpar(hdr, 'CYLO')-1	; CYLO is for a 1-based array
				skyhdr.bd	= fxpar(hdr, 'D_LO')
			ENDELSE

			mtime = TimeOp(skyhdr.stime,skyhdr.etime,/mean)
			mgain = smei_camera_gain(mtime,camera=skyhdr.camera)

			CASE 0 OF
			adus_map: gain = 1.0
			fixgain : gain = 1.0
			ELSE: BEGIN
				gain = mgain
				IF silent LE 0 THEN message, /info, $
					'camera'+strcompress(skyhdr.camera) + $
					' gain correction of'+strcompress(gain)
			END
			ENDCASE

			IF IsType(map,/undefined) THEN BEGIN

				cra  = skyhdr.cx
				cdec = skyhdr.cy
				bd	 = skyhdr.bd

				ndim = fxpar(hdr, 'NAXIS*')

				xgrid = (gridgen(ndim[0])-cra )/bd/dpm
				ygrid = (gridgen(ndim[1])-cdec)/bd/dpm

				xedge = (([-0.5d0,ndim[0]-0.5d0])-cra )/bd/dpm
				yedge = (([-0.5d0,ndim[1]-0.5d0])-cdec)/bd/dpm

				rgrid = fltarr(2,ndim[0]*ndim[1])
				rgrid[0,*] = xgrid#replicate(1,ndim[1])
				rgrid[1,*] = replicate(1,ndim[0])#ygrid

			ENDIF

			IF adus_map THEN BEGIN

				IF rmbkgnd THEN BEGIN
					message, /info, 'sidereal bkgnd subtraction skipped'
				ENDIF

				IF rmzld THEN BEGIN
					message, /info, 'subtract zodiacal light distribution'
					n = where(sky NE 0)
					IF n[0] NE -1 THEN $
						sky[n] -= smei_zld_model(mtime,rgrid[*,n],/equatorial,degrees=degrees)*mgain
				ENDIF

			ENDIF

			IF IsType(map,/undefined) THEN BEGIN
				n = size(sky,/dim)
				map = make_array(type=IsType(sky   ), dim=n)
				num = make_array(type=IsType(camera), dim=n)
			ENDIF

			n = where(sky NE 0)
			IF n[0] NE -1 THEN BEGIN
				map[n] = sky[n]/gain
				num[n] = camera
			ENDIF

			bad_data = fxpar(main_hdr,'BAD_DATA')

		END

		ENDCASE

	END

	eclsun_map: BEGIN		; Extract one of the helioecliptic maps

		IF silent LE 2 THEN message, /info, hide_env(skyfile)

		CASE (file_search(skyfile))[0] EQ '' OF
		0: sky = readfits(skyfile, hdr, silent=silent+1, exten_no=exten_no)
		1: BEGIN
			IF silent LE 1 THEN message, /info, 'no '+hide_env(skyfile)
			sky = -1.0
		END
		ENDCASE

		CASE IsType(sky,/array) OF

		0: sky = float(sky)

		1: BEGIN

			IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

			IF angle_map THEN sky /= dpm

			skyhdr.file  = skyfile
			skyhdr.camera= camera
			skyhdr.mode  = mode

			main_hdr = headfits(skyfile, silent=silent)

			skyhdr.stime = TimeSet(fxpar(main_hdr,'STIME'))
			skyhdr.etime = TimeSet(fxpar(main_hdr,'ETIME'))

			skyhdr.oneorbit = fxpar(main_hdr,'ONEORBIT')
			skyhdr.orbit	= fxpar(main_hdr,'ORBIT')
			skyhdr.map		= fxpar(hdr, 'MAP')

			skyhdr.cx		= fxpar(hdr, 'CRPIX1')
			IF !err GE 0 THEN BEGIN
				skyhdr.cx  -= 1
				skyhdr.cy	= fxpar(hdr, 'CRPIX2')-1
				skyhdr.bd	= 1.0/fxpar(hdr, 'CDELT2')
			ENDIF ELSE BEGIN
				skyhdr.cx	= fxpar(hdr, 'CXLO')-1	; CXLO is for a 1-based array
				skyhdr.cy	= fxpar(hdr, 'CYLO')-1	; CYLO is for a 1-based array
				skyhdr.bd	= fxpar(hdr, 'D_LO')
			ENDELSE

			mtime = TimeOp(skyhdr.stime,skyhdr.etime,/mean)
			mgain = smei_camera_gain(mtime,camera=skyhdr.camera)

			CASE 0 OF
			adus_map: gain = 1.0
			fixgain : gain = 1.0
			ELSE: BEGIN
				gain = mgain
				IF silent LE 0 THEN message, /info, $
					'camera'+strcompress(skyhdr.camera) + $
					' gain correction of'+strcompress(gain)
			END
			ENDCASE

			IF IsType(map,/undefined) THEN BEGIN

				cra  = skyhdr.cx
				cdec = skyhdr.cy
				bd	 = skyhdr.bd

				ndim = fxpar(hdr, 'NAXIS*')

				xgrid = (gridgen(ndim[0])-cra )/bd/dpm
				ygrid = (gridgen(ndim[1])-cdec)/bd/dpm

				xedge = (([-0.5d0,ndim[0]-0.5d0])-cra )/bd/dpm
				yedge = (([-0.5d0,ndim[1]-0.5d0])-cdec)/bd/dpm

				rgrid = fltarr(2,ndim[0]*ndim[1])
				rgrid[0,*] = xgrid#replicate(1,ndim[1])
				rgrid[1,*] = replicate(1,ndim[0])#ygrid

			ENDIF

			IF adus_map THEN BEGIN

				IF rmbkgnd THEN BEGIN
					message, /info, 'sidereal bkgnd subtraction skipped'
				ENDIF

				IF rmzld THEN BEGIN
					message, /info, 'zodiacal light subtraction skipped'
				ENDIF

			ENDIF

			IF IsType(map,/undefined) THEN BEGIN
				n = size(sky,/dim)
				map = make_array(type=IsType(sky   ), dim=n)
				num = make_array(type=IsType(camera), dim=n)
			ENDIF

			n = where(sky NE 0)
			IF n[0] NE -1 THEN BEGIN
				map[n] = sky[n]/gain
				num[n] = camera
			ENDIF

			bad_data = fxpar(main_hdr,'BAD_DATA')

		END

		ENDCASE

	END

	frame_map: BEGIN

		IF silent LE 2 THEN message, /info, hide_env(skyfile)

		CASE (file_search(skyfile))[0] EQ '' OF
		0: sky = readfits(skyfile, hdr, exten_no=exten_no, silent=silent+1)
		1: BEGIN
			IF silent LE 1 THEN message, /info, 'no '+hide_env(skyfile)
			sky = -1.0
		END
		ENDCASE

		CASE IsType(sky,/array) OF

		0: sky = float(sky)

		1: BEGIN

			IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

			skyhdr.file  = skyfile
			skyhdr.camera= camera
			skyhdr.mode  = mode

			main_hdr = headfits(skyfile, silent=silent)

			skyhdr.stime = TimeSet(fxpar(main_hdr,'STIME'))
			skyhdr.etime = TimeSet(fxpar(main_hdr,'ETIME'))

			skyhdr.oneorbit = fxpar(main_hdr,'ONEORBIT')
			skyhdr.orbit	= fxpar(main_hdr,'ORBIT')
			skyhdr.map		= fxpar(hdr,'MAP')

			skyhdr.cx	 = -1.0
			skyhdr.cy	 = -1.0
			skyhdr.bd	 = -1.0

			IF IsType(map,/undefined) THEN BEGIN
				n = size(sky,/dim)
				map = make_array(type=IsType(sky   ), dim=n)
				num = make_array(type=IsType(camera), dim=n)

				ndim = fxpar(hdr, 'NAXIS*')
			ENDIF

			n = where(sky NE 0)
			IF n[0] NE -1 THEN BEGIN
				map[n] = sky[n]
				num[n] = camera
			ENDIF

			bad_data = fxpar(main_hdr,'BAD_DATA')

		END

		ENDCASE

	END

	ELSE: BEGIN 			; Combine high-resolution maps into single map

		IF silent LE 2 THEN message, /info, hide_env(skyfile)

		CASE (file_search(skyfile))[0] EQ '' OF
		0: BEGIN
			ext = use_flat*3+0		; Read equatorial map
			sky = readfits(skyfile, hdr, ext=ext, silent=silent+1)
		END
		1: BEGIN
			IF silent LE 1 THEN message, /info, 'no '+hide_env(skyfile)
			sky = -1.0
		END
		ENDCASE

		CASE IsType(sky,/array) OF

		0: BEGIN
			sky = float(sky)
			;message, /info, hide_env(skyfile)+' error reading equatorial map (Fits ext 0)'
		END

		1: BEGIN

			IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

			skyhdr.file  = skyfile
			skyhdr.camera= camera
			skyhdr.mode  = mode

			skyhdr.stime = TimeSet(fxpar(hdr,'STIME'))
			skyhdr.etime = TimeSet(fxpar(hdr,'ETIME'))

			skyhdr.oneorbit = fxpar(hdr,'ONEORBIT')
			skyhdr.orbit	= fxpar(hdr,'ORBIT')

			skyhdr.map		= 'Equatorial map RA=[0,360],DEC=[-90,+90]'

			skyhdr.cx = cra_out
			skyhdr.cy = cdec_out
			skyhdr.bd = 1.0d0/dra_out

			mtime = TimeOp(skyhdr.stime,skyhdr.etime,/mean)
			mgain = smei_camera_gain(mtime,camera=skyhdr.camera)

			bad_data = fxpar(hdr,'BAD_DATA')

			IF rmbkgnd THEN BEGIN
				message, /info, 'subtract sidereal background'
				bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=0, /silent)
				n = where(sky NE 0)
				IF n[0] NE -1 THEN sky[n] -= bkgnd[n]*mgain
			ENDIF

			cam_map = fltarr(nra_out,ndec_out)

			ndim	= fxpar(hdr, 'NAXIS*')

			cra_in = fxpar(hdr, 'CRPIX1')
			IF !err GE 0 THEN BEGIN
				cra_in -= 1
				cdec_in = fxpar(hdr, 'CRPIX2')-1
				dd		= 1.0/fxpar(hdr, 'CDELT2')
			ENDIF ELSE BEGIN
				cra_in	= fxpar(hdr, 'CXEQ'  )-1	; CXEQ is for a 1-based array
				cdec_in = fxpar(hdr, 'CYEQ'  )-1	; CYEQ is for a 1-based array
				dd		= fxpar(hdr, 'D_EQ'  )
			ENDELSE

			; Sanity check (probably redundant)

			;IF ndim[0] NE nra_in OR ndim[1] NE mdec_in OR	$
			;	dd*dra_in NE 1.0d0 OR dd*ddec_in NE 1.0d0 THEN $
			;		message, 'unexpected dimensions in equatorial map'
			IF dd*dra_in NE 1.0d0 OR dd*ddec_in NE 1.0d0 THEN $
				message, 'unexpected resolution in equatorial map'

			cnt = float(sky NE 0.0)

			; Rebin to output resolution by averaging over nra_bins x ndec_bins bins.

			sdim = ndim/[nra_bins,ndec_bins]

			IF nbin*sdim[0] NE ndim[0] THEN BEGIN
				sky = sky[0:sdim[0]*nbin-1,*]
				cnt = cnt[0:sdim[0]*nbin-1,*]
			ENDIF
			IF nbin*sdim[1] NE ndim[1] THEN BEGIN
				sky = sky[*,0:sdim[1]*nbin-1]
				cnt = cnt[*,0:sdim[1]*nbin-1]
			ENDIF
			ndim = nbin*sdim

			sky = reform(sky,nra_bins,sdim[0],ndec_bins,sdim[1],/overwrite)
			cnt = reform(cnt,nra_bins,sdim[0],ndec_bins,sdim[1],/overwrite)

			sky = total(total(sky,3),1)
			cnt = total(total(cnt,3),1)

			n = where(cnt NE 0)
			IF n[0] NE -1 THEN sky[n] /= cnt[n]

			; To find out where to insert array sky into array map we need to
			; know where bin (0,0) of the input grid is put in the output grid

			ra0  = round( cra_out +((0-cra_in )*dra_in )/dra_out  )
			dec0 = round( cdec_out+((0-cdec_in)*ddec_in)/ddec_out )

			cam_map[ra0:ra0+sdim[0]-1,dec0:dec0+sdim[1]-1] = sky
			;cam_map[*,dec0:dec0+mdec_out-1] = sky

			IF dec0 NE 0 THEN BEGIN									; Map not filled to the south

				ext = use_flat*3+1
				sky = readfits(skyfile, hdr, ext=ext, silent=silent+1); Read map of north pole

				CASE IsType(sky,/array) OF

				0: BEGIN
					sky = float(sky)
					message, /info, 'no N-pole map (ext'+ $
						strcompress(ext)+'), '+hide_env(skyfile)
				END

				1: BEGIN

					IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

					ndim = fxpar(hdr, 'NAXIS*')

					cx = fxpar(hdr, 'CRPIX1')
					IF !err GE 0 THEN BEGIN
						cx -= 1
						cy  = fxpar(hdr, 'CRPIX2')-1
						dd  = 1.0/fxpar(hdr, 'CDELT2')
					ENDIF ELSE BEGIN
						cx  = fxpar(hdr, 'CXPL')-1
						cy  = fxpar(hdr, 'CYPL')-1
						dd  = fxpar(hdr, 'D_PL')
					ENDELSE

					IF rmbkgnd THEN BEGIN
						bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=1, /silent)
						n = where(sky NE 0)
						IF n[0] NE -1 THEN sky[n] -= bkgnd[n]*mgain
					ENDIF

					ra  = (gridgen(nra_out)-cra_out)*dra_out

					; Latitude at top edge of equatorial map (at index mdec_in-0.5d0)
					; and number of latitudes to be filled at north pole.

					idec = (mdec_in-0.5d0-cdec_in)*ddec_in

					idec = cdec_out+idec/ddec_out
					idec = ceil(idec)
					mdec = ndec_out-idec

					dec = (idec+gridgen(mdec)-cdec_out)*ddec_out

					ra  = ra#replicate(1,mdec)
					dec = replicate(1,nra_out)#dec

					ra  = ra*(!dpi/180.0d0)
					dec = dd*(90.0d0-dec)

					xx  = round(cx+dec*cos(ra))
					yy  = round(cy+dec*sin(ra))

					dx  = round(dra_out*cos(!dpi/3)*dd/2.0d0)
					dy  = dx

					pac = fltarr(nra_out,mdec,2*dx+1,2*dy+1)
					FOR i=0,2*dx DO FOR j=0,2*dy DO pac[*,*,i,j] = sky[i-dx+xx,j-dy+yy]
					pac = reform(pac,nra_out,mdec,(2*dx+1)*(2*dy+1),/overwrite)

					cnt = float(pac NE 0.0) 
					cnt = reform(cnt,nra_out,mdec,(2*dx+1)*(2*dy+1),/overwrite)

					pac = total(pac,3)
					cnt = total(cnt,3)

					n = where(cnt NE 0)
					IF n[0] NE -1 THEN pac[n] /= cnt[n]

					cam_map[*,idec:*] = pac

				END

				ENDCASE

			ENDIF

			IF dec0+mdec_out NE ndec_out THEN BEGIN 				; Map not filled in the north

				ext = use_flat*3+2
				sky = readfits(skyfile, hdr, ext=ext, silent=silent+1); Read map of south pole

				CASE IsType(sky,/array) OF

				0: BEGIN											; Read map of south pole
					sky = float(sky)
					message, /info, 'no S-pole map (ext'+ $
						strcompress(ext)+'), '+hide_env(skyfile)
				END

				1: BEGIN

					IF silent LE 3 THEN message, /info, fxpar(hdr,'MAP')

					ndim = fxpar(hdr, 'NAXIS*')

					cx = fxpar(hdr, 'CRPIX1')
					IF !err GE 0 THEN BEGIN
						cx -= 1
						cy	= fxpar(hdr, 'CRPIX2')-1
						dd	= 1.0/fxpar(hdr, 'CDELT2')
					ENDIF ELSE BEGIN
						cx	= fxpar(hdr, 'CXPL')-1
						cy	= fxpar(hdr, 'CYPL')-1
						dd	= fxpar(hdr, 'D_PL')
					ENDELSE

					IF rmbkgnd THEN BEGIN
						bkgnd = readfits(bkgnd_map, bkgnd_hdr, exten_no=2, /silent)
						n = where(sky NE 0)
						IF n[0] NE -1 THEN sky[n] -= bkgnd[n]*mgain
					ENDIF

					ra  = (gridgen(nra_out)-cra_out)*dra_out

					; Latitude at bottom edge of equatorial map (at index -0.5d0)

					idec = (-0.5d0-cdec_in)*ddec_in

					idec = cdec_out+idec/ddec_out
					idec = floor(idec)
					mdec = idec+1

					dec = (gridgen(mdec)-cdec_out)*ddec_out

					ra  = ra#replicate(1,mdec)
					dec = replicate(1,nra_out)#dec

					ra  = ra*(!dpi/180.0d0)
					dec = dd*(dec+90.0d0)

					xx  = round(cx+dec*cos(ra))
					yy  = round(cy+dec*sin(ra))

					dx  = round(dra_out*cos(!dpi/3)*dd/2.0d0)
					dy  = dx

					pac = fltarr(nra_out,mdec,2*dx+1,2*dy+1)
					FOR i=0,2*dx DO FOR j=0,2*dy DO pac[*,*,i,j] = sky[i-dx+xx,j-dy+yy]
					pac = reform(pac,nra_out,mdec,(2*dx+1)*(2*dy+1),/overwrite)

					cnt = float(pac NE 0.0) 
					cnt = reform(cnt,nra_out,mdec,(2*dx+1)*(2*dy+1),/overwrite)

					pac = total(pac,3)
					cnt = total(cnt,3)

					n = where(cnt NE 0)
					IF n[0] NE -1 THEN pac[n] /= cnt[n]

					cam_map[*,0:idec] = pac

				END

				ENDCASE

			ENDIF


			CASE 0 OF
			fixgain: gain = 1.0
			ELSE: BEGIN
				gain = mgain
				IF silent LE 0 THEN message, /info, $
					'camera'+strcompress(skyhdr.camera) + $
					' gain correction of'+strcompress(gain)
			END
			ENDCASE

			IF IsType(map, /undefined) THEN BEGIN

				xgrid = (gridgen(nra_out )-cra_out )*dra_out /dpm
				ygrid = (gridgen(ndec_out)-cdec_out)*ddec_out/dpm

				xedge = (([-0.5d0,nra_out -0.5d0])-cra_out )*dra_out /dpm
				yedge = (([-0.5d0,ndec_out-0.5d0])-cdec_out)*ddec_out/dpm

				rgrid = fltarr(2,nra_out*ndec_out)
				rgrid[0,*] = xgrid#replicate(1,ndec_out)
				rgrid[1,*] = replicate(1,nra_out)#ygrid

			ENDIF

			IF rmzld THEN BEGIN
				message, /info, 'subtract zodiacal light distribution'
				n = where(cam_map NE 0.0)
				IF n[0] NE -1 THEN $
					cam_map[n] -= smei_zld_model(mtime,rgrid[*,n],/equatorial,degrees=degrees)*mgain
			ENDIF

			CASE IsType(map, /defined) OF
			0: BEGIN
				n = size(cam_map,/dim)
				map = make_array(type=IsType(cam_map), dim=n)
				num = make_array(type=IsType(camera ), dim=n)
			END
			1: BEGIN
				n = where(map NE 0.0 AND cam_map NE 0.0, cnt)
				IF n[0] NE -1 THEN		$
					message, /info, 'cam'+strcompress(camera) + $
						' overlaps in'+strcompress(cnt) 			+ $
						' pix with mean diff of'+strcompress(mean(cam_map[n]-map[n]))
			END
			ENDCASE

			n = where(cam_map NE 0.0)
			IF n[0] NE -1 THEN BEGIN
				map[n] = cam_map[n]/gain
				num[n] = camera
			ENDIF

		END

		ENDCASE

	END

	ENDCASE

	boost, hdrs, skyhdr

ENDFOR

CASE 1 OF

IsType(map,/undefined): map = -1

hires_map: BEGIN

	ndim = size(map,/dim)

	IF nbin NE 1 THEN BEGIN

		cnt = float(map NE 0)

		; Rebin to output dimension by averaging over nbin x nbin bins.

		sdim = ndim/nbin
		IF nbin*sdim[0] NE ndim[0] THEN BEGIN
			map = map[0:sdim[0]*nbin-1,*]
			num = num[0:sdim[0]*nbin-1,*]
			cnt = cnt[0:sdim[0]*nbin-1,*]
		ENDIF
		IF nbin*sdim[1] NE ndim[1] THEN BEGIN
			map = map[*,0:sdim[1]*nbin-1]
			num = num[*,0:sdim[1]*nbin-1]
			cnt = cnt[*,0:sdim[1]*nbin-1]
		ENDIF
		ndim = nbin*sdim

		map = reform(map,nbin,sdim[0],nbin,sdim[1],/overwrite)
		num = reform(num,nbin,sdim[0],nbin,sdim[1],/overwrite)
		cnt = reform(cnt,nbin,sdim[0],nbin,sdim[1],/overwrite)

		map = total(total(map,3),1)
		num = total(total(num,3),1)
		cnt = total(total(cnt,3),1)

		n = where(cnt NE 0)
		IF n[0] NE -1 THEN BEGIN
			map[n] /= cnt[n]
			num[n] /= cnt[n]
		ENDIF

		n = where(cnt EQ 0)
		IF n[0] NE -1 THEN BEGIN
			map[n] = BadValue(map)
			num[n] = 0
		ENDIF

	ENDIF

	; Indices of output array in terms of index in input array

	xgrid = 0.5d0*(nbin-1)+nbin*gridgen(ndim[0]/nbin)
	ygrid = 0.5d0*(nbin-1)+nbin*gridgen(ndim[1]/nbin)

	xgrid = (xgrid-cx)/bd/dpm
	ygrid = (ygrid-cy)/bd/dpm

	xedge = 0.5d0*(nbin-1)+nbin*[-0.5d0,(ndim[0]-1)/nbin+0.5d0]
	yedge = 0.5d0*(nbin-1)+nbin*[-0.5d0,(ndim[1]-1)/nbin+0.5d0]

	xedge = (xedge-cx)/bd/dpm
	yedge = (yedge-cy)/bd/dpm

END

lores_map:

frame_map:

ELSE:

ENDCASE

IF size(map,/n_dim) NE 0 THEN BEGIN
	n = where(map EQ bad_data)
	IF n[0] NE -1 THEN BEGIN
		map[n] = BadValue(map)
		num[n] = 0
	ENDIF
ENDIF

RETURN, map  &  END