code: def sigma_clip(arrin, niter=4, nsig=4, extra={}, verbose=False):
      
      arr = numpy.array(arrin, ndmin=1, copy=False)
  
      index = numpy.arange( arr.size )
  
      for i in numpy.arange(niter):
          m = arr[index].mean()
          s = arr[index].std()
  
          if verbose:
              sys.stdout.write('iter %s\tnuse: %s\tmean %s\tstdev %s\n' % \
                  (i+1, index.size,m,s))
  
          clip = nsig*s
  
          w, = numpy.where( (numpy.abs(arr[index]) - m) < clip )
  
          if w.size == 0:
              sys.stderr.write("nsig too small. Everything clipped on iteration %d" % i+1)
              return m,s
  
  
          index = index[w]
  
      # Calculate final stats
      amean = arr[index].mean()
      asig = arr[index].std()
  
      extra['index'] = index
      return amean, asig
docs: NAME:
        sigma_clip()
        
      PURPOSE:
        Calculate the mean/stdev of an array with sigma clipping. Iterate
        niter times, removing elements that are outside nsig, and recalculating
        mean/stdev.
  
      CALLING SEQUENCE:
        mean,stdev = sigma_clip(arr, niter=4, nsig=4, extra={})
      
      INPUTS:
        arr: A numpy array or a sequence that can be converted.
  
      OPTIONAL INPUTS:
        niter: number of iterations, defaults to 4
        nsig: number of sigma, defaults to 4
  
      OUTPUTS:
        mean,stdev: A tuple containing mean and standard deviation.
      OPTIONAL OUTPUTS
        extra={}: Dictionary containing the array of used indices in
           extra['index']
  
      REVISION HISTORY:
        Converted from IDL: 2006-10-23. Erin Sheldon, NYU