"""
This script tests the Hankel transform used for generating the Fourier
transform of a model image. It compares it with the previous im3shape
approch where an upsampled real-space model image is generated first
and then Fourier transformed. The interface for using one of these
approaches is the 'use_hankel_transform' option. Here, we test it for
the generation of sersics model images in 'i3_sersics_model_image'
function in 'i3_sersics.c'.
"""

import im3shape
import numpy as np 

# Neede for plotting
try:
    import pylab
except ImportError:
    pylab = None

def main(argv):
    """
    Script to test Hankel transform used for model image generation.
    Here we generate a sersic galaxy model with a Moffat PSF and
    compare the generated model images.
    """

    # Create i3_image of certain stamp size
    stamp_size = 100
    i3_galaxy_legacy = im3shape.I3_image(stamp_size, stamp_size)
    i3_galaxy_hankel = im3shape.I3_image(stamp_size, stamp_size)

    # Create i3_options with default values
    i3o = im3shape.I3_options()
    i3o.stamp_size = stamp_size

    # Create i3_parameter_set for sersics model
    i3p = im3shape.I3_parameter_set('sersics')
    i3p.x0 = 50.0
    i3p.y0 = 50.0
    i3p.e1 = 0.10
    i3p.e2 = 0.20
    i3p.radius = 20.0
    i3p.radius_ratio = 1.0
    i3p.bulge_A = 0.1
    i3p.disc_A = 1.0
    i3p.bulge_index = 4.0
    i3p.disc_index = 1.0
    i3p.delta_e_bulge = 0.0
    i3p.delta_theta_bulge = 0.0

    # Create PSF from Moffat parameters
    i3m = im3shape.I3_moffat_psf()
    i3m.beta = 5
    i3m.fwhm = 1.0
    i3m.e1 = 0.0
    i3m.e2 = 0.0
    i3m.x = 50
    i3m.y = 50

    # Create sersics
    i3o.use_hankel_transform = False
    i3_galaxy_legacy.add_real_space_sersic_with_moffat(i3p, i3o, i3m)

    i3o.use_hankel_transform = True
    i3_galaxy_hankel.add_real_space_sersic_with_moffat(i3p, i3o, i3m)

    import pdb
    pdb.set_trace()
    # Compare both images
    if pylab:
        pylab.figure(1)
        pylab.subplot(131)
        pylab.imshow(i3_galaxy_legacy)
        pylab.title('Sersic model generated old-school way')
        pylab.subplot(132)
        pylab.imshow(i3_galaxy_hankel)
        pylab.title('Sersic model with Hankel transform')
        pylab.subplot(133)
        pylab.imshow(i3_galaxy_hankel-i3_galaxy_legacy)
        pylab.title('Difference image')
 
if __name__ == "__main__":
    main(sys.argv)