C+
C NAME:
C	smei_sky_read
C PURPOSE:
C	Read one of the real*4 arrays in an orbital sky map
C CATEGORY:
C	ucsd/camera/for/lib
C CALLING SEQUENCE:
	subroutine smei_sky_read(cFile, mode, nx, ny, skymap, cx, cy, dd, ctype)
C INPUTS:
C	cFile*(*)	character	fully qualified Fits file name
C
C	Sky maps:
C	mode		integer		selects Fits extension
C					= 0: Equatorial map, RA=[0,360],DEC=[-60,+60]
C					= 1: Polar projection of north pole, DEC=[+50,+90]
C					= 2: Polar projection of south pole, DEC=[-50,-90]
C					= 3: (smei_sky_iread) Counts for equatorial map, RA=[0,360],DEC=[-60,+60]
C					= 4: (smei_sky_iread) Counts for map of north pole, DEC=[+50,+90]
C					= 5: (smei_sky_iread) Counts for map of south pole, DEC=[-50,-90]
C					= 6: Dirty sky brightness RA=[0,360],DEC=[-90,+90]
C					= 7: Time (fraction of orbit) RA=[0,360],DEC=[-90,+90]
C					= 8: PSF position angle from north RA=[0,360],DEC=[-90,+90]
C					= 9: PSF position angle from equinox RA=[0,360],DEC=[-90,+90]
C					=10: Discarded response map RA=[0,360],DEC=[-90,+90]
C					=11: Discarded triangles RA=[0,360],DEC=[-90,+90]
C					=12: Theta-x for glare removal
C					=13: Theta-y for glare removal
C					=14: Time (secs since TORIGIN)
C					=15: Direction cosine angle-x
C					=16: Contributing pixel count
C
C	Equ maps:
C	mode		integer		selects Fits extension
C					= 0: Equatorial map, RA=[0,360],DEC=[-60,+60]
C					= 1: Polar projection of north pole, DEC=[+50,+90]
C					= 2: Polar projection of south pole, DEC=[-50,-90]
C					= 3: Equatorial map, RA=[0,360],DEC=[-60,+60] (flat)
C					= 4: Polar projection of north pole, DEC=[+50,+90] (flat)
C					= 5: Polar projection of south pole, DEC=[-50,-90] (flat)
C					= 6: Stars removed from equator, RA=[0,360],DEC=[-60,+60]
C					= 7: Stars removed from north pole, DEC=[+50,+90]
C					= 8: Stars removed from south pole, DEC=[+50,+90]
C					= 9: Time (fraction of orbit) RA=[0,360],DEC=[-90,+90]
C					=10: Time (sec since TORIGIN) RA=[0,360],DEC=[-90,+90]
C					=11: Direction cosine angle-x RA=[0,360],DEC=[-90,+90]
C					=12: Contributing pixel count RA=[0,360],DEC=[-90,+90]
C
C					The equatorial map is 3600x1200; the polar maps 800x800.
C					all others are 720x360 arrays.
C OUTPUTS:
C	nx,ny		integer		dimensions of sky map
C					If there is an open error (usually because the file does not
C					exist) then nx=0 and ny=0 is returned.
C	skymap(nx,ny)	real		sky map
C	cx,cy,dd	double precision
C					constants relating array index to RA and dec.
C					For north polar map:
C					   i=cx+dd*RA*cos(90-decl); j=cy+dd*RA*sin(90-decl)
C					For south polar map:
C					   i=cx+dd*RA*cos(90+decl); j=cy+dd*RA*sin(90+decl)
C					For all other (equatorial) maps:
C					   i=cx+dd*RA; j=cy+dd*decl
C	ctype		character*(*)	string describing data
C INCLUDE:
	include		'smei_frm_layout.h'
	include		'filparts.h'
	include		'ftspar.h'
C CALLS:
C	iGetLun, iFreeLun, say_fts, Say
C	FTNOPN, FTMAHD, FTGISZ, FTGPVE, FTGKYD, FTCLOS, FTGKYS
C SEE ALSO:
C	smei_htm_fts
C PROCEDURE:
C MODIFICATION HISTORY:
C	JAN-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    character		cFile*(*)
	    integer		mode
	    integer		nx
	    integer		ny
	    real		skymap(*)
	    double precision	cx
	    double precision	cy
	    double precision	dd
	    character		ctype*(*)

	character	cSay*13			/'smei_sky_read'/
	character	cSeverity		/'E'/
	character	cStr*(FIL__LENGTH)

	logical		anyf
	integer		nAxes(2)

	iU = iGetLun(cFile)

	istat = 0

	call FTNOPN(iU, cFile, FTS__READONLY, istat)
	if (istat .ne. 0) then
	    iU = iFreeLun(iU)
	    call Say(cSay,'E',cFile,'not found, or open error')
	end if

	! Advance to the requested map

	if (mode .ne. 0) call FTMAHD(iU, mode+1, ihdutype, istat)

	call FTGISZ(iU, 2, nAxes, istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	nx = nAxes(1)
	ny = nAxes(2)

	anyf = .FALSE.
	call FTGPVE(iU, 1, 1, nx*ny, 0.0, skymap, anyf, istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	! Replace zero by bad values

	!call ArrR4Copy(SMEI__FRM_NPIX,pattern,pat_mode_0)
	!call ArrR4Mask(SMEI__FRM_NPIX,pattern,BadR4(),0.0,0.0,0.0,1.0,pat_mode_0)

	! One of the keywords CRPIX1, CXEQ, CXPL or CXLO must be present in the header
	! (at some point only the WCS keywords will be present, i.e. CRPIX1 will be present)
 
	call FTGKYD(iU, 'CRPIX1', cx, cStr, istat)

	if (istat .eq. 0) then

	    call FTGKYD(iU, 'CRPIX2', cy, cStr, istat)
	    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)
	    call FTGKYD(iU, 'CDELT2', dd, cStr, istat)
	    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)
	    dd = 1.0d0/dd

	else

	    istat = 0
	    call FTGKYD(iU, 'CXEQ', cx, cStr, istat)

	    if (istat .eq. 0) then			! Equatorial map

		call FTGKYD(iU, 'CYEQ', cy, cStr, istat)
		if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)
		call FTGKYD(iU, 'D_EQ', dd, cStr, istat)
		if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	    else

		istat = 0
		call FTGKYD(iU, 'CXPL', cx, cStr, istat)

		if (istat .eq. 0) then			! Polar map

		    call FTGKYD(iU, 'CYPL', cy, cStr, istat)
		    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)
		    call FTGKYD(iU, 'D_PL', dd, cStr, istat)
		    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

		else					! Low resolution map

		    istat = 0
		    call FTGKYD(iU, 'CXLO', cx, cStr, istat)
		    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

		    call FTGKYD(iU, 'CYLO', cy, cStr, istat)
		    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)
		    call FTGKYD(iU, 'D_LO', dd, cStr, istat)
		    if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

		end if

	    end if

	end if

	call FTGKYS(iU, 'MAP' , cType, cStr, istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	! Close fts file

	call FTCLOS(iU, istat)
	if (istat .ne. 0) call say_fts(cSay,cSeverity,istat)

	iU = iFreeLun(iU)

	return
	end