Source code for fourier

""" Functions for fourier far-field calculations.  Much copied or
adapted from Joe Green's Matlab code @ JPL.

"""

from numpy import array, abs, angle, blackman, concatenate, cos, dot, exp, fft, floor, ones, pi, real, sin, sinc


def MakeRectangle(x,y,parms, *arg, **kw):
    """
    [rectangle] = MakeRectangle (x,y,parms, window_fraction);
    ------------------------------------------------------------------
    
    Make a bandlimited representation of an rectangle function by taking
    the Fourier transform of an appropriately oriented rect function
    
    Inputs:
    -------
    x      = x-coodinate system (made with meshgrid)
    y      = y-coodinate system (made with meshgrid)
    parms[0]  = x-offset
    parms[1]  = y-offset
    parms[2]  = x-diameter
    parms[3]  = y-diameter
    parms[4]  = rotation angle in radians
    
    window_fraction = NOT CURRENTLY WORKING!
                      fraction of fourier domain to be windowed to
                      reduce 'ringing'
    
    Outputs:
    --------
    rectangle   = NxM grided representation of specified rectangle
    where NxM are the size of the x,y input meshgrids
    
    NOTE:  The intergal of the output function is set to equal the 
           area of the spcified rectangle in the units of the x,y 
           meshgrid
           """   

# ------------------------------------------------------------------
# Version: 2008-12-22 14:10 IJC: Converted Matlab to Python
# Version: 26-Jul-2005 - IJC.  Fixed windowing in fourier domain
# Version: 27-Jun-2005 - Ian.J.Crossfield@jpl.nasa.gov
# Version: 30-Dec-2004 - Joseph.J.Green@jpl.nasa.gov
# ------------------------------------------------------------------

    x = array(x, copy=True)
    y = array(y, copy=True)

    xo    = parms[0]
    yo    = parms[1]
    sig_x = parms[2]
    sig_y = parms[3]
    theta = parms[4]

    defaults = dict(verbose=False)
    for key in defaults:
        if (not kw.has_key(key)):
            kw[key] = defaults[key]

    verbose = kw['verbose']
    if len(arg) == 0:
        win_size = 0
    else:
        win_size = arg[0]

    (npix, npiy) = x.shape

    dx = x[0,1]-x[0,0];
    dy = y[1,0]-y[0,0];

    r = abs(x+1j*y);
    a = angle(x+1j*y);

    xn = (r*cos(a + theta)/dx) * sig_y/dy/npix;
    yn = (r*sin(a + theta)/dx) * sig_x/dx/npiy;

    if win_size>0:
        nwx = floor(npix*win_size)
        w = blackman(nwx).reshape(1,nwx)
        w2x = concatenate((w[:,0:nwx/2], ones((1,npix-nwx)), w[:,nwx/2:nwx]), 1)
        nwy = floor(npiy*win_size)
        w = blackman(nwy).reshape(1,nwy)
        w2y = concatenate((w[:,0:nwy/2], ones((1,npiy-nwy)), w[:,nwy/2:nwy]), 1)
        win = dot(w2x.transpose(), w2y)
        if verbose:
            print "nwx>>" + str(nwx)
            print "nwy>>" + str(nwy)
            print "win.shape>>" + str(win.shape)
    else:
        win = 1
   
    if verbose:
        print "npix>>" + str(npix)
        print "xn.shape>>" + str(xn.shape)
 
    rectangle = real(fft.ifftshift(fft.fft2(fft.fftshift(sinc(xn) * sinc(yn) * 
                                            exp((2j*pi*xn)* xo/sig_y) * 
                                            exp((2j*pi*yn)* yo/sig_x) * win))))

    rectangle = rectangle / rectangle.sum();
    rectangle = rectangle * sig_x * sig_y /dx /dy;

    return rectangle