docs: def createGWB(psr, Amp, gam, noCorr=False, seed=None, turnover=False,
                clm=[N.sqrt(4.0*N.pi)], lmax=0, f0=1e-9, beta=1,
                power=1, userSpec=None, npts=600, howml=10):
      """
      Function to create GW-induced residuals from a stochastic GWB as defined
      in Chamberlin, Creighton, Demorest, et al. (2014).
      
      :param psr: pulsar object for single pulsar
      :param Amp: Amplitude of red noise in GW units 
      :param gam: Red noise power law spectral index
      :param noCorr: Add red noise with no spatial correlations
      :param seed: Random number seed
      :param turnover: Produce spectrum with turnover at frequency f0
      :param clm: coefficients of spherical harmonic decomposition of GW power
      :param lmax: maximum multipole of GW power decomposition
      :param f0: Frequency of spectrum turnover
      :param beta: Spectral index of power spectram for f << f0
      :param power: Fudge factor for flatness of spectrum turnover
      :param userSpec: User-supplied characteristic strain spectrum 
                       (first column is freqs, second is spectrum)
      :param npts: Number of points used in interpolation
      :param howml: Lowest frequency is 1/(howml * T) 
      
      :returns: list of residuals for each pulsar    
      """
  
      if seed is not None:
          N.random.seed(seed)
  
      # number of pulsars
      Npulsars = len(psr)
  
      # gw start and end times for entire data set
      start = N.min([p.toas().min()*86400 for p in psr]) - 86400
      stop = N.max([p.toas().max()*86400 for p in psr]) + 86400
          
      # duration of the signal
      dur = stop - start
      
      # get maximum number of points
      if npts is None:
          # default to cadence of 2 weeks
          npts = dur/(86400*14)
  
      # make a vector of evenly sampled data points
      ut = N.linspace(start, stop, npts)
  
      # time resolution in days
      dt = dur/npts
  
      # compute the overlap reduction function
      if noCorr:
          ORF = N.diag(N.ones(Npulsars)*2)
      else:
          psrlocs = N.zeros((Npulsars,2))
          
          for ii in range(Npulsars):
              if 'RAJ' and 'DECJ' in psr[ii].pars():
                  psrlocs[ii] = N.double(psr[ii]['RAJ'].val), N.double(psr[ii]['DECJ'].val)
              elif 'ELONG' and 'ELAT' in psr[ii].pars():
                  fac = 180./N.pi
                  # check for B name
                  if 'B' in psr[ii].name:
                      epoch = '1950'
                  else:
                      epoch = '2000'
                  coords = ephem.Equatorial(ephem.Ecliptic(str(psr[ii]['ELONG'].val*fac),
                                                           str(psr[ii]['ELAT'].val*fac)),
                                                           epoch=epoch)
                  psrlocs[ii] = float(repr(coords.ra)), float(repr(coords.dec))
  
          psrlocs[:,1] = N.pi/2. - psrlocs[:,1]
          anisbasis = N.array(anis.CorrBasis(psrlocs,lmax)) 
          ORF = sum(clm[kk]*anisbasis[kk] for kk in range(len(anisbasis)))
          ORF *= 2.0
  
      # Define frequencies spanning from DC to Nyquist. 
      # This is a vector spanning these frequencies in increments of 1/(dur*howml).
      f = N.arange(0, 1/(2*dt), 1/(dur*howml))
      f[0] = f[1] # avoid divide by 0 warning
      Nf = len(f)
  
      # Use Cholesky transform to take 'square root' of ORF
      M = N.linalg.cholesky(ORF)
  
      # Create random frequency series from zero mean, unit variance, Gaussian distributions
      w = N.zeros((Npulsars, Nf), complex)
      for ll in range(Npulsars):
          w[ll,:] = N.random.randn(Nf) + 1j*N.random.randn(Nf)
  
      # strain amplitude
      if userSpec is None:
          
          f1yr = 1/3.16e7
          alpha = -0.5 * (gam-3)
          hcf = Amp * (f/f1yr)**(alpha)
          if turnover:
              si = alpha - beta
              hcf /= (1+(f/f0)**(power*si))**(1/power)
  
      elif userSpec is not None:
          
          freqs = userSpec[:,0]
          if len(userSpec[:,0]) != len(freqs):
              raise ValueError("Number of supplied spectral points does not match number of frequencies!")
          else:
              fspec_in = interp.interp1d(N.log10(freqs), N.log10(userSpec[:,1]), kind='linear')
              fspec_ex = extrap1d(fspec_in)
              hcf = 10.0**fspec_ex(N.log10(f))
  
      C = 1 / 96 / N.pi**2 * hcf**2 / f**3 * dur * howml
  
      ### injection residuals in the frequency domain
      Res_f = N.dot(M, w)
      for ll in range(Npulsars):
          Res_f[ll] = Res_f[ll] * C**(0.5)    # rescale by frequency dependent factor
          Res_f[ll,0] = 0			    # set DC bin to zero to avoid infinities
          Res_f[ll,-1] = 0		    # set Nyquist bin to zero also
  
      # Now fill in bins after Nyquist (for fft data packing) and take inverse FT
      Res_f2 = N.zeros((Npulsars, 2*Nf-2), complex)    
      Res_t = N.zeros((Npulsars, 2*Nf-2))
      Res_f2[:,0:Nf] = Res_f[:,0:Nf]
      Res_f2[:, Nf:(2*Nf-2)] = N.conj(Res_f[:,(Nf-2):0:-1])
      Res_t = N.real(N.fft.ifft(Res_f2)/dt)
  
      # shorten data and interpolate onto TOAs
      Res = N.zeros((Npulsars, npts))
      res_gw = []
      for ll in range(Npulsars):
          Res[ll,:] = Res_t[ll, 10:(npts+10)]
          f = interp.interp1d(ut, Res[ll,:], kind='linear')
          res_gw.append(f(psr[ll].toas()*86400))
  
      #return res_gw
      ct = 0
      for p in psr:
          p.stoas[:] += res_gw[ct]/86400.0
          ct += 1