function SkyGridSub, x, rData, rClose, fov, toskygrid=ToSkyGrid, eqtransform=EqTransform
;+
; NAME:
;	SkyGridSub
; PURPOSE:
;	Internal use. Should be called by SkyGrid only.
; CATEGORY:
; CALLING SEQUENCE:
;	Result = SkyGridSub(x,rData,rClose, fov, /toskygrid, eqtransform=EqTransform
; INPUTS:
;	x			/toskygrid NOT SET:
;				float array[4,n^3]
;					Initialized in calling program SkyGrid with:
;						x = fltarr(4,n*n*n, /nozero)
;						x[3,*] = 1.		; Needed for transformation with !p.t
;
;				/toskygrid SET:
;				float array[3,*]		Cartesian coordinates in the X,Y,Z-coordinate system
;										(Origin at viewer location; Z-axis along line of sight
;										through center of field of view; +Y axis is 'right';
;										+X axis is 'up')
;	rData
;	rClose
;	fov
; OPTIONAL INPUT PARAMETERS:
;	/toskygrid
;	eqtransform = EqTransform
;					int scalar		1,2 or 3
;						1: preserve distance to line of sight through center of image
;						2: preserve angular distance to line of sight through center of image
;						3: preserve area on sky
; OUTPUTS:
;	Result		/toskygrid NOT SET:
;				float array[3,n^3]		Cartesian coordinates in X,Y,Z-coordinate system
;										where interpolated densities are required for the display
;										density matrix.
;				/toskygrid SET:
;				float array[3,*]		x,y,z IDL graphics coordinates
; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;	FEB-1998, Paul Hick (UCSD/CASS, pphick@ucsd.edu)
;-
on_error, 2

ToSkyGrid	= keyword_set(ToSkyGrid)
fovedge = [ tan(fov), fov, sqrt(2*(1-cos(fov))) ]

if ToSkyGrid then begin
	s = (sqrt(total(x*x,1))-rClose)/rData; Distance to ViewLoc (-rClose, normalize to rData)

	case EqTransform of
	0: begin
		tanfov = tan(fov)				; fov should be in radians
										; Put negative z at positive z outside fov
		x[2  ,*] = x[2,*] > ( ( abs(x[0,*]) < abs(x[1,*]))/(2*tanfov) )
										; Range inside fov: [-1,1]
		x[0:1,*] = ( x[0:1,*] / ([1.,1.]#x[2,*]) )/fovedge(EqTransform)

	end
	1: begin
		xy = sqrt(total(x[0:1,*]*x[0:1,*],1))

		i = where(xy ne 0)
		if i[0] ne -1 then xy[i] = atan(xy[i],x[2,i])/xy[i]
		i = where(xy eq 0)
		if i[0] ne -1 then xy[i] = 1.

		x[0:1,*] = ( x[0:1,*]*([1.,1.]#[xy]) )/fovedge(EqTransform)

	end
	2: begin
		xy = sqrt(total(x[0:1,*]*x[0:1,*],1))

		i = where(xy ne 0)
		if i[0] ne -1 then xy[i] = sqrt(2*(1.-x[2,i]/s[i]))/xy[i]
		i = where(xy eq 0)
		if i[0] ne -1 then xy[i] = 1.

		x[0:1,*] = ( x[0:1,*]*([1.,1.]#[xy]) )/fovedge(EqTransform)

	end
	endcase

	x[0:1,*] = 1.+x[0:1,*]			; Shift x,y range to [0,2]
	x[2  ,*] = 2.-s					; Shift z range to [0,2]; reverse z-axis

	return, x

endif

r = size(x,/dim)
nDimX = r[1]
nDimY = r[2]
nDimZ = r[3]

OneX = replicate(1.,nDimX)
OneY = replicate(1.,nDimY)
OneZ = replicate(1.,nDimZ)

Zero2TwoX = findgen(nDimX)/(nDimX-1)*2.
Zero2TwoY = findgen(nDimY)/(nDimY-1)*2.
Zero2TwoZ = findgen(nDimZ)/(nDimZ-1)*2.

; The field of view is covered by a regular grid of dimension nDimX x nDimY
; Create r[2,nDimX,nDimY] array. r[0,*,*] and r[1,*,*] contain the
; x- and y-coordinates, respectively, across the field of view

r = fltarr(2,nDimX,nDimY)
r[0,*,*] = OneX#(-1.+Zero2TwoY)
r[1,*,*] = (-1.+Zero2TwoX)#OneY

r = fovedge[EqTransform]*r					; 2 x (nDimX x nDimY) grid

r = reform(r,2,nDimX*nDimY, /overwrite)		; 2 x (nDimX*nDimY) array

; The radial distance between rClose and rClose+2*rData is covered by an
; equally spaced grid of nDim elements. Each grid point in the
; field of view combines with each radial element to define a 3-D grid.

; Convert to the Cartesian x"-y"-z" coordinate system (this amounts to
; reversing a perspective transformation).

; The z"-component is calculated first.
; The x"- and y"-component scale proportionally with the z"-component
; The unit of length is the unit used for rData and rClose.

rr = reverse(rClose+Zero2TwoZ*rData)

;===========================================================================
; The Do-loop is faster in a low memory configuration

x = reform(x,4,nDimX*nDimY,nDimZ, /overwrite)

case EqTransform of
0: begin
	ss = 1./sqrt(1+total(r*r,1))

	for i=0,nDimZ-1 do begin
		s = ss*rr[i]
		x[0,*,i] = s*reform(r[0,*])
		x[1,*,i] = s*reform(r[1,*])
		x[2,*,i] = s
	endfor

;===================================
;	x = fltarr(4,nDimX*nDimY*nDimZ, /nozero)
;	s = (1./sqrt(1+total(r*r,1)))#rr
;	x[0,*] = s*(reform(r[0,*])#OneZ)	; All [nDimX*nDimY,nDimZ] arrays
;	x[1,*] = s*(reform(r[1,*])#OneZ)
;	x[2,*] = s

end
1: begin
	s = sqrt(total(r*r,1))

	ss = s
	i = where(s ne 0)
	if i[0] ne -1 then ss[i] = sin(s[i])/s[i]
	i = where(s eq 0)
	if i[0] ne -1 then ss[i] = 1.

	s = cos(s)

	for i=0,nDimZ-1 do begin
		x[0,*,i] = rr[i]*ss*reform(r[0,*])
		x[1,*,i] = rr[i]*ss*reform(r[1,*])
		x[2,*,i] = rr[i]*s
	endfor

end
2: begin
	s = total(r*r,1)
	ss = sqrt(1.-s/4.)
	s = 1.-s/2.

	for i=0,nDimZ-1 do begin
		x[0,*,i] = rr[i]*ss*reform(r[0,*])
		x[1,*,i] = rr[i]*ss*reform(r[1,*])
		x[2,*,i] = rr[i]*s
	endfor

end
endcase

x = reform(x,4,nDimX*nDimY*nDimZ)

return, x  &  end