;+
;
; Project   : Densitycube, reconstruction, 3D, and s3drs software
;                   
; Name      : cube_sphere
;               
; Purpose   : used to manipulate 3D densitycubes and do 3D reconstructions
;               
; Explanation: creates a 3D spherical mask cube for a given radius
;               
; Use       : standalone or via s3drs GUI
;    
; Inputs    : N, r_inner, r_outer, center
;               
; Outputs   : returns the 3D boolean mask cube
;
; Keywords  : 
;               
; Common    : none
;               
; Restrictions: none
;               
; Side effects: none
;               
; Category    : 3d, density, cube
;               
; Prev. Hist. : None.
;
; Written     : Sandy Antunes, NRC, March-April 2009
;               
;-            
; creates a spherical mask/cube for a given size N with region
; 'radius'.
; usage:
;   mask3d = cube_sphere(64,4)

; to mask out an existing ne density cube, try something like this:
;   N=64; e.g.
;   m=p_getmask('COR2',N,rpixels=rpixels)
;   mask3D = cube_sphere(N,rpixels[0],r_outer=rpixels[1])
;   img=img*mask3D
;
; To carve out a typical pixon cube, you can also use:
;   mask3D=cube_sphere(pxwrapper.sr.ncube,r_occ/pxwrapper.sr.L0tau[0],$
;          r_outer=pxwrapper.sr.ncube/2)
;

FUNCTION cube_sphere,N,r_inner,r_outer=r_outer,center=center

  if (n_elements(N) eq 0) then begin
      print,'usage is: mask3D = cube_los(N,r_inner,[r_outer=r_outer],[center=center])'
      return,-1
  end
  if (n_elements(center) eq 0) then center=[(N-1)/2.0,(N-1)/2.0,(N-1)/2.0]

  ; do outside radius as a recursive call and invert
  if (n_elements(r_outer) ne 0) then begin
      mask3D = ~ (cube_sphere(N,r_outer,center=center))
  end else begin
      mask3D = bytarr(N,N,N) + 1
  end

  ; MAIN LOOP
  if (n_elements(r_inner) ne 0) then begin

      limits = [center - r_inner,center + r_inner]
      for i=0,2 do limits[i]=MAX([limits[i],0])
      for i=3,5 do limits[i]=MIN([limits[i],N])

      for x=limits[0],limits[3]-1 do $
        for y=limits[1],limits[4]-1 do $
        for z=limits[2],limits[5]-1 do $
        if (sqrt((x-center[0])^2+(y-center[1])^2+(z-center[2])^2) le r_inner) $
        then mask3d[x,y,z]=0

  end

  return,mask3D

END