pro VertexSkirt, Radius, Height, Vertex, Polygon, angles=Angles,	$
		closetop=CloseTop, closebottom=CloseBottom, inside=Inside
;+
; NAME:
;	VertexSkirt
; PURPOSE:
;	Set up vertex and polygon array for a cylindrically symmetric surface
; CATEGORY:
; CALLING SEQUENCE:
;	VertexSkirt, Radius, Height, Vertex, Polygon, angles=Angles,	$
;		closetop=CloseTop, closebottom=CloseBottom, inside=Inside
; INPUTS:
;	Radius		scalar, array[n]
;					radii of surface at various heights
;	Height		scalar, array[n]
;					set of heights along the symmetry axis
;					Radius and Height should be arrays of equal size
;					if Radius and Height are scalars then
;						the arrays [Radius,Radius] and [0,Height] are used
;					if Radius is a scalar and Height an array then
;						Radius = replicate(Radius, n_elements(Height)
;					if Height is a scalar and Radius and array then
;						Height = Height*findgen(n)/(n-1)
;						(i.e. n equal steps between 0 and Height)
;					if the arrays are not of equal size then the longer
;						array is truncated.
; OPTIONAL INPUT PARAMETERS:
;	angles=Angles	array[k]
;					angles (radians)
;					if not set then Angles = [0,1,2,...,360] degrees is used;
;					at each Radius/Height combination vertices are put at
;					locations [Radius*cos(Angle),Radius*sin(Angle),Height]
;	/closetop		if set, a polygon for the top of the surface
;					at Height(n-1),Radius(n-1) is added
;	/closebottom	if set, a polygon for the bottom of the surface
;					at Height(0),Radius(0) is added
;	/inside			if set, the polygons describe the inside of the
;					the surface (i.e. the vertices are ordered in clockwise
;					rather than counter-clockwise direction)
; OUTPUTS:
;	Vertex		float array[3,m1]
;	Polygon		float array[m2]
;					vertex and polygon arrays as required by the Polyshade
;					procedure. The symmetry axis of the surface is used as
;					the z-axis.
; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
;>	At each height vertices are taken at
;	[Radius*cos(Angle),Radius*sin(Angle),Height]
;>	Pairs of vertices at neigboring heights are combined to form rectangular
;	polygons, or triangular polygons (if Radius=0 at one of the heights).
;>	Note that it is possible to produce non-cylindrically symmetric surfaces,
;	e.g. using a scalar Radius=Height and Angle=[0,90,180,270,360], will
;	produce a cube.
; MODIFICATION HISTORY:
;	June 1997, Paul Hick (UCSD)
;-
on_error, 2

Inside  = fix( keyword_set(Inside) )	; Don't remove fix call !!!!!
Outside = 1-Inside

if n_elements(Angles) le 1 then A = 2*!pi*(findgen(361)/360) else A = Angles
nA = n_elements(A)

nH = n_elements(Height)
nR = n_elements(Radius)

if nH*nR eq 1 then begin
	R = [Radius(0),Radius(0)]
	H = [0,Height(0)]
	nR = 2
endif else if nH eq 1 then begin
	R = Radius
	H = Height(0)*findgen(nR)/(nR-1)
endif else if nR eq 1 then begin
	R = replicate(Radius(0),nH)
	H = Height
	nR = nH
endif else begin
	R = Radius*replicate(1,nH)
	H = Height*replicate(1,nR)
	nR = nR < nH
endelse
nH = nR										; Now nR = nH

CloseTop	= keyword_set(CloseTop)		and R(nR-1) ne 0
CloseBottom = keyword_set(CloseBottom)	and R(0   ) ne 0

nT = 0
if CloseBottom then nT = nT+1				; Polygon to close top
if CloseTop    then nT = nT+1				; Polygon to close bottom

; The 'skirt' is circumscribed by a cylinder with axis of symmetry
; in the z-direction. The 'skirt' makes contact with the cylinder
; at points in the x-y plane with phase angles given in array Angles.
; Each polygon is a vertical rectangle

n3 = 0
n4 = 0
nV = 0										; # vertices

for iR=0,nR-1 do begin						; Loop over nR-1 pairs of radii
	case R(iR) of
	0: begin         						; Zero radius: R(iR)=0
		nV = nV+1   						; Only one additional vertex
		if iR+1 lt nR then if R(iR+1) ne 0 then n3 = n3+1  ; Row of triangular polygons
	end
	else: begin								; Non-zero R(iR)
		nV = nV+nA							; nA vertices
		if iR+1 lt nR then begin
 	   		case R(iR+1) of
    		0:		n3 = n3+1       		; Row of triangular polygons
    		else:	n4 = n4+1				; Row of rectangular polygons
    		endcase
    	endif
	end
	endcase
endfor

n3 = n3*(nA-1)								; # triangular polygons
n4 = n4*(nA-1) 								; # rectangular polygons

Vertex  = fltarr(3,nV)						; List of all vertices
Polygon = lonarr(4*n3+5*n4+(1+nA)*nT)

cA = cos(A)									; Auxilliary arrays
sA = sin(A)
rA = replicate(1,nA)
lA1 = lindgen(nA-1)
lA = [lA1,nA-1]

iV = -1  	        	               		; Last vertex added (second index of Vertex)
iP = -1     	                            ; Last index filled in Polygon
iR = 0										; Next radius to be processed
AddVertex = 1B

while iR lt nR-1 do begin					; Loop over all pairs of neighboring radii
	if R(iR) eq 0 then begin
		if R(iR+1) ne 0 then begin			; Zero followed by non-zero radius

			if AddVertex then begin
        		iV = iV+1
        		Vertex(*,iV) = [0,0,H(iR)]	; Add vertex on z-axis
        	endif

			i = iR+1
        	Vertex(*,iV+1:iV+nA) = transpose([ [R(i)*cA], [R(i)*sA], [H(i)*rA] ])

        	i = iP+1+4*lA1
			Polygon(i  ) = 3				; # pnts in triangular polygon
			Polygon(i+1) = iV				; Points to vertex on z-axis
			Polygon(i+2) = iV+1+Outside+lA1
			Polygon(i+3) = iV+1+Inside +lA1

			iV = iV+nA                  	; Added nA vertices
			iP = iP+4*(nA-1)				; Added nA-1 triangular polygons

		endif

		iR = iR+1
		AddVertex = R(iR) eq 0

	endif else begin						; Non-zero R(iR)

		i0 = (where(R(iR:*) eq 0))(0)		; Find next zero radius
		if i0 eq 1 then begin            	; Non-zero followed by zero radius

			if AddVertex then begin
  	      		Vertex(*,iV+1:iV+nA) = transpose([ [R(iR)*cA], [R(iR)*sA], [H(iR)*rA] ])
  	      		iV = iV+nA
        	endif
        	iV = iV+1
        	Vertex(*,iV) = [0,0,H(iR+1)]	; Add vertex on z-axis

        	i = iP+1+4*lA1
			Polygon(i  ) = 3				; # pnts in triangular polygon
			Polygon(i+1) = iV				; Points to vertex on z-axis
			Polygon(i+2) = iV-nA+Inside +lA1
			Polygon(i+3) = iV-nA+Outside+lA1

			iP = iP+4*(nA-1)				; Added nA-1 triangular polygons

		endif else begin
			if i0 eq -1 then i0 = nR-iR
			i0 = i0-1						; Always i0 > 0; iR:iR+i0 is range of non-zero radii

            if AddVertex then i = 0 else i = 1	; Drop first row of vertices

			n = nA*(i0-i+1)
			Vertex(0,iV+1:iV+n) = cA#R(iR+i:iR+i0)
			Vertex(1,iV+1:iV+n) = sA#R(iR+i:iR+i0)
			Vertex(2,iV+1:iV+n) = rA#H(iR+i:iR+i0)

			iV = iV-i*nA

			i = iP+1+5*lindgen( (nA-1)*i0 )
			Polygon(i) = 4					; 4 vertices per rectangular polygon
			i = i+1
			Polygon(i  ) = iV+1+replicate(nA,nA-1)#lindgen(i0)+lA1#replicate(1,i0)
			Polygon(i+1) = nA*Inside +1*Outside+Polygon(i)
			Polygon(i+2) = nA        +1        +Polygon(i)
			Polygon(i+3) = nA*Outside+1*Inside +Polygon(i)

			iV = iV+nA*(i0+1)
			iP = iP+5*(nA-1)*i0

		endelse

		iR = iR+i0
		AddVertex = 0B

	endelse

endwhile

if CloseBottom or CloseTop then lA = (nA-1)*Outside+(Inside-Outside)*lA

if CloseBottom then begin					; Close bottom (R(0)) of surface
	i = iP+1
	Polygon(i) = nA
	Polygon(i+1:i+nA) = lA
	iP = iP+1+nA
endif

if CloseTop then begin						; Close top (R(nR-1)) of surface
	i = iP+1
	Polygon(i) = nA
	Polygon(i+1:i+nA) = nV-1-lA
endif

return  &  end