;+
;NAME:
;genxyzref:  generate an xyz spherical grid with random radial errors
;SYNTAX:n=genrefxyz(gridlen,cen,radius,diameter,rms,xyzN,grid=grid,nodata=nodata,radialErr=radialErr)
;ARGS:
;gridlen : long  number of points in x or y dimension
;cen[3]  : dblarr center of curvature for sphere
;radius  : double radius of curvature of sphere
;diameter: double diameter across top part of spherical section.. for primary this would be 1000 ft
;rms     : double rms for radial errors. if 0 then no noise entered.
;                      must be in same units as cen,radius, maxx
;KEYWORDS:
; grid   :        if set then return an  (3,gridlen,gridlen) array. 
;                 points outside the sphere will contain nodata in the z value
;                also return radialerr(gridlen,gridlen) will be the offset from the reference
;                radius and center. they will also have nodata in the locations outside of the sphere
;                 outside of the sphere will contain nodata=nodata value
;nodata :double   if /grid supplied then points outside of the sphere will contain
;                 the no data value in their z component, def. -999
;RETURNS:
;n        : long  number of points in grid. the grid will be limited to points that sit on the sphere
;xyzN[3,n]: dblarr   grid of points. 
;radialErr[gridlen,gridlen] dblarr if /grid is supplied and rms ne 0 then this will hold the radial error
;                    for each x,y point. if outside the sphere it will contain the nodata value 
;
;DESCRIPTION:
;  This is used to generate model data representing the primary reflector.
;The primary has
; 1000 ft diameter
; 870  ft radius of curvature
; 0,0,870 : center of sphere 
;
;To duplicate this data you would use
; 1000ft diameter --> diameter =1000
;Typically the data will be centered at ao9
;  this would put the center of curvature at 0,0,870 (or maybe some small offset from ao9)
;gridlen: number of points across the diameter
;if at ao9 then cen=[0,0,870]
;
;rms  ne will add random noise to the radius with an rms of rms.
;
; This routine will work in ft or meters.. but you must be consistant with all 
; the arguments.
;-
function genrefxyz,gridlen,cen,radius,diameter,rms,xyzN,grid=grid,nodata=nodata,radialErr=radialErr

	nodata=(n_elements(nodata) eq 1)?nodata: -999d
	usegrid = keyword_set(grid)
	d=(dindgen(gridlen)/gridLen - .5)*diameter
	xyzN=dblarr(3,gridlen,gridlen)
	for i=0,gridlen-1 do begin &$
		xyzN[0,*,i]=d &$
		xyzN[1,i,*]=d &$
	endfor
	nn=(1l*gridlen)*gridlen
	xyzN=reform(xyzN,3,nn)
;
; limit points to xyradius if not usegrid
;
	xyRad=reform(sqrt(xyzN[0,*]^2 + xyzN[1,*]^2))
	maxX=diameter/2d
	ii= where(xyRad le maxX,nn)
	if not usegrid then begin
		xyzN=xyzN[*,ii]
	endif
;
;   now compute z for each point
;                                         since zcen above z this makes the sqrt sign work
;    use R^2 = (x -xcen)^2 + (y-ycen)^2 + (zcen-z)^2
;
;    zcen-z=sqrt(r^2 - (x-xcen)^2 - (y-ycen)^2)
;
	xyzN[2,*]=cen[2]  -  sqrt(Radius^2 - ((reform(xyzN[0,*])-cen[0])^2 + (reform(xyzN[1,*])-cen[1])^2)) 
   
;
;   see if they want noise 
	if (rms ne 0) then begin
;  go to spherical coord, add noise to radius then come back
;
		rtp=blmxyztosph(xyzN,center=cen)
		rtp[0,*]+=(randomN(0,nn,/double,/normal)*rms)
		xyzN=blmsphtoxyz(rtp,cent=cen)
	endif
	if usegrid then begin
		if rms eq 0. then rtp=blmxyztosph(xyzN,center=cen)
	    ii= where(xyRad gt maxX,nn)
		xyzN[2,ii]=nodata
		xyzN=reform(xyzN,3,gridlen,gridlen)
		radialErr=reform(rtp[0,*] - radius)
		radialErr[ii]=nodata
		radialErr=reform(radialErr,gridlen,gridlen)
	endif
	return,nn
end