#!/usr/bin/python
##! /Users/dhw/anaconda/bin/python
# gaussft.py -- compute and plot the FT of a Gaussian function
# gaussft.py  n loc1 sigma1 
#   n = length of array
#   loc1 = location of Gaussian in interval 0-n
#   sigma1 = dispersion of Gaussian 1, in pixels

import sys
import numpy as np
import matplotlib.pyplot as plt
from math import *

n=int(sys.argv[1])
loc1=float(sys.argv[2])
sigma1=float(sys.argv[3])

if (loc1<0 or loc1>n):
    sys.exit("ERROR: locations must be in the range 0-n")

fac=1./sqrt(2.*pi)
x=np.arange(n)

# compute minimum distances assuming a periodic boundary
# This array-based solution was suggesed by Michael Fausnaugh
deltax=np.abs(x-loc1)
mask1 = deltax > n/2
deltax[mask1] = n-deltax[mask1]

# It serves the same function as this loop
#deltax=np.abs(x-loc1)
#for i in range(deltax.size):
#    if (deltax[i]>n/2):
#	    deltax[i]=n-deltax[i]

y=(fac/sigma1)*np.exp(-deltax**2/(2*sigma1**2))

z=np.fft.fft(y)/n
f=np.arange(z.size)
zr=np.real(z)
zi=np.imag(z)

# fper will have the true (negative) frequencies for elements > n/2
fper=np.concatenate((f[f<=n/2],f[f>n/2]-n))
y1=np.exp(-(fper**2)*2*pi**2*(sigma1/n)**2)/n

plt.subplot(2,1,1)
plt.plot(x,y)
plt.plot(x,y,'bo')

plt.subplot(2,1,2)
plt.plot(f,zr,'ro',f,zi,'b^')

# this is the expected FT ONLY if loc1=0
plt.plot(f,y1,'r')

plt.show()