	SUBROUTINE NESTED_PIXNO(EE,IRES,IPIX)
C
	REAL*4 EE(3),E(3)	! input unit vector
	INTEGER IRES		! resolution
	INTEGER IPIX		! pixel number in `nested' scheme
C
C	written by Ned Wright 22-Jul-97 from Hinshaw's verbal description
C	of Gorski's scheme
C
C	NPIX = number of pixels = 12*4**IRES
C	range of IRES 0..13
C	range of IPIX 0..(NPIX-1)
C
C	sphere is cut into an equatorial band and four triangles on the top
C	and 4 on the bottom.  Band is cut into 4 squares and each is cut into
C	four triangles.  total 24 triangles now reassembled into 12 diamonds
C	for equal area need edge of equatorial band at sin(lat) = 2/3 or
C	lat = 41.81^o
C	y = sin(lat) up to this point
C	y_max = 4/3
C	in polar regions, y > 2/3
C	and	d(area) = dy*2*pi*(4/3 -y)*(3/2) = 2*pi*cos(lat)*d(lat)
C	thus 2*(y-2/3)-0.75*(y^2-4/9) = [sin(lat)-2/3]
C	and y = (4/3) - sqrt((4/3)*(1-sin(lat)))
C
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
	PARAMETER (R2V3 = 0.66666666666667)
C	table of powers of 2
	INTEGER*4 IT14(14)/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,
	1	8192/
C	table of powers of 4
	INTEGER*4 IT28(14)/1,4,16,64,256,1024,4096,16384,65536,262144,
	1	1048576,4194304,16777216,67108864/
C	lower left corner X coordinates in units of pi/4
	INTEGER IXLLC(12)/1,3,5,7,0,2,4,6,1,3,5,7/
C	lower left corner Y coordinates in units of 2/3
	INTEGER IYLLC(12)/0,0,0,0,-1,-1,-1,-1,-2,-2,-2,-2/
C	table of integers with bits separated made by SPRCNPX_SET
	INTEGER*4 IX,IY
	COMMON /SPRBIT/ IX(128),IY(128)
	IF (IX(128).EQ.0) CALL SPRBIT_SET
	IF (IRES.LT.0 .OR. IRES.GT.13) STOP 'IRES out of range in PIXEL'
C	need E to be normalized
	EM = EE(1)**2 + EE(2)**2 + EE(3)**2
	IF(ABS(EM-1).GT.5.E-3) THEN
	  F = 1./SQRT(EM)
	ELSE
	  F = 8./(3.+EM*(6.-EM))
	ENDIF
	DO I=1,3
	  E(I) = F*EE(I)
	ENDDO
C
	Z = ABS(E(3))
	X = 0.
	IF (E(2).NE.0. .OR. E(1).NE.0.) X = ATAN2(E(2),E(1))
	IF (X.LT.0.) X=X+TWOPI
	N = X/PIOVER2
	IF (N.GT.3) N=3
	IF (Z.LE.R2V3) THEN
	  YY = E(3)
	  XX = X
	ELSE
	  Y = (4.-SQRT(12.*(E(1)**2+E(2)**2)/(1+Z)))/3.
	  YY = SIGN(Y,E(3))
	  XX = X*(2.-1.5*Y)+(N+0.5)*PIOVER2*(1.5*Y-1.)
	ENDIF
C	face finding logic
	IF (E(3).GT.R2V3) THEN
	  NFACE = N
	ELSE IF (E(3).LT.-R2V3) THEN
	  NFACE = N+8
	ELSE
	  M = X/PIOVER2 + 0.5
	  DX = X/PIOVER2 - M
C	  IF (DX.GT.2.) DX = DX-4.
	  IF (Z.LT.(R2V3-1.33333333*ABS(DX))) THEN
	    NFACE = X/PIOVER2+4.5
	    IF (NFACE.GE.8) NFACE = NFACE-4
	  ELSE
	    NFACE = N
	    IF (E(3).LT.0.) NFACE = NFACE+8
	  ENDIF
	ENDIF
	X = XX - IXLLC(NFACE+1)*PIOVER4
	IF (X.GT.PI) X = X-TWOPI
	Y = YY - IYLLC(NFACE+1)/1.5
C	resolve the diamond into a square
	XP = (X/PIOVER4+1.5*Y)/2.0
	YP = (1.5*Y-X/PIOVER4)/2.0
C	write (*,*) xx,yy,nface
C	write (*,*) x,y
C	write (*,*) xp,yp
	I=IT14(IRES+1)*XP
	J=IT14(IRES+1)*YP
	I=MIN(IT14(IRES+1)-1,I)			!Roundoff protection for cases very
	J=MIN(IT14(IRES+1)-1,J)			!near an edge of the cube
	IL=MOD(I,128)
	IH=I/128
	JL=MOD(J,128)
	JH=J/128
	NINFACE = IX(IL+1)+IY(JL+1)+16384*(IX(IH+1)+IY(JH+1))
	IPIX = NFACE*IT28(IRES+1) + NINFACE
	RETURN
	END

	SUBROUTINE NESTED_CENPIX(IPIX,IRES,E)
	INTEGER IPIX	! nested scheme pixel number, input, 0..(12*4**IRES -1)
	INTEGER IRES	! resolution, range 0..13
C
C	convert "nested" scheme pixel number into a pixel center
C	in unit vector E
C
	REAL*4 E(3)
	INTEGER*2 IX,IY
	COMMON /SPRCNPX/ IX(0:1023),IY(0:1023)
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
C	table of powers of 2
	INTEGER*4 IT14(14)/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,
	1	8192/
C	table of powers of 4
	INTEGER*4 IT28(14)/1,4,16,64,256,1024,4096,16384,65536,262144,
	1	1048576,4194304,16777216,67108864/
C	lower left corner X coordinates in units of pi/4
	INTEGER IXLLC(12)/1,3,5,7,0,2,4,6,1,3,5,7/
C	lower left corner Y coordinates in units of 2/3
	INTEGER IYLLC(12)/0,0,0,0,-1,-1,-1,-1,-2,-2,-2,-2/
C	central longitude index for the triangular pieces
	INTEGER NOFFACE(12)/0,1,2,3,-1,-1,-1,-1,0,1,2,3/
	IF (IY(1023).EQ.0) CALL SPRCNPX_SET
	IF (IRES.LT.0 .OR. IRES.GT.13) STOP 'IRES out of range in CENPIX'
	NPIX = 12*IT28(IRES+1)
	IF (IPIX.LT.0 .OR. IPIX.GE.NPIX) STOP 'IPIX out of range in CENPIX'
	NFACE = IPIX/IT28(IRES+1)
	N=IPIX-NFACE*IT28(IRES+1)
	I=MOD(N,1024)
	N=N/1024
	J=MOD(N,1024)
	K=N/1024
	JX=1024*IX(K)+32*IX(J)+IX(I)
	JY=1024*IY(K)+32*IY(J)+IY(I)
	XP = (JX+0.5)/IT14(IRES+1)
	YP = (JY+0.5)/IT14(IRES+1)
	X = PIOVER4*(XP-YP)
	Y = (XP+YP)/1.5
	YY = Y + IYLLC(NFACE+1)/1.5
	XX = X + IXLLC(NFACE+1)*PIOVER4
	IF (XX.LT.0.) XX = XX+TWOPI
	N = NOFFACE(NFACE+1)
	IF (N.GE.4) N=3
	YA = ABS(YY)
	IF (YA.GT.(2./3.)) THEN
	  SINLAT = 2*YA-0.75*YA*YA-1./3.
	  COSLAT = SQRT((1.+SINLAT)/12.)*(4-3.*YA)
	  SINLAT = SIGN(SINLAT,YY)
	  IF (N.GE.0) THEN
	    DIV = 2.-1.5*YA
	    IF (DIV.GT.0) THEN
	      ALONG = (XX-(N+0.5)*PIOVER2*(1.5*YA-1.))/(2.-1.5*YA)
	    ELSE
	      ALONG = (N+0.5)*PIOVER2
	    ENDIF
	  ELSE
	    ALONG = XX
	  ENDIF
	ELSE
	  SINLAT = YY
	  COSLAT = SQRT(1.-SINLAT**2)
	  ALONG = XX
	ENDIF
	E(3) = SINLAT
	E(1) = COSLAT*COS(ALONG)
	E(2) = COSLAT*SIN(ALONG)
	RETURN
	END

	SUBROUTINE NESTED_EDGES(IPIX,IRES,NP,E)
	INTEGER IPIX	! nested scheme pixel number, input, 0..(12*4**IRES -1)
	INTEGER IRES	! resolution, range 0..13
C
C	convert "nested" scheme pixel number into unit vectors
C	in array E which trace the deges of the pixel. NP points
C	per edge. E must have second dimension GE 4*NP
C
	REAL*4 E(3,*)
	INTEGER*2 IX,IY
	COMMON /SPRCNPX/ IX(0:1023),IY(0:1023)
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
C	lower left corner X coordinates in units of pi/4
	INTEGER IXLLC(12)/1,3,5,7,0,2,4,6,1,3,5,7/
C	lower left corner Y coordinates in units of 2/3
	INTEGER IYLLC(12)/0,0,0,0,-1,-1,-1,-1,-2,-2,-2,-2/
C	central longitude index for the triangular pieces
	INTEGER NOFFACE(12)/0,1,2,3,-1,-1,-1,-1,0,1,2,3/
	INTEGER IDX(4) /0,1,1,0/
	INTEGER IDY(4) /0,0,1,1/
C	table of powers of 2
	INTEGER*4 IT14(14)/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,
	1	8192/
C	table of powers of 4
	INTEGER*4 IT28(14)/1,4,16,64,256,1024,4096,16384,65536,262144,
	1	1048576,4194304,16777216,67108864/
	IF (IY(1023).EQ.0) CALL SPRCNPX_SET
	IF (IRES.LT.0 .OR. IRES.GT.13) STOP 'IRES out of range in EDGES'
	NPIX = 12*IT28(IRES+1)
	IF (IPIX.LT.0 .OR. IPIX.GE.NPIX) STOP 'IPIX out of range in EDGES'
	NFACE = IPIX/IT28(IRES+1)
	N=MOD(IPIX,IT28(IRES+1))
	I=MOD(N,1024)
	N=N/1024
	J=MOD(N,1024)
	K=N/1024
	JX=1024*IX(K)+32*IX(J)+IX(I)
	JY=1024*IY(K)+32*IY(J)+IY(I)
	XO = FLOAT(JX+IDX(1))/IT14(IRES+1)
	YO = FLOAT(JY+IDY(1))/IT14(IRES+1)
	JP = 0
	DO JC=1,4
	  IC = MOD(JC,4)+1
	  XN = FLOAT(JX+IDX(IC))/IT14(IRES+1)
	  YN = FLOAT(JY+IDY(IC))/IT14(IRES+1)
	  DO IP=1,NP
	    F = (IP-1.)/NP
	    XP = (1.-F)*XO + F*XN
	    YP = (1.-F)*YO + F*YN
	    X = PIOVER4*(XP-YP)
	    Y = (XP+YP)/1.5
	    YY = Y + IYLLC(NFACE+1)/1.5
	    XX = X + IXLLC(NFACE+1)*PIOVER4
	    IF (XX.LT.0.) XX = XX+TWOPI
	    N = NOFFACE(NFACE+1)
	    YA = ABS(YY)
	    IF (YA.GT.(2./3.)) THEN
	      SINLAT = 2*YA-0.75*YA*YA-1./3.
	      COSLAT = SQRT((1.+SINLAT)/12.)*(4-3.*YA)
	      SINLAT = SIGN(SINLAT,YY)
	      IF (N.GE.0) THEN
	        DIV = 2.-1.5*YA
	        IF (DIV.GT.0) THEN
	          ALONG = (XX-(N+0.5)*PIOVER2*(1.5*YA-1.))/(2.-1.5*YA)
	        ELSE
	          ALONG = (N+0.5)*PIOVER2
	        ENDIF
	      ELSE
	        ALONG = XX
	      ENDIF
	    ELSE
	      SINLAT = YY
	      COSLAT = SQRT(1.-SINLAT**2)
	      ALONG = XX
	    ENDIF
	    JP = JP+1
	    E(3,JP) = SINLAT
	    E(1,JP) = COSLAT*COS(ALONG)
	    E(2,JP) = COSLAT*SIN(ALONG)
	  ENDDO
	  XO = XN
	  YO = YN
	ENDDO
	RETURN
	END

	SUBROUTINE NESTED_TO_ORDERED(IPIX,IRES,ILAT,ILONG,IHALF,NLONG,JPIX)
	INTEGER IPIX,IRES	! input
	INTEGER ILAT,ILONG,IHALF,JPIX	! output
C
C	convert nested pixel number into lat, long ordered pixel number
C	order of JPIX: 0-3 for northernmost row, 4-11 for next row...
C	N = 2**IRES is the number of pixels for side of each rhombus
C	NPIX = 12*4**IRES = 12*N**2	range of IPIX and JPIX 0..NPIX-1
C	range of IRES 0..13
C	ILAT = 0 for northernmost row to 4*N-2 for southernmost row
C	NLONG is number of longitudes in this row: range 4 to 4*N
C	ILONG = 0..NLONG-1
C	IHALF = 0 or 1
C	longitude = 360*(ILONG+0.5*IHALF)/NLONG degrees
C
	INTEGER*2 KX,KY
	COMMON /SPRCNPX/ KX(0:1023),KY(0:1023)
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
C	lower left corner X coordinates in units of pi/4
	INTEGER IXLLC(12)/1,3,5,7,0,2,4,6,1,3,5,7/
C	lower left corner Y coordinates in units of 2/3
	INTEGER IYLLC(12)/0,0,0,0,-1,-1,-1,-1,-2,-2,-2,-2/
C	table of powers of 2
	INTEGER*4 IT14(14)/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,
	1	8192/
C	table of powers of 4
	INTEGER*4 IT28(14)/1,4,16,64,256,1024,4096,16384,65536,262144,
	1	1048576,4194304,16777216,67108864/
	IF (KY(1023).EQ.0) CALL SPRCNPX_SET
	IF (IRES.LT.0 .OR. IRES.GT.13) 
	1	STOP 'IRES out of range in NESTED_TO_ORDERED'
	NPIX = 12*IT28(IRES+1)
	IF (IPIX.LT.0 .OR. IPIX.GE.NPIX) STOP
	1	'IPIX out of range in IPIX_2_JPIX'
	NFACE = IPIX/IT28(IRES+1)
	N=MOD(IPIX,IT28(IRES+1))
	I=MOD(N,1024)
	N=N/1024
	J=MOD(N,1024)
	K=N/1024
	JX=1024*KX(K)+32*KX(J)+KX(I)
	JY=1024*KY(K)+32*KY(J)+KY(I)
C
	N = IT14(IRES+1)
C	range of JX and JY is 0..N-1
	IY = JX+JY
	IX = JX-JY
	IF (NFACE.LE.3) THEN
	  ILAT = 2*N-2-IY
	  IF (IY.GT.N-1) THEN
	    ILONG = (ILAT+1)*NFACE + (IX+ILAT)/2
	    IHALF = 1
	    NLONG = 4*(ILAT+1)
	  ELSE
	    ILONG = N*NFACE + (N+IX)/2
	    IHALF = MOD(ILAT,2)
	    NLONG = 4*N
	  ENDIF
	ELSE IF (NFACE.LE.7) THEN
	  NLONG = 4*N
	  ILAT = 3*N-2-IY
	  ILONG = N*(NFACE-4) + (2*N+IX)/2 - N
	  IF (ILONG.LE.0) ILONG = ILONG+4*N
	  IF (ILONG.GE.4*N) ILONG = ILONG-4*N
	  IHALF = MOD(ILAT,2)
	ELSE
	  ILAT = 4*N-2-IY
	  IF (IY.LE.N-1) THEN
	    ILONG = (IY+1)*(NFACE-8) + (IX+IY)/2
	    IHALF = 1
	    NLONG = 4*(IY+1)
	  ELSE
	    ILONG = N*(NFACE-8) + (IX+N)/2
	    IHALF = MOD(ILAT,2)
	    NLONG = 4*N
	  ENDIF
	ENDIF
C	compute JPIX
	IF (ILAT.LE.N) THEN
	  JPIX = 2*ILAT*(ILAT+1)+ILONG
	ELSE IF (ILAT.LE.3*N) THEN
	  JPIX = 4*N*(ILAT-N) + 2*N*(N+1) + ILONG
	ELSE
	  JPIX = N*(2+10*N) + 2*(ILAT-3*N)*(N-1 + (N-(ILAT-3*N))) + ILONG
	ENDIF
C
	RETURN
	END

	SUBROUTINE ORDERED_CENPIX(JPIX,IRES,E)
C
C	convert lat-long ordered pixel number JPIX into pixel center
C	unit vector E
C	NPIX = number of pixels = 12*4**IRES
C	range of JPIX: 0..NPIX-1
C	range of IRES: 0..13
C
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
	INTEGER JPIX
	INTEGER IRES
	REAL*4 E(3)
	DATA JRES /0/
	COMMON /JPIX2IPIX/ JRES,N,ILT(0:32768)
	IF (IRES.LT.0 .OR. IRES.GT.13) STOP 
	1	'IRES out of range in ORDERED_CENPIX'
	IF (IRES.NE.JRES) THEN
C	  compute ILT for this res
	  N = 2**IRES
	  DO ILAT = 0,4*N-1
	    IF (ILAT.LE.N) THEN
	      ILT(ILAT) = 2*ILAT*(ILAT+1)
	    ELSE IF (ILAT.LE.3*N) THEN
	      ILT(ILAT) = 4*N*(ILAT-N) + 2*N*(N+1)
	    ELSE
	      ILT(ILAT) = N*(2+10*N) + 2*(ILAT-3*N)*(N-1 + (N-(ILAT-3*N)))
	    ENDIF
	  ENDDO
	  JRES = IRES
	ENDIF
C	normal entry point
	NPIX = 12*N**2
	IF (JPIX.LT.0 .OR. JPIX.GE.NPIX) STOP 
	1	'JPIX out of range in ORDERED_CENPIX'
	ILL = 0
	ILH = 4*N-2
C	binary search for table entry
100	IT = (ILL+ILH)/2
	IF (JPIX.LT.ILT(IT)) THEN
	  ILH = IT-1
	  GO TO 100
	ENDIF
	IF (JPIX.GE.ILT(IT+1)) THEN
	  ILL = IT+1
	  GO TO 100
	ENDIF
	ILAT = IT
	ILONG = JPIX-ILT(ILAT)
	IF (ILAT.LT.N) THEN
	  NLONG = 4*(ILAT+1)
	  IHALF = 1
	ELSE IF (ILAT.LT.3*N) THEN
	  NLONG = 4*N
	  IHALF = MOD(ILAT,2)
	ELSE
	  NLONG = 4*(4*N-(ILAT+1))
	  IHALF = 1
	ENDIF
	Y = (4.-8.*(ILAT+1.0)/(4*N))/3.
	ALONG = TWOPI*(ILONG + 0.5*IHALF)/NLONG
	YA = ABS(Y)
	IF (YA.LE.(2./3.)) THEN
	  E(3) = Y
	  COSLAT = SQRT(1.-E(3)**2)
	ELSE
	  SINLAT = 2*YA-0.75*YA**2-1./3.
	  COSLAT = SQRT((1.+SINLAT)/12.)*(4-3.*YA)
	  E(3) = SIGN(SINLAT,Y)
	ENDIF
	E(1) = COSLAT*COS(ALONG)
	E(2) = COSLAT*SIN(ALONG)
	RETURN
	END

	SUBROUTINE V_TO_XY(E,XX,YY)
C	convert unit vector into x,y coordinates
	REAL*4 E(3)
	REAL XX,YY
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
	CALL NORM(E)
	Z = ABS(E(3))
	X = 0.
	IF (E(2).NE.0. .OR. E(1).NE.0.) X = ATAN2(E(2),E(1))
	IF (X.LT.0.) X=X+TWOPI
	N = X/PIOVER2
	IF (N.GT.3) N=3
	IF (Z.LE.(2./3.)) THEN
	  YY = E(3)
	  XX = X
	ELSE
	  Y = (4.-SQRT(12.*(E(1)**2+E(2)**2)/(1+Z)))/3.
	  YY = SIGN(Y,E(3))
	  XX = X*(2.-1.5*Y)+(N+0.5)*PIOVER2*(1.5*Y-1.)
	ENDIF
	RETURN
	END

	SUBROUTINE XY_TO_V(XX,YY,E)
C	convert XX, YY on equal area plane into unit vector
C	XX in 0..2*pi
C	YY in -4/3 to +4/3
	REAL*4 E(3)
	REAL XX,YY
	PARAMETER (TWOPI = 6.28318530717958, PI = 3.14159265358979)
	PARAMETER (PIOVER2 = 1.57079632679490, PIOVER4 = 0.785398163397448)
C	check inputs - allow extended range
	IF (ABS(YY).GE.3.) STOP 'YY out of range in XY_TO_V'
	IF (XX.LT.-4. .OR. XX.GT.8.)
	1	STOP 'XX out of range in XY_TO_V'
C	fold extended range into correct range
	Y = YY
	X = XX
	IF (Y.GT.1.3333333) THEN
	  Y = 2.666666667-Y
	  X = X+PI
	ENDIF
	IF (Y.LT.-1.3333333) THEN
	  Y = -2.666666667-Y
	  X = X+PI
	ENDIF
	IF (X.LT.0.) X=X+TWOPI
	IF (X.GT.TWOPI) X=X-TWOPI
	N = X/PIOVER2
	IF (N.GE.4) N=3
	YA = ABS(YY)
	IF (YA.GT.(2./3.)) THEN
	  SINLAT = 2*YA-0.75*YA*YA-1./3.
	  COSLAT = SQRT((1.+SINLAT)/12.)*(4-3.*YA)
	  SINLAT = SIGN(SINLAT,Y)
	  DIV = 2.-1.5*YA
	  IF (DIV.GT.0) THEN
	    ALONG = (X-(N+0.5)*PIOVER2*(1.5*YA-1.))/(2.-1.5*YA)
	  ELSE
	    ALONG = (N+0.5)*PIOVER2
	  ENDIF
	ELSE
	  SINLAT = Y
	  COSLAT = SQRT(1.-SINLAT**2)
	  ALONG = X
	ENDIF
	E(3) = SINLAT
	E(2) = COSLAT*SIN(ALONG)
	E(1) = COSLAT*COS(ALONG)
	RETURN
	END

	SUBROUTINE SPRBIT_SET
C
C	sets up bit table
C
	INTEGER*4 IX,IY,K,IP
	SAVE /SPRBIT/
	COMMON /SPRBIT/ IX(128),IY(128)
	DATA IX(128) /0/
	WRITE (*,*) 'SPRBIT_SET called'
	IF(IX(128).EQ.0) THEN			!Set up the lookup tables
	   DO I=1,128				!for converting x,y into
	      J=I-1				!pixel numbers
	      K=0
	      IP=1
10	      IF(J.EQ.0) THEN
		IX(I)=K
		IY(I)=2*K
	      ELSE
		ID=MOD(J,2)
		J=J/2
		K=IP*ID+K
		IP=IP*4
		GO TO 10
	      ENDIF
	   ENDDO
	ENDIF
	RETURN
	END

	SUBROUTINE SPRCNPX_SET
C
	INTEGER*2 JX,JY
	COMMON /SPRCNPX/ JX(0:1023),JY(0:1023)
	DATA JX(1023)/0/
C
	SAVE /SPRCNPX/
	write (*,*) 'SPRCNPX_SET called'
	DO KPIX=0,1023
	  JPIX=KPIX
	  IX=0
	  IY=0
	  IP=1
	  DO WHILE (JPIX.NE.0)			!Break up the pixel number
	    ID=MOD(JPIX,2)			!By even and odd bits to get
	    JPIX=JPIX/2				!IX and IY
	    IX=ID*IP+IX
	    ID=MOD(JPIX,2)
	    JPIX=JPIX/2
	    IY=ID*IP+IY
	    IP=2*IP
	  ENDDO
	  JX(KPIX)=IX
	  JY(KPIX)=IY
	ENDDO
	RETURN
	END
