C+
C NAME:
C	smei_skyd_pixel
C PURPOSE:
C	Construct list of nodes (skybins) to which a given pixel contributes 
C CATEGORY:
C	smei/camera/for/sky
C CALLING SEQUENCE:
	integer function smei_skyd_pixel(level,ra_pix,dec_pix,max_count,node_list)
C INPUTS:
C	level		integer			determines resolution of pool of nodes
C						(1 <= level <= SKYD__LEVEL=3
C	ra_pix (0:4)	double precision	RA  for center and four corners of a pixel
C	dec_pix(0:4)	double precision	dec for center and four corners of a pixel
C
C	max_count	integer			Max allowed number of nodes
C OUTPUTS:
C	smei_skyd_pixel	integer			number of nodes on node_list
C	node_list(*)	integer			node id
C CALLS:
C	Say
C INCLUDE:
	include		'smei_skyd_dim.h'
C PROCEDURE:
C	The name of a node is related to the linear array index into one
C	of the skymaps (equatorial, south and north pole).
C	The equatorial map has NEQ = (level*SKYD__NXEQ)*(level*SKYD__NYEQ) bins.
C	The two polar maps have NPL = (level*SKYD__NXPL)*(level*SKYD__NYPL) bins.
C
C	Each node corresponds to a skybin in these maps.
C	An equatorial node has node id 1 ... NEQ
C	An north polar node has node id NEQ+1 ... NEQ+NPL
C	A south polar node has node id NEQ+NPL+1 ... NEQ+2*NPL
C
C	Program execution is aborted if there are more than max_count nodes
C MODIFICATION HISTORY:
C	FEB-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		level
	    double precision	ra_pix (0:*)
	    double precision	dec_pix(0:*)
	    integer		max_count
	    integer		node_list(*)

	double precision	xeq_center
	double precision	yeq_center
	double precision	xeq_corner(4)
	double precision	yeq_corner(4)

	double precision	xpl_center
	double precision	ypl_center
	double precision	xpl_corner(4)
	double precision	ypl_corner(4)

	double precision	binsize
	double precision	xorig
	double precision	yorig

	double precision	dval
	double precision	rx
	double precision	ry

	double precision	dsind			! Needed for GNU compiler
	double precision	dcosd

	character		cInt2Str*14
	logical			inside_area
	integer			count

	include			'smei_skyd_fnc.h'

	! RA =   0+delta -> i = nint(   0.5+10*delta) = 1
	! RA = 360-delta -> i = nint(3600.5-10*delta) = 3600
	! RA =   0 -> i = nint(   0.5) = 1
	! RA = 360 -> i = nint(3600.5) = 3601
	! So as long as 0 <= ra < 360 -> 1 <= i <= 3600

	! Dec -60+delta -> j = nint(   0.5+10*delta) = 1
	! Dec +60-delta -> j = nint(1200.5-10*delta) = 1200

	!j = nint(origin(SKYD__NYEQ,level)+dec_pix(0)*SKYD__D_EQ*level)

	count = 0					! Counter for nodes

	binsize	= level*SKYD__D_EQ
	xorig	= origin(level,0)
	yorig	= origin(level,SKYD__NYEQ)
	ixrange = level*SKYD__NXEQ
	iyrange = level*SKYD__NYEQ

	xeq_center = xorig+ra_pix (0)*binsize
	yeq_center = yorig+dec_pix(0)*binsize

	if (1.0d0-binsize .le. yeq_center .and. yeq_center .le. iyrange+binsize) then

	    do k=1,4

		rx = ra_pix (k)
		ry = dec_pix(k)

		! Corners might lie on the other side of the 0/360 deg meridian.
		! Make sure the corner values are close to the center values.

		dval = rx-ra_pix(0)
		if (abs(dval) .gt. 180.0d0) rx = rx-dsign(360.0d0,dval)

		xeq_corner(k) = xorig+rx*binsize
		yeq_corner(k) = yorig+ry*binsize

		i = nint(xeq_corner(k))
		j = nint(yeq_corner(k))

		if (k .eq. 1) then
		    ix_min = i
		    ix_max = i
		    jy_min = j
		    jy_max = j
		else
		    ix_min = min(ix_min,i)
		    ix_max = max(ix_max,i)
		    jy_min = min(jy_min,j)
		    jy_max = max(jy_max,j)
		end if

	    end do

	    ! Inside equatorial map
	    ! Loop over all sky bins that the pixel overlaps
	    ! If the center of the bin is inside the pixel then put it on node_list

	    do j=max(1,jy_min),min(jy_max,iyrange)

		do i=ix_min,ix_max

		    if (inside_area(dble(i),dble(j),xeq_center,yeq_center,4,xeq_corner,yeq_corner)) then

			count = count+1

			if (count .le. max_count) then
			    ii = i
			    if (i .lt.       1) ii = i+ixrange
			    if (i .gt. ixrange) ii = i-ixrange
			    node_list(count) = ip_eq(level,ii,j)
			end if

		    end if

		end do

	    end do

	end if

	! k=1: North Polar plot (>  33 degrees)
	! k=2: South Polar plot (< -33 degrees)

	! In the polar maps the maximum polar distance is reached in the
	! corners of the array. Find the corner farthest away from the pole
	! and calculate its polar distance.

	binsize	= level*SKYD__D_PL
	xorig	= origin(level,SKYD__NXPL)
	yorig	= origin(level,SKYD__NYPL)
	ixrange = level*SKYD__NXPL
	iyrange = level*SKYD__NYPL

	rx = xorig-0.5d0
	ry = yorig-0.5d0

	dval = (90.0d0-dabs(dec_pix(0)))*binsize	! Polar distance

	if (dval .lt. dsqrt(rx*rx+ry*ry)+binsize) then

	    rx = ra_pix(0)
	    xpl_center = xorig+dval*dcosd(rx)
	    ypl_center = yorig+dval*dsind(rx)

	    do k=1,4

		dval = (90.0d0-dabs(dec_pix(k)))*binsize

		rx = ra_pix(k)
		xpl_corner(k) = xorig+dval*dcosd(rx)
		ypl_corner(k) = yorig+dval*dsind(rx)

		i = nint(xpl_corner(k))
		j = nint(ypl_corner(k))

		if (k .eq. 1) then
		    ix_min = i
		    ix_max = i
		    jy_min = j
		    jy_max = j
		else
		    ix_min = min(ix_min,i)
		    ix_max = max(ix_max,i)
		    jy_min = min(jy_min,j)
		    jy_max = max(jy_max,j)
		end if

	    end do

	    k = (3-sign(1,int(dec_pix(0))))/2		! 1 for North, 2 for South

	    do j=max(1,jy_min),min(jy_max,iyrange)

		do i=max(1,ix_min),min(ix_max,ixrange)

		    if (inside_area(dble(i),dble(j),xpl_center,ypl_center,4,xpl_corner,ypl_corner)) then

			count = count+1
			if (count .le. max_count) node_list(count) = ip_pl(level,i,j,k)

		    end if

		end do

	    end do

	end if

	if (count .eq. 0) then
	    print *, 'RA=',ra_pix(0),' Dec=',dec_pix(0),' level=',level
	    print *, 'ra =[',(ra_pix (i),',',i=1,3),ra_pix (4),']'
	    print *, 'dec=[',(dec_pix(i),',',i=1,3),dec_pix(4),']'
	    print *, 'equator'
	    print *, 'x0=',xeq_center
	    print *, 'y0=',yeq_center
	    print *, 'x=[',(sngl(xeq_corner(i)),',',i=1,3),sngl(xeq_corner(4)),']'
	    print *, 'y=[',(sngl(yeq_corner(i)),',',i=1,3),sngl(yeq_corner(4)),']'
	    print *, 'pole'
	    print *, 'x0 = ',xpl_center
	    print *, 'y0 = ',ypl_center
	    print *, 'x = [',(sngl(xpl_corner(i)),',',i=1,3),sngl(xpl_corner(4)),']'
	    print *, 'y = [',(sngl(ypl_corner(i)),',',i=1,3),sngl(ypl_corner(4)),']'
	    stop 'no votes for pixel?'			! Not supposed to happen
	end if

	if (count .gt. max_count) call Say('skyd_pixel','E','max_count='//cInt2Str(max_count),
     &		'is too small#need at least '//cInt2Str(count))

	smei_skyd_pixel = count

	return
	end