;+
; NAME:
;  sft
; PURPOSE: (one line)
;  slow Fourier transform
; DESCRIPTION:
;  Decompose a function into terms
;  f_t = [f_0, f_1, ...], : complex[n_term+1]
;   s.t., f(p) = f_0 + f_1 exp(i*p) + f_2 exp(2*i*p) + ...
; CATEGORY:
;  Math
; CALLING SEQUENCE:
;  f_t = sft(p, f, nfreq)
; INPUT
;  p = phase, radian, array of length n_phase
;  f = function (real), array of length n_phase
;  nfreq = number of non-zero freq
; OUTPUT
;  return constants f_t[0..nfreq] s.t.
;   f = f_t[0] + f_t[1] * exp(i p) + f_t[2] * exp(i 2 p) + ...
; OPTIONAL OUTPUT
;  f_j = array [nfreq+1, n_phase] s.t.
;    f_j[j, *] = Re( f_t[j] * exp(i j p) )
;  2012 Sep 29. Written Leslie Young SwRI
;  2016 Apr 11 LAY. Modified for inclusion in math IDL library.
;-

function sft, p, f, nfreq, f_j

    _i = complex(0,1)

    n_phase = n_elements(p)

    f_t = complexarr(nfreq+1)  ; the amplitudes
    f_t[0] = mean(f)
    for j = 1, nfreq do f_t[j] = 2*total(f * exp(- _i * j * p ) ) / n_phase
    f_j = fltarr(nfreq+1, n_phase)   ; the values
    for j=0,nfreq do f_j[j,*] = float(f_t[j] * exp(_i*j*p) )

    return, f_t

 end