## -*-SM-*- Set SM-mode in emacs
#
# Macros to make dealing with complex numbers for FFT's easier
# Assumes that complex vector `name' is represented by two vectors
# called name_r and name_i
#
fourier		# Help for entire fourier package
  echo Macros to make dealing with complex numbers for FFT's easier
  echo Assumes that complex vector `name' is represented by two vectors
  echo called name_r and name_i. Also assumes that you said
  echo      overload fft 1
  echo
  echo fft	Direct FFT: fft name-x name-k
  echo ifft	Inverse FFT: ifft name-k name-x
  echo cadd	Add complex numbers: \$1 = \$2 + \$3
  echo carg	Argument: \$1 = arg(\$2) (if one arg, return the modulus)
  echo cdiv	Divide complex numbers: \$1 = \$2/\$3
  echo cmod	Modulus: \$1 = |\$2| (if one arg, return the modulus)
  echo cmult	Multiply complex numbers: \$1 = \$2*\$3
  echo cset     Set a complex number: \$1 = \$2 or \$1 = \$2 + i\$3
  echo csub	Subtract complex numbers: \$1 = \$2 - \$3
  echo Im	Imaginary part: Im(\$1)
  echo imag	Imaginary part: \$1 = Im(\$2) (if one argument, return imag part)
  echo Re	Real part: Re(\$1)
  echo real	Real part: \$1 = Re(\$2) (if one argument, return real part)
  echo vcenter	an alias for vcentre
  echo vcentre	shift a vector so that its 0th element appears in the middle
  echo vshift	shift a vector by a specified result
fft	2	# Direct FFT: fft name-x name-k
		# Needs overloading to work: OVERLOAD FFT 1
		if(!is_vector($1_r) && !is_vector($1_i) && is_vector($1)) {
		   local set $1_r = $1  local set $1_i = 0*$1
		}
		FFT 1 $1_r $1_i $2_r $2_i
ifft	2	# Inverse FFT: ifft name-k name-x
		if(!is_vector($1_r) && !is_vector($1_i) && is_vector($1)) {
		   local set $1_r = $1  local set $1_i = 0*$1
		}
		FFT -1 $1_r $1_i $2_r $2_i
cadd    3       # Add complex numbers: $1 = $2 + $3
		SET $1_r = $2_r + $3_r
		SET $1_i = $2_i + $3_i
carg	12	# Argument: $1 = arg($2) (if one arg, return the modulus)
		if(!$?2) {
		   define 2 $1 define 1 $0
		}
		SET $1 = atan2($2_i, $2_r)
cdiv    3       # Divide complex numbers: $1 = $2/$3
		SET _$0 = $3_r**2 + $3_i**2
		SET $1_r = ($2_r*$3_r + $2_i*$3_i)/_$0
		SET $1_i = ($2_i*$3_r - $2_r*$3_i)/_$0
		DELETE _$0
cmod	12	# Modulus: $1 = |$2| (if one arg, return the modulus)
		if(!$?2) {
		   define 2 $1 define 1 $0
		}
		SET $1 = SQRT($2_r**2 + $2_i**2)
cmult	3	# Multiply complex numbers: $1 = $2*$3
		SET $1_r = $2_r*$3_r - $2_i*$3_i
		SET $1_i = $2_i*$3_r + $2_r*$3_i
cset	23	# Set one complex number equal to another or a constant
		if($?3) {
		   SET $1_r = $2 + 0*$3
		   SET $1_i = 0*$2 + $3
		} else {
		   SET $1_r = $2_r
		   SET $1_i = $2_i
		}
csub    3       # Subtract complex numbers: $1 = $2 - $3
		SET $1_r = $2_r - $3_r
		SET $1_i = $2_i - $3_i
Im	11	# an alias for imag
		set $0 = imag($1)
imag	12	# Imaginary part: $1 = Im($2) (if one argument, return imag part)
		if(!$?2) {
		   define 2 $1 define 1 $0
		}
		SET $1 = $2_i
Re	11	# an alias for real
		set $0 = real($1)
real	12	# Real part: $1 = Re($2) (if one argument, return real part)
		if(!$?2) {
		   define 2 $1 define 1 $0
		}
		SET $1 = $2_r
vcenter	11	## an alias for vcentre
		set $0 = vcentre($1)
vcentre	11	# shift a vector so that its 0th element appears in the middle;
		# this is nice for looking at power spectra: e.g.
		#  fft 1 y fy   cmod p fy   set p=p*p
		#  lim x p box con x (vcentre(p))
		local set i=0,dimen($1)-1
		local set rx=$1 if(i < dimen($1)/2)
		local set lx=$1 if(i >= dimen($1)/2)
		set $0=lx concat rx
		#
vshift	22	# shift a vector $1 by $2 indices; like vcentre but
		# no periodic boundary conditions
		local set i=0,dimen($1)-1
		local set rx=$1 if(i < $2)
		set $0= 0*do(1,$2) concat rx
		#
