code: def create_fourier_design_matrix(t, nmodes, freq=False, Tspan=None,
                                   logf=False, fmin=None, fmax=None):
      
  
      N = len(t)
      F = np.zeros((N, 2 * nmodes))
  
      if Tspan is not None:
          T = Tspan
      else:
          T = t.max() - t.min()
  
      # define sampling frequencies
      if fmin is not None and fmax is not None:
          f = np.linspace(fmin, fmax, nmodes)
      else:
          f = np.linspace(1 / T, nmodes / T, nmodes)
      if logf:
          f = np.logspace(np.log10(1 / T), np.log10(nmodes / T), nmodes)
  
      Ffreqs = np.zeros(2 * nmodes)
      Ffreqs[0::2] = f
      Ffreqs[1::2] = f
  
      F[:,::2] = np.sin(2*np.pi*t[:,None]*f[None,:])
      F[:,1::2] = np.cos(2*np.pi*t[:,None]*f[None,:])
  
      if freq:
          return F, Ffreqs
      else:
          return F
docs: Construct fourier design matrix from eq 11 of Lentati et al, 2013
  
      :param t: vector of time series in seconds
      :param nmodes: number of fourier coefficients to use
      :param freq: option to output frequencies
      :param Tspan: option to some other Tspan
      :param logf: use log frequency spacing
      :param fmin: lower sampling frequency
      :param fmax: upper sampling frequency
  
      :return: F: fourier design matrix
      :return: f: Sampling frequencies (if freq=True)