test_chromatic.py

# Copyright (c) 2012-2026 by the GalSim developers team on GitHub
# https://github.com/GalSim-developers
#
# This file is part of GalSim: The modular galaxy image simulation toolkit.
# https://github.com/GalSim-developers/GalSim
#
# GalSim is free software: redistribution and use in source and binary forms,
# with or without modification, are permitted provided that the following
# conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
#    list of conditions, and the disclaimer given in the accompanying LICENSE
#    file.
# 2. Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions, and the disclaimer given in the documentation
#    and/or other materials provided with the distribution.
#

import copy
import os
import copy
import numpy as np

import galsim
from galsim_test_helpers import *

bppath = os.path.join(galsim.meta_data.share_dir, "bandpasses")
sedpath = os.path.join(galsim.meta_data.share_dir, "SEDs")

# from pylab import *
# def plotme(image):
#     imshow(image.array)
#     show()

# liberal use of globals here...
zenith_angle = 20 * galsim.degrees
R500 = galsim.dcr.get_refraction(500.0, zenith_angle) # normalize refraction to 500nm

# some profile parameters to test with
bulge_n = 4.0
bulge_hlr = 0.5
bulge_e1 = 0.2
bulge_e2 = 0.2

disk_n = 1.0
disk_hlr = 1.0
disk_e1 = 0.4
disk_e2 = 0.2

PSF_hlr = 0.3
PSF_beta = 3.0
PSF_e1 = 0.01
PSF_e2 = 0.06

shear_g1 = 0.01
shear_g2 = 0.02

# load a filter
bandpass = (galsim.Bandpass(os.path.join(bppath, 'LSST_r.dat'), 'nm'))
bandpass_g = (galsim.Bandpass(os.path.join(bppath, 'LSST_g.dat'), 'nm'))
bandpass_z = (galsim.Bandpass(os.path.join(bppath, 'LSST_z.dat'), 'nm'))

# load some spectra
bulge_SED = (galsim.SED(os.path.join(sedpath, 'CWW_E_ext.sed'), wave_type='ang',
                        flux_type='flambda')
             .thin(rel_err=1e-3)
             .withFluxDensity(target_flux_density=0.3, wavelength=500.0))

disk_SED = (galsim.SED(os.path.join(sedpath, 'CWW_Sbc_ext.sed'), wave_type='ang',
                       flux_type='flambda')
            .thin(rel_err=1e-3)
            .withFluxDensity(target_flux_density=0.3, wavelength=500.0))

# define the directory containing some reference images
refdir = os.path.join(".", "chromatic_reference_images") # Directory containing the reference


@timer
def test_draw_add_commutativity():
    """Compare two chromatic images, one generated by adding up GSObject profiles before drawing,
    and one generated (via galsim.chromatic) by drawing image summands wavelength-by-wavelength
    while updating the profile and adding as you go.
    """
    import time

    stamp_size = 32
    pixel_scale = 0.2

    #------------------------------------------------------------------------------
    # Use galsim.base functions to generate chromaticity by creating an effective
    # PSF by adding together weighted monochromatic PSFs.
    # Profiles are added together before drawing.
    #------------------------------------------------------------------------------

    # make galaxy
    GS_gal = galsim.Sersic(n=bulge_n, half_light_radius=bulge_hlr)
    GS_gal = GS_gal.shear(e1=bulge_e1, e2=bulge_e2)
    GS_gal = GS_gal.shear(g1=shear_g1, g2=shear_g2)

    # make effective PSF with Riemann sum midpoint rule
    mPSFs = [] # list of flux-scaled monochromatic PSFs
    N = 50
    h = (bandpass.red_limit * 1.0 - bandpass.blue_limit) / N
    ws = [bandpass.blue_limit + h*(i+0.5) for i in range(N)]
    shift_fn = lambda w:(0, ((galsim.dcr.get_refraction(w, zenith_angle) - R500)
                             * (galsim.radians / galsim.arcsec)))
    dilate_fn = lambda w:(w/500.0)**(-0.2)
    for w in ws:
        flux = bulge_SED(w) * bandpass(w) * h
        mPSF = galsim.Moffat(flux=flux, beta=PSF_beta, half_light_radius=PSF_hlr*dilate_fn(w))
        mPSF = mPSF.withGSParams(maxk_threshold=1.e-4)
        mPSF = mPSF.shear(e1=PSF_e1, e2=PSF_e2)
        mPSF = mPSF.shift(shift_fn(w))
        mPSFs.append(mPSF)
    PSF = galsim.Add(mPSFs)

    # final profile
    final = galsim.Convolve([GS_gal, PSF])
    GS_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    GS_kimage = galsim.ImageCD(stamp_size, stamp_size, scale=pixel_scale)
    t2 = time.time()
    GS_image = final.drawImage(image=GS_image)
    GS_kimage = final.drawKImage(image=GS_kimage)
    t3 = time.time()
    print('GS_object drawImage, drawKImage took {0} seconds.'.format(t3-t2))
    # plotme(GS_image)

    #------------------------------------------------------------------------------
    # Use galsim.chromatic to generate chromaticity.  Internally, this module draws
    # the result at each wavelength and adds the results together.  I.e., drawing
    # and adding happen in the reverse order of the above.
    #------------------------------------------------------------------------------

    # make galaxy
    mono_gal = galsim.Sersic(n=bulge_n, half_light_radius=bulge_hlr)
    chromatic_gal = mono_gal * bulge_SED
    check_pickle(bulge_SED)
    check_pickle(chromatic_gal, lambda x: x.drawImage(bandpass, method='no_pixel',
                                                   nx=10, ny=10, scale=1))
    check_pickle(chromatic_gal)

    # Shear object
    chromatic_gal = chromatic_gal.shear(e1=bulge_e1, e2=bulge_e2)
    chromatic_gal = chromatic_gal.shear(g1=shear_g1, g2=shear_g2)
    check_pickle(chromatic_gal, lambda x: x.drawImage(bandpass, method='no_pixel',
                                                   nx=10, ny=10, scale=1))
    check_pickle(chromatic_gal)

    # make chromatic PSF
    mono_PSF = galsim.Moffat(beta=PSF_beta, half_light_radius=PSF_hlr)
    mono_PSF = mono_PSF.withGSParams(maxk_threshold=1.e-4)
    mono_PSF = mono_PSF.shear(e1=PSF_e1, e2=PSF_e2)
    chromatic_PSF = galsim.ChromaticTransformation(mono_PSF, flux_ratio=1.0)
    check_pickle(chromatic_PSF, lambda x: (x.evaluateAtWavelength(bandpass.effective_wavelength)
                                         .drawImage(method='no_pixel', nx=10, ny=10, scale=1)))
    check_pickle(chromatic_PSF)
    chromatic_PSF = chromatic_PSF.dilate(dilate_fn)
    chromatic_PSF = chromatic_PSF.shift(shift_fn)

    # final profile
    chromatic_final = galsim.Convolve([chromatic_gal, chromatic_PSF])
    chromatic_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    chromatic_kimage = galsim.ImageCD(stamp_size, stamp_size, scale=pixel_scale)
    # use chromatic parent class to draw without ChromaticConvolution acceleration...
    t4 = time.time()
    integrator = galsim.integ.ContinuousIntegrator(galsim.integ.midptRule, N=N, use_endpoints=False)
    # NB. You cannot use ChromaticObject.drawImage() here, since it will automatically farm out to
    #     the ChromaticConvolution version of drawImage rather than respecting the
    #     ChromaticObject specification.  Using super() doesn't seem to work either.  So I just
    #     went ahead and converted this statement to the new format.  There are a couple other
    #     similar times in the test suite where we want to force it to use the base class
    #     implementation, so those had to be switched as well.
    galsim.ChromaticObject.drawImage(chromatic_final, bandpass, image=chromatic_image,
                                     integrator=integrator, add_to_image=True)
    galsim.ChromaticObject.drawKImage(chromatic_final, bandpass, image=chromatic_kimage,
                                      integrator=integrator)
    t5 = time.time()
    print('ChromaticObject drawImage, drawKImage took {0} seconds.'.format(t5-t4))
    # plotme(chromatic_image)

    peak = chromatic_image.array.max()
    printval(GS_image, chromatic_image)
    np.testing.assert_array_almost_equal(
        chromatic_image.array/peak, GS_image.array/peak, 6,
        err_msg="Directly computed chromatic image disagrees with image created using "
                +"galsim.chromatic")

    kpeak = chromatic_kimage.array.real.max()
    np.testing.assert_array_almost_equal(
        chromatic_kimage.array/kpeak, GS_kimage.array/kpeak, 6,
        err_msg="Directly computed chromatic kimage disagrees with kimage created using "
                +"galsim.chromatic")

    # Repeat with multiple inseparable profiles.
    delta = galsim.ChromaticObject(galsim.DeltaFunction()).rotate(lambda wave: wave*galsim.degrees)
    chromatic_final2 = galsim.Convolve(chromatic_gal, chromatic_PSF, delta)
    chromatic_final2.drawImage(bandpass, image=chromatic_image, integrator=integrator)
    chromatic_final2.drawKImage(bandpass, image=chromatic_kimage, integrator=integrator)
    # Note: fft vs real space differences now, so only accurate to 1.e-3
    np.testing.assert_array_almost_equal(chromatic_image.array/peak, GS_image.array/peak, 3)
    np.testing.assert_array_almost_equal(chromatic_kimage.array/kpeak, GS_kimage.array/kpeak, 6)

    # Check error handling of too few sample points
    integrator = galsim.integ.ContinuousIntegrator(galsim.integ.midptRule, N=1, use_endpoints=False)
    with assert_raises(ValueError):
        chromatic_final.drawImage(bandpass, integrator=integrator)
    integrator = galsim.integ.ContinuousIntegrator(galsim.integ.trapzRule, N=1, use_endpoints=False)
    with assert_raises(ValueError):
        chromatic_final.drawImage(bandpass, integrator=integrator)

    # As an aside, check for appropriate tests of 'integrator' argument.
    assert_raises(ValueError, chromatic_final.drawImage, bandpass, method='no_pixel',
                  integrator='midp') # minor misspelling
    assert_raises(ValueError, chromatic_final.drawKImage, bandpass,
                  integrator='midp') # minor misspelling
    assert_raises(TypeError, chromatic_final.drawImage, bandpass, method='no_pixel',
                  integrator=galsim.integ.midptRule)
    assert_raises(TypeError, chromatic_final.drawKImage, bandpass,
                  integrator=galsim.integ.midptRule)

    # Can't use base class directly.
    assert_raises(NotImplementedError, galsim.integ.ImageIntegrator)


@timer
def test_ChromaticConvolution_InterpolatedImage():
    """Check that we can interchange the order of integrating over wavelength and convolving for
    separable ChromaticObjects.  This involves storing the results of integrating first in an
    InterpolatedImage.
    """
    pixel_scale = 0.2
    stamp_size = 32

    # stars are fundamentally delta-fns with an SED
    star = galsim.Gaussian(fwhm=1.e-8) * bulge_SED
    mono_PSF = galsim.Gaussian(half_light_radius=PSF_hlr)
    PSF = galsim.ChromaticAtmosphere(mono_PSF, base_wavelength=500.0,
                                     zenith_angle=zenith_angle)

    final = galsim.Convolve(star, PSF)
    image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)

    check_pickle(star)
    check_pickle(PSF)
    check_pickle(final)

    # draw image using speed tricks in ChromaticConvolution.draw
    # For this particular test, need to set iimult=4 in order to pass.
    II_image = final.drawImage(bandpass, image=image, iimult=4)
    II_flux = II_image.array.sum()

    image2 = image.copy()
    # draw image without any speed tricks using ChromaticObject.drawImage
    D_image = galsim.ChromaticObject.drawImage(final, bandpass, image=image2)
    D_flux = D_image.array.sum()

    #compare
    print('Flux when integrating first, convolving second: {0}'.format(II_flux))
    print('Flux when convolving first, integrating second: {0}'.format(D_flux))
    printval(II_image, D_image)
    # This used to work with decimal=5, but is apparently sensitive to the particular details of the
    # bandpass thinning used.  decimal=4 is still good though.
    np.testing.assert_array_almost_equal(
        II_image.array, D_image.array, 4,
        err_msg="ChromaticConvolution draw not equivalent to regular draw")

    # Check flux scaling
    II_image2 = (final * 2.).drawImage(bandpass, image=image, iimult=4)
    II_flux2 = II_image2.array.sum()
    np.testing.assert_array_almost_equal(
        II_flux2, 2.*II_flux, 5,
        err_msg="ChromaticConvolution * 2 resulted in wrong flux.")


@timer
def test_chromatic_add():
    """Test the `+` operator on ChromaticObjects"""
    stamp_size = 32
    pixel_scale = 0.2

    # create galaxy profiles
    mono_bulge = galsim.Sersic(n=bulge_n, half_light_radius=bulge_hlr)
    bulge = mono_bulge * bulge_SED
    bulge = bulge.shear(e1=bulge_e1, e2=bulge_e2)

    mono_disk = galsim.Sersic(n=disk_n, half_light_radius=disk_hlr)
    disk = mono_disk * disk_SED
    disk = disk.shear(e1=disk_e1, e2=disk_e2)

    # test `+` operator
    bdgal = bulge + disk
    bdgal = bdgal.shear(g1=shear_g1, g2=shear_g2)

    # now shear the indiv profiles
    bulge = bulge.shear(g1=shear_g1, g2=shear_g2)
    disk = disk.shear(g1=shear_g1, g2=shear_g2)

    # create PSF
    mono_PSF = galsim.Moffat(beta=PSF_beta, half_light_radius=PSF_hlr)
    mono_PSF = mono_PSF.shear(e1=PSF_e1, e2=PSF_e2)
    chromatic_PSF = galsim.ChromaticAtmosphere(mono_PSF, base_wavelength=500.0,
                                               zenith_angle=zenith_angle)

    # create final profile
    final = galsim.Convolve(bdgal, chromatic_PSF)
    image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    image = final.drawImage(bandpass, image=image)

    check_pickle(bulge)
    check_pickle(disk)
    check_pickle(bdgal)
    check_pickle(chromatic_PSF)
    check_pickle(final)

    bulge_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    bulge_part = galsim.Convolve([bulge, chromatic_PSF])
    bulge_image = bulge_part.drawImage(bandpass, image=bulge_image)
    disk_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    disk_part = galsim.Convolve([disk, chromatic_PSF])
    disk_image = disk_part.drawImage(bandpass, image=disk_image)

    piecewise_image = bulge_image + disk_image
    print('bulge image flux: {0}'.format(bulge_image.array.sum()))
    print('disk image flux: {0}'.format(disk_image.array.sum()))
    print('piecewise image flux: {0}'.format(piecewise_image.array.sum()))
    print('bdimage flux: {0}'.format(image.array.sum()))
    printval(image, piecewise_image)
    np.testing.assert_array_almost_equal(
            image.array, piecewise_image.array, 6,
            err_msg="`+` operator doesn't match manual image addition")

    # Check flux scaling
    flux = image.array.sum()
    image = (final * 2.).drawImage(bandpass, image=image)
    flux2 = image.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="ChromaticConvolution with sum * 2 resulted in wrong flux.")

    # apply flux scaling to ChromaticSum
    final2 = galsim.Convolve(bdgal*2, chromatic_PSF)
    image = final2.drawImage(bandpass, image=image)
    flux2 = image.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="ChromaticSum * 2 resulted in wrong flux.")

    # also check that a - b  ==  a + (-b)
    c = bulge_part - disk_part
    d = bulge_part + (-disk_part)
    assert c == d

    cimage = c.drawImage(bandpass, nx=stamp_size, ny=stamp_size, scale=pixel_scale)
    dimage = d.drawImage(bandpass, cimage.copy())
    np.testing.assert_equal(cimage, dimage, "chromatic a-b != a+(-b)")

@timer
def test_dcr_moments():
    """Check that zenith-direction surface brightness distribution first and second moments obey
    expected behavior for differential chromatic refraction when comparing objects drawn with
    different SEDs."""
    stamp_size = 256
    pixel_scale = 0.025

    # stars are fundamentally delta-fns with an SED
    star1 = galsim.Gaussian(fwhm=1.e-8) * bulge_SED
    star2 = galsim.Gaussian(fwhm=1.e-8) * disk_SED

    shift_fn = lambda w:(0, ((galsim.dcr.get_refraction(w, zenith_angle) - R500)
                             * (galsim.radians / galsim.arcsec)))
    mono_PSF = galsim.Moffat(beta=PSF_beta, half_light_radius=PSF_hlr)
    PSF = galsim.ChromaticObject(mono_PSF).shift(shift_fn)

    final1 = galsim.Convolve([star1, PSF])
    final2 = galsim.Convolve([star2, PSF])

    image1 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    image2 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)

    image1 = final1.drawImage(bandpass, image=image1)
    image2 = final2.drawImage(bandpass, image=image2)
    # plotme(image1)

    mom1 = galsim.utilities.unweighted_moments(image1)
    mom2 = galsim.utilities.unweighted_moments(image2)
    dR_image = (mom1['My'] - mom2['My']) * pixel_scale
    dV_image = (mom1['Myy'] - mom2['Myy']) * (pixel_scale)**2

    # analytic moment differences
    R_bulge, V_bulge = bulge_SED.calculateDCRMomentShifts(bandpass, zenith_angle=zenith_angle)
    R_disk, V_disk = disk_SED.calculateDCRMomentShifts(bandpass, zenith_angle=zenith_angle)
    dR_analytic = (R_bulge[1] - R_disk[1]) * 180.0/np.pi * 3600
    dV_analytic = (V_bulge[1,1] - V_disk[1,1]) * (180.0/np.pi * 3600)**2

    # also compute dR_analytic using ChromaticObject.calculateCentroid()
    centroid1 = final1.calculateCentroid(bandpass)
    centroid2 = final2.calculateCentroid(bandpass)
    dR_centroid = (centroid1 - centroid2).y

    print('image delta R:    {0}'.format(dR_image))
    print('analytic delta R: {0}'.format(dR_analytic))
    print('centroid delta R: {0}'.format(dR_centroid))
    print('image delta V:    {0}'.format(dV_image))
    print('analytic delta V: {0}'.format(dV_analytic))
    np.testing.assert_almost_equal(dR_image, dR_analytic, 5,
                                   err_msg="dRbar Shift from DCR doesn't match analytic formula")
    np.testing.assert_almost_equal(dR_analytic, dR_centroid, 10,
                                   err_msg="direct dRbar calculation doesn't match"
                                           +" ChromaticObject.calculateCentroid()")
    np.testing.assert_almost_equal(dV_image, dV_analytic, 5,
                                   err_msg="dV Shift from DCR doesn't match analytic formula")


@timer
def test_chromatic_seeing_moments():
    """Check that surface brightness distribution second moments obey expected behavior
    for chromatic seeing when comparing stars drawn with different SEDs."""
    pixel_scale = 0.0075
    stamp_size = 1024

    # stars are fundamentally delta-fns with an SED
    star1 = galsim.Gaussian(fwhm=1e-8) * bulge_SED
    star2 = galsim.Gaussian(fwhm=1e-8) * disk_SED

    indices = [-0.2, 0.6, 1.0]
    for index in indices:

        mono_PSF = galsim.Gaussian(half_light_radius=PSF_hlr)
        PSF = galsim.ChromaticObject(mono_PSF).dilate(lambda w:(w/500.0)**index)

        final1 = galsim.Convolve([star1, PSF])
        final2 = galsim.Convolve([star2, PSF])

        image1 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
        image2 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)

        image1 = final1.drawImage(bandpass, image=image1)
        image2 = final2.drawImage(bandpass, image=image2)

        shape1 = galsim.utilities.unweighted_shape(image1)
        shape2 = galsim.utilities.unweighted_shape(image2)
        dr2byr2_image = (shape1['rsqr'] - shape2['rsqr']) / shape1['rsqr']

        # analytic moment differences
        r2_1 = bulge_SED.calculateSeeingMomentRatio(bandpass, alpha=index)
        r2_2 = disk_SED.calculateSeeingMomentRatio(bandpass, alpha=index)
        dr2byr2_analytic = (r2_1 - r2_2) / r2_1

        np.testing.assert_almost_equal(dr2byr2_image, dr2byr2_analytic, 5,
                                       err_msg="Moment Shift from chromatic seeing doesn't"+
                                               " match analytic formula")

        print('image delta(r^2) / r^2:    {0}'.format(dr2byr2_image))
        print('analytic delta(r^2) / r^2: {0}'.format(dr2byr2_analytic))


@timer
def test_monochromatic_filter():
    """Check that ChromaticObject drawn through a very narrow band filter matches analogous
    GSObject.
    """
    pixel_scale = 0.2
    stamp_size = 32

    chromatic_gal = galsim.Gaussian(fwhm=1.0) * bulge_SED
    GS_gal = galsim.Gaussian(fwhm=1.0)

    shift_fn = lambda w:(0, ((galsim.dcr.get_refraction(w, zenith_angle) - R500)
                             * (galsim.radians / galsim.arcsec)))
    dilate_fn = lambda wave: (wave/500.0)**(-0.2)
    mono_PSF = galsim.Gaussian(half_light_radius=PSF_hlr)
    mono_PSF = mono_PSF.shear(e1=PSF_e1, e2=PSF_e2)
    chromatic_PSF = galsim.ChromaticObject(mono_PSF).dilate(dilate_fn).shift(shift_fn)

    chromatic_final = galsim.Convolve([chromatic_gal, chromatic_PSF])

    fws = [350, 475, 625, 750, 875, 975] # approximate ugrizy filter central wavelengths
    for fw in fws:
        chromatic_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
        narrow_filter = galsim.Bandpass(galsim.LookupTable([fw-0.01, fw, fw+0.01],
                                                           [1.0, 1.0, 1.0],
                                                           interpolant='linear'), 'nm')
        chromatic_image = chromatic_final.drawImage(narrow_filter, image=chromatic_image)
        # take out normalization
        chromatic_image /= 0.02
        chromatic_image /= bulge_SED(fw)

        # now do non-chromatic version
        GS_PSF = galsim.Gaussian(half_light_radius=PSF_hlr)
        GS_PSF = GS_PSF.shear(e1=PSF_e1, e2=PSF_e2)
        GS_PSF = GS_PSF.dilate(dilate_fn(fw))
        GS_PSF = GS_PSF.shift(shift_fn(fw))
        GS_final = galsim.Convolve([GS_gal, GS_PSF])
        GS_image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
        GS_final.drawImage(image=GS_image)
        # plotme(GS_image)

        printval(chromatic_image, GS_image)
        np.testing.assert_array_almost_equal(chromatic_image.array, GS_image.array, 5,
                err_msg="ChromaticObject.drawImage() with monochromatic filter doesn't match"+
                        "GSObject.drawImage()")

@timer
def test_monochromatic_sed(run_slow):
    # Similar to the above test, but this time with a broad bandpass and a narrow sed.

    bandpass = galsim.Bandpass(galsim.LookupTable([500,1000], [1,1], 'linear'), wave_type='nm')
    flux = 1.e6
    rng = galsim.BaseDeviate(1234)
    pixel_scale = 0.01

    gal_achrom = galsim.Sersic(n=2.8, half_light_radius=0.03, flux=flux)

    diam = 3.1 # meters
    obscuration = 0.11
    nstruts = 5
    aberrations = np.array([0,0,0,0, 0.02, -0.05, -0.15, -0.02, 0.13, 0.06, -0.09, 0.11])

    if run_slow:
        wave_list = [515, 690, 900]
    else:
        wave_list = [515]
    for wave in wave_list:
        # First do the achromatic version at the given wavelength.
        psf_achrom = galsim.OpticalPSF(lam=wave, diam=diam,
                                       aberrations=aberrations, obscuration=obscuration,
                                       nstruts=nstruts)
        obj_achrom = galsim.Convolve(gal_achrom, psf_achrom)
        im1 = galsim.ImageD(50, 50, scale=pixel_scale)
        print('obj_achrom = ',obj_achrom)
        print('im1 = ',im1)
        obj_achrom.drawImage(image=im1)
        print('im1.max,sum = ', im1.array.max(), im1.array.sum())

        # Next do this chromatically using an emission line SED.
        psf_chrom = galsim.ChromaticOpticalPSF(lam=wave, diam=diam,
                                               aberrations=aberrations, obscuration=obscuration,
                                               nstruts=nstruts)
        sed = galsim.EmissionLine(wave)
        gal_chrom = (gal_achrom * sed).withFlux(flux, bandpass=bandpass)
        obj_chrom = galsim.Convolve(gal_chrom, psf_chrom)
        im2 = obj_chrom.drawImage(bandpass, image=im1.copy())
        print('im2.max,sum = ', im2.array.max(), im2.array.sum())
        print('max diff/flux = ',np.max(np.abs(im1.array-im2.array)/flux))
        np.testing.assert_allclose(im2.array/flux, im1.array/flux, rtol=1.e-4, atol=1.e-5)

        # Now achromatic with phot
        im3 = obj_achrom.drawImage(image=im1.copy(), method='phot', rng=rng)
        print('im3.max,sum = ', im3.array.max(), im3.array.sum())
        print('max diff/flux = ',np.max(np.abs(im1.array-im3.array)/flux))
        np.testing.assert_allclose(im3.array/flux, im1.array/flux, atol=3.e-4)

        # Finally, chromatic with phot
        im4 = obj_chrom.drawImage(bandpass, image=im1.copy(), method='phot', rng=rng)
        print('im4.max,sum = ', im4.array.max(), im4.array.sum())
        print('max diff/flux = ',np.max(np.abs(im1.array-im4.array)/flux))
        np.testing.assert_allclose(im4.array/flux, im1.array/flux, atol=3.e-4)

@timer
def test_chromatic_flux():
    """Test that the total drawn flux is equal to the integral of bandpass * sed over wavelength.
    """
    pixel_scale = 0.5
    stamp_size = 64

    # stars are fundamentally delta-fns with an SED
    star = galsim.Gaussian(fwhm=1e-8) * bulge_SED
    mono_PSF = galsim.Gaussian(half_light_radius=PSF_hlr)
    PSF = galsim.ChromaticAtmosphere(mono_PSF, base_wavelength=500,
                                     zenith_angle=zenith_angle)

    final = galsim.Convolve([star, PSF])
    image = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    image2 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)

    final.drawImage(bandpass, image=image)
    ChromaticConvolve_flux = image.array.sum()

    galsim.ChromaticObject.drawImage(final, bandpass, image=image2)
    ChromaticObject_flux = image2.array.sum()

    # analytic integral...
    analytic_flux = bulge_SED.calculateFlux(bandpass)

    printval(image, image2)
    np.testing.assert_almost_equal(ChromaticObject_flux/analytic_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticObject flux doesn't match " +
                                   "analytic prediction")
    np.testing.assert_almost_equal(ChromaticConvolve_flux/analytic_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "analytic prediction")

    # Also check that the flux is okay and the image fairly consistent when using interpolation
    # for the ChromaticAtmosphere.
    PSF = PSF.interpolate(waves=np.linspace(bandpass.blue_limit, bandpass.red_limit, 30))
    final_int = galsim.Convolve([star, PSF])
    image3 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    final_int.drawImage(bandpass, image=image3)
    int_flux = image3.array.sum()
    # Be *slightly* less stringent in this test given that we did use interpolation.
    printval(image, image3)
    np.testing.assert_almost_equal(
        int_flux/analytic_flux, 1.0, 3,
        err_msg="Drawn ChromaticConvolve flux (interpolated) doesn't match analytic prediction")
    # As an aside, check for appropriate tests of 'integrator' argument.
    assert_raises(ValueError, final_int.drawImage, bandpass, integrator='midp') # minor misspelling
    assert_raises(ValueError, final_int.drawImage, bandpass, integrator=galsim.integ.midptRule)
    check_pickle(PSF)

    # Check option to not use exact SED
    PSF = PSF.deinterpolated
    PSF = PSF * 1.0
    PSF = PSF.interpolate(waves=np.linspace(bandpass.blue_limit, bandpass.red_limit, 30),
                          use_exact_sed=False)
    final_int = galsim.Convolve([star, PSF])
    image3 = galsim.ImageD(stamp_size, stamp_size, scale=pixel_scale)
    final_int.drawImage(bandpass, image=image3)
    int_flux = image3.array.sum()
    # Be *slightly* less stringent in this test given that we did use interpolation.
    printval(image, image3)
    np.testing.assert_almost_equal(
        int_flux/analytic_flux, 1.0, 3,
        err_msg="Drawn ChromaticConvolve flux (interpolated) doesn't match analytic prediction")
    # As an aside, check for appropriate tests of 'integrator' argument.
    assert_raises(ValueError, final_int.drawImage, bandpass, integrator='midp') # minor misspelling
    assert_raises(ValueError, final_int.drawImage, bandpass, integrator=galsim.integ.midptRule)
    check_pickle(PSF)

    # Go back to no interpolation (this will effect the PSFs that are used below).
    PSF = PSF.deinterpolated

    # Try adjusting flux to something else.
    target_flux = 2.63
    bulge_SED2 = bulge_SED.withFlux(target_flux, bandpass)
    star2 = galsim.Gaussian(fwhm=1e-8) * bulge_SED2
    final = galsim.Convolve([star2, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using SED.withFlux()")

    # Use flux_ratio instead.
    flux_ratio = target_flux / analytic_flux
    bulge_SED3 = bulge_SED * flux_ratio
    star3 = galsim.Gaussian(fwhm=1e-8) * bulge_SED3
    final = galsim.Convolve([star3, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using SED * flux_ratio")

    # This should be equivalent.
    bulge_SED3 = flux_ratio * bulge_SED
    star3 = galsim.Gaussian(fwhm=1e-8) * bulge_SED3
    final = galsim.Convolve([star3, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using flux_ratio * SED")

    # Use flux_ratio on the chromatic object instead.
    star4 = star * flux_ratio
    final = galsim.Convolve([star4, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using ChromaticObject * flux_ratio")

    # This should be equivalent.
    star4 = flux_ratio * star
    final = galsim.Convolve([star4, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using flux_ratio * ChromaticObject")

    # As should this.
    star4 = star.withScaledFlux(lambda wave: flux_ratio)
    final = galsim.Convolve([star4, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match " +
                                   "using ChromaticObject.withScaledFlux(flux_ratio)")
    # Can't scale GSObject by function (just SED)
    with assert_raises(TypeError):
        galsim.Gaussian(fwhm=1e-8).withScaledFlux(lambda wave: flux)
    with assert_raises(TypeError):
        galsim.Gaussian(fwhm=1e-8) * (lambda wave: flux)

    # Test ChromaticObject.withFlux
    star5 = star.withFlux(1.0, bandpass)
    final = galsim.Convolve([star5, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum(), 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match "
                                   "using ChromaticObject.withFlux(1.0)")

    # Test withMagnitude.  The zeropoint is the magnitude at which the object produces
    # 1 photon/sec/cm^2 in the given filter.  So drawing with mag == zeropoint should yield a flux
    # of 1.0
    bandpass2 = bandpass.withZeropoint(25.0)
    star6 = star.withMagnitude(25.0, bandpass2)
    final = galsim.Convolve([star6, PSF])
    final.drawImage(bandpass2, image=image)
    np.testing.assert_almost_equal(image.array.sum(), 1.0, 5,
                                   err_msg="Drawn ChromaticConvolve flux doesn't match "
                                   "using ChromaticObject.withMagnitude(0.0)")
    assert_raises(galsim.GalSimError, star.withMagnitude, 25.0, bandpass)

    # Some very simple tests of withFluxDensity.
    star7 = star.withFluxDensity(5.0, 500)
    final = galsim.Convolve([star7, PSF])
    final.evaluateAtWavelength(500).drawImage(image=image)
    np.testing.assert_almost_equal(image.array.sum(), 5.0, 4,
                                   err_msg="Drawn ChromaticConvolve flux density doesn't match "
                                   "using ChromaticObject.withFluxDensity(5.0, 500)")
    np.testing.assert_almost_equal(5.0, final.sed(500), 7,
                                   err_msg="ChromaticObject.sed(500) doesn't match "
                                   "withFluxDensity.")
    from astropy import units
    star8 = star.withFluxDensity(5.0, 5000*units.AA)
    assert star7 == star8
    star9 = star.withFluxDensity(0.5*units.astrophys.photon/(units.s*units.cm**2*units.AA), 500)
    assert star7 == star9

    # Test with non-unit initial flux
    star9 = galsim.Gaussian(fwhm=1e-8, flux=flux_ratio) * bulge_SED
    final = galsim.Convolve([star9, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum()/target_flux, 1.0, 5,
                                   err_msg="obj * SED doesn't respect obj.flux")

    # Check withFlux again
    print('star9 flux = ',star9.calculateFlux(bandpass))
    star10 = star9.withFlux(1.0, bandpass)
    print('star10 flux = ',star10.calculateFlux(bandpass))
    final = galsim.Convolve([star10, PSF])
    final.drawImage(bandpass, image=image)
    np.testing.assert_almost_equal(image.array.sum(), 1.0, 5,
                                   err_msg="withFlux doesn't work when obj has flux.")


@timer
def test_double_ChromaticSum():
    ''' Test logic section of ChromaticConvolve that splits apart ChromaticSums for the case that
    more than one ChromaticSum's are convolved together.
    '''
    a = galsim.Gaussian(fwhm=1.0) * bulge_SED
    b = galsim.Gaussian(fwhm=2.0) * bulge_SED
    c = galsim.Gaussian(fwhm=3.0)
    d = galsim.Gaussian(fwhm=4.0)

    image = galsim.ImageD(16, 16, scale=0.2)
    obj = galsim.Convolve(a+b, c+d)
    obj.drawImage(bandpass, image=image, method='no_pixel')

    check_pickle(obj)

    image_a = galsim.ImageD(16, 16, scale=0.2)
    image_b = galsim.ImageD(16, 16, scale=0.2)
    obj_a = galsim.Convolve(a, c+d)
    obj_b = galsim.Convolve(b, c+d)
    obj_a.drawImage(bandpass, image = image_a, method='no_pixel')
    obj_b.drawImage(bandpass, image = image_b, method='no_pixel')
    printval(image, image_a+image_b)

    np.testing.assert_almost_equal(image.array, (image_a+image_b).array, 5,
                                   err_msg="Convolving two ChromaticSums failed")


@timer
def test_ChromaticSum_nphot():
    """Test n_photons parameter with ChromaticSum photon shooting
    """
    # in response to issue #1156
    rng = galsim.BaseDeviate(12345)
    sed1 = galsim.SED("CWW_E_ext.sed", wave_type='A', flux_type='flambda')
    sed2 = galsim.SED("CWW_Sbc_ext.sed", wave_type='A', flux_type='flambda')
    bandpass = galsim.Bandpass("LSST_r.dat", wave_type="nm")

    flux1 = 990
    flux2 = 20
    flux = flux1 + flux2
    obj1 = (galsim.Gaussian(fwhm=0.6) * sed1).withFlux(flux1, bandpass)
    obj2 = (galsim.Gaussian(fwhm=0.3) * sed2).withFlux(flux2, bandpass)
    obj = obj1 + obj2

    class Counter(galsim.PhotonOp):
        def __init__(self):
            self.nphot = []
            self.meanflux = []

        def applyTo(self, photon_array, local_wcs=None, rng=None):
            self.nphot.append(len(photon_array))
            self.meanflux.append(np.mean(photon_array.flux))

    # Looks okay when n_photons unspecified
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng, poisson_flux=False,
        photon_ops=[counter]
    )
    print("n_photons = 0, poisson_flux=False:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.isclose(np.sum(counter.nphot), flux)
    assert np.isclose(img.array.sum(), flux)
    assert np.isclose(img.added_flux, flux)
    assert np.isclose(img.array.sum(), np.sum(counter.nphot))
    assert np.isclose(counter.nphot[0], flux1)
    assert np.isclose(counter.nphot[1], flux2)
    np.testing.assert_allclose(counter.meanflux, 1.0)

    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng, poisson_flux=True,
        photon_ops=[counter]
    )
    print("n_photons = 0, poisson_flux=True:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.isclose(np.sum(counter.nphot), flux, rtol=0.1)
    assert np.isclose(img.array.sum(), flux, rtol=0.1)
    assert np.isclose(img.added_flux, flux, rtol=0.1)
    assert np.isclose(img.array.sum(), np.sum(counter.nphot))
    assert np.isclose(counter.nphot[1]/counter.nphot[0], flux2/flux1, rtol=0.3)
    np.testing.assert_allclose(counter.meanflux, 1.0)

    # When n_photons is explicit, used to shoot too many photons.  cf. #1156
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng,
        photon_ops=[counter], n_photons=101, poisson_flux=False,
    )
    print("n_photons = 101, poisson_flux=False:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.sum(counter.nphot) == 101
    assert np.isclose(img.array.sum(), flux)
    assert np.isclose(img.added_flux, flux)
    assert np.isclose(counter.nphot[1]/counter.nphot[0], flux2/flux1, rtol=1)
    np.testing.assert_allclose(counter.meanflux, np.mean(counter.meanflux))

    # Repeat with poisson_flux=True
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng,
        photon_ops=[counter], n_photons=101, poisson_flux=True,
    )
    print("n_photons = 101, poisson_flux=True:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.sum(counter.nphot) == 101
    assert np.isclose(img.array.sum(), flux, rtol=0.1)
    assert np.isclose(img.added_flux, flux, rtol=0.1)
    assert np.isclose(counter.nphot[1]/counter.nphot[0], flux2/flux1, rtol=1)
    np.testing.assert_allclose(counter.meanflux, np.mean(counter.meanflux))

    # Do few enough that one component gets no photons.
    rng = galsim.BaseDeviate(1234)
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng,
        photon_ops=[counter], n_photons=11, poisson_flux=False,
    )
    print("n_photons = 11, poisson_flux=False:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.sum(counter.nphot) == 11
    assert len(counter.nphot) == 1  # Not a priori required, but works for this rng.
    assert np.isclose(img.array.sum(), flux)
    assert np.isclose(img.added_flux, flux)

    print('Start #1170 test')
    # If multiple objects at the end has zero flux, it used to cause a ZeroDivisionError.
    # (cf. #1170)
    # Note: they need to all have different seds to avoid ChromaticSum collapsing them.
    obj1 = (galsim.Gaussian(fwhm=0.6) * sed1).withFlux(1.0, bandpass)
    sed3 = galsim.SED("CWW_Im_ext.sed", wave_type='A', flux_type='flambda')
    obj3 = (galsim.Gaussian(fwhm=0.3) * sed3).withFlux(0, bandpass)
    obj = galsim.ChromaticSum([obj1, obj2 * 0, obj3])
    rng = galsim.BaseDeviate(1234)
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng,
        photon_ops=[counter], n_photons=11, poisson_flux=False,
    )
    print("3 objects, two with 0 flux, n_photons=11, poisson_flux=False:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.sum(counter.nphot) == 11
    assert np.isclose(img.array.sum(), 1.0)
    assert np.isclose(img.added_flux, 1.0)

    # Also check if the first two have zero flux.
    obj = galsim.ChromaticSum([obj3, obj1, obj2 * 0])
    counter = Counter()
    img = obj.drawImage(
        bandpass, nx=24, ny=24, scale=0.2, method='phot', rng=rng,
        photon_ops=[counter], n_photons=11, poisson_flux=False,
    )
    print("3 objects, last two with 0 flux, n_photons=11, poisson_flux=False:")
    print("counter.nphot = ",counter.nphot)
    print("counter.meanflux = ",counter.meanflux)
    print("img.array.sum() = ",img.array.sum())
    print("img.added_flux = ",img.added_flux)
    assert np.sum(counter.nphot) == 11
    assert np.isclose(img.array.sum(), 1.0)



@timer
def test_ChromaticConvolution_of_ChromaticConvolution():
    """Check that the __init__ of ChromaticConvolution properly expands arguments that are already
    ChromaticConvolutions.
    """
    a = galsim.Gaussian(fwhm=1.0) * bulge_SED
    b = galsim.Gaussian(fwhm=2.0)
    c = galsim.Gaussian(fwhm=3.0)
    d = galsim.Gaussian(fwhm=4.0)

    e = galsim.Convolve(a, b)
    f = galsim.Convolve(c, d)
    g = galsim.Convolve(e, f)
    if any(isinstance(h, galsim.ChromaticConvolution) for h in g.obj_list):
        raise AssertionError("ChromaticConvolution did not expand ChromaticConvolution argument")

    assert_raises(TypeError, galsim.ChromaticConvolution)
    assert_raises(TypeError, galsim.ChromaticConvolution, bulge_SED)
    assert_raises(TypeError, galsim.ChromaticConvolution, [a,b], invalid=True)
    assert_raises(NotImplementedError, galsim.ChromaticConvolution, [a,b], real_space=True)


@timer
def test_ChromaticAutoConvolution():
    a = galsim.Gaussian(fwhm=1.0)
    b = galsim.Gaussian(fwhm=2.0) * bulge_SED
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    c = galsim.Convolve(a, a, b)
    c.drawImage(bandpass, image=im1, method='no_pixel')
    d = galsim.Convolve(galsim.AutoConvolve(a), b)
    d.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(im1.array, im2.array, 5,
                                         "ChromaticAutoConvolution(a) not equal to "
                                         "ChromaticConvolution(a,a)")

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (d * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="ChromaticAutoConvolution * 2 resulted in wrong flux.")


@timer
def test_ChromaticAutoCorrelation():
    a = galsim.Gaussian(fwhm=1.0)
    b = galsim.Gaussian(fwhm=2.0) * bulge_SED
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    c = galsim.Convolve(a, a.rotate(180.0 * galsim.degrees), b)
    c.drawImage(bandpass, image=im1, method='no_pixel')
    d = galsim.Convolve(galsim.AutoCorrelate(a), b)
    d.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(im1.array, im2.array, 5,
                                         "ChromaticAutoCorrelate(a) not equal to "
                                         "ChromaticConvolution(a,a.rotate(180.0*galsim.degrees)")

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (d * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2/(2.*flux), 1.0, 5,
        err_msg="ChromaticAutoCorrelation * 2 resulted in wrong flux.")


@timer
def test_ChromaticObject_expand():
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    a = galsim.Gaussian(fwhm=1.0).expand(1.1) * bulge_SED
    b = (galsim.Gaussian(fwhm=1.0) * bulge_SED).expand(1.1)

    a.drawImage(bandpass, image=im1, method='no_pixel')
    b.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(im1.array, im2.array, 5,
                                         "ChromaticObject.expand not equal to Chromatic.expand")

    check_pickle(a)
    check_pickle(b)

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (b * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="expanded ChromaticObject * 2 resulted in wrong flux.")

    # Test the wavelength-dependent expansion option
    pixel_scale = 0.3
    sigma = 1.7
    gal = galsim.Gaussian(sigma=sigma, flux=132)
    cgal = galsim.ChromaticObject(gal) * galsim.SED('1', 'nm', 'fphotons')
    im0 = gal.drawImage(scale=pixel_scale, dtype=float, method='no_pixel')
    # Initial mxx,myy should be sigma**2 in units of arcsec^2, so (sigma/pixel_scale)**2 in pix.
    mom = galsim.utilities.unweighted_moments(im0)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, 1.0, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, 1.0, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # First a very simple case with no actual wavelength dependence, but using the
    # functional syntax.
    gal1 = cgal.expand(lambda w: 1.2)
    # Use a simple bandpass so we can do the integral below analytically
    bp = galsim.Bandpass(lambda w: 1. - 0.12*(w-600)**2/100**2, 'nm', 500, 700)
    im1 = gal1.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im1)
    print('simple growth = ',mom['Mxx']/(sigma/pixel_scale)**2, 1.2**2)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, 1.2**2, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, 1.2**2, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # Use an expansion that varies quadratically within the range 500-700
    expansion = lambda w: 1.2 + 0.11*(w-600)/100 - 0.36*(w-600)**2/100**2
    gal2 = cgal.expand(expansion)
    im2 = gal2.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    # The test here relies on linearity to swap the order of the integrals.
    # The moments test is essentially int(dxdy x^2 int(dlambda I(x,y,lam) F(lam)))
    # But we can swap the order of the integrals to do int(dlambda int(dxdy x^2 I(x,y,lam) F(lam)))
    # and do the integral analytically.
    mom = galsim.utilities.unweighted_moments(im2)
    # Now the analytic integral:
    # mxx = int(dlam (sigma * expansion)**2 * expansion**2 * F(lam))
    #         / int(dlam expansion**2 * F(lam))
    # Note that two of the factors of expansion are the flux scaling.  The other two, only in the
    # numerator, are the size scaling.
    #     = sigma**2 * int(dt (1.2 + 0.11 t - 0.36 t**2)**4 (1 - 0.12 t**2), t=-1..1)
    #            / int(dt (1.2 + 0.11 t - 0.36 t**2)**2 (1 - 0.12 t**2), t=-1..1)
    #     = sigma**2 * 1.243224162 (according to Maple)
    growth_factor = galsim.integ.int1d(lambda w: expansion(w)**4 * bp(w),500,700)
    growth_factor /= galsim.integ.int1d(lambda w: expansion(w)**2 * bp(w),500,700)
    print('growth factor = ',mom['Mxx']/(sigma/pixel_scale)**2, growth_factor)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # Repeat using transform rather than expand
    gal3 = cgal.transform(expansion, 0.0, 0.0, expansion)
    im3 = gal3.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im3)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # Repeat using magnify
    gal4 = cgal.magnify(lambda w: expansion(w)**2)
    im4 = gal4.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # Repeat using lens
    gal5 = cgal.lens(0., 0., lambda w: expansion(w)**2)
    im5 = gal5.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im5)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # Dilate isn't quite the same, since it doesn't have the extra flux factor:
    gal5 = cgal.dilate(expansion)
    im5 = gal5.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im5)
    dilate_growth_factor = galsim.integ.int1d(lambda w: expansion(w)**2 * bp(w),500,700)
    dilate_growth_factor /= galsim.integ.int1d(lambda w: bp(w),500,700)
    print('dilate_growth factor = ',mom['Mxx']/(sigma/pixel_scale)**2, dilate_growth_factor)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, dilate_growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, dilate_growth_factor, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)


@timer
def test_ChromaticObject_rotate():
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    a = (galsim.Gaussian(fwhm=1.0)
         .shear(eta=0.1, beta=0 * galsim.degrees)
         .rotate(1.1 * galsim.radians)) * bulge_SED
    b = (((galsim.Gaussian(fwhm=1.0)
           .shear(eta=0.1, beta=0 * galsim.degrees)) * bulge_SED)
           .rotate(1.1 * galsim.radians))

    a.drawImage(bandpass, image=im1, method='no_pixel')
    b.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(
            im1.array, im2.array, 5,
            "ChromaticObject.rotate not equal to ChromaticTransformation.rotate")

    check_pickle(a)
    check_pickle(b)

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (b * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="rotated ChromaticObject * 2 resulted in wrong flux.")

    # Test the wavelength-dependent rotation option
    pixel_scale = 0.3
    sigma = 1.7
    gal = galsim.Gaussian(sigma=sigma, flux=132).shear(e1=0.3)
    cgal = galsim.ChromaticObject(gal) * galsim.SED('1', 'nm', 'fphotons')
    im0 = gal.drawImage(scale=pixel_scale, dtype=float, method='no_pixel')
    # Initial distortion should be (e1,0).
    mom = galsim.utilities.unweighted_moments(im0)
    print('e1 = ',(mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']))
    print('e2 = ',(2*mom['Mxy'])/(mom['Mxx']+mom['Myy']))
    print('(mxx+myy)/sigma^2 = ',(mom['Mxx']+mom['Myy'])/(sigma/pixel_scale)**2)
    print('2/(mxx+myy)/sigma^2 = ',2./((mom['Mxx']+mom['Myy'])/(sigma/pixel_scale)**2))
    print('sqrt(1-0.3**2) = ',np.sqrt(1.-0.3**2))
    fact = np.sqrt(1.-0.3**2)
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), 0.3, decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), 0.0, decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, (1+0.3)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, (1-0.3)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0, decimal=4)

    # First a very simple case with no actual wavelength dependence, but using the
    # functional syntax.
    gal1 = cgal.rotate(lambda w: 0.4 * galsim.radians)
    bp = galsim.Bandpass(lambda w: 1. - 0.12*(w-600)**2/100**2, 'nm', 500, 700)
    im1 = gal1.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im1)
    print('simple angle = ',(np.arctan2(2.*mom['Mxy'],mom['Mxx']-mom['Myy'])/2.), 0.4)
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), 0.3*np.cos(2.*0.4), decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), 0.3*np.sin(2.*0.4), decimal=4)
    rot_e1 = 0.3*np.cos(2.*0.4)
    rot_e2 = 0.3*np.sin(2.*0.4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, (1+rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, (1-rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, rot_e2/fact, decimal=4)

    # Use a rotation that varies quadratically within the range 500-700
    rotation = lambda w: (0.4 + 0.11*(w-600)/100 - 0.36*(w-600)**2/100**2) * galsim.radians
    gal2 = cgal.rotate(rotation)
    im2 = gal2.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im2)
    # The analytic integrals:
    # mxx = int(dw (1 + 0.3 cos(2theta(w))) * F(w))  / int(dw F(w))
    # myy = int(dw (1 - 0.3 cos(2theta(w))) * F(w))  / int(dw F(w))
    # mxy = int(dw 0.3 sin(2theta(w)) * F(w))  / int(dw F(w))
    rot_mxx = galsim.integ.int1d(lambda w: (1.+0.3*np.cos(2*rotation(w))) * bp(w),500,700)
    rot_mxx /= galsim.integ.int1d(lambda w: bp(w),500,700)
    rot_myy = galsim.integ.int1d(lambda w: (1.-0.3*np.cos(2*rotation(w))) * bp(w),500,700)
    rot_myy /= galsim.integ.int1d(lambda w: bp(w),500,700)
    rot_mxy = galsim.integ.int1d(lambda w: (0.3*np.sin(2*rotation(w))) * bp(w),500,700)
    rot_mxy /= galsim.integ.int1d(lambda w: bp(w),500,700)
    rot_e1 = (rot_mxx-rot_myy)/(rot_mxx+rot_myy)
    rot_e2 = (2*rot_mxy)/(rot_mxx+rot_myy)
    print('rot e1 = ',(mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), rot_e1)
    print('rot e2 = ',(2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), rot_e2)
    print('rot e = ',np.sqrt(rot_e1**2+rot_e2**2))
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), rot_e1, decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), rot_e2, decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, (1+rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, (1-rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, rot_e2/fact, decimal=4)

    # Repeat using transform rather than rotate
    gal3 = cgal.transform(
                lambda w: np.cos(rotation(w)),
                lambda w: -np.sin(rotation(w)),
                lambda w: np.sin(rotation(w)),
                lambda w: np.cos(rotation(w)) )
    im3 = gal3.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im3)
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), rot_e1, decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), rot_e2, decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, (1+rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, (1-rot_e1)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, rot_e2/fact, decimal=4)


@timer
def test_ChromaticObject_shear():
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    shear = galsim.Shear(g1=0.1, g2=0.1)
    a = galsim.Gaussian(fwhm=1.0).shear(shear) * bulge_SED
    b = (galsim.Gaussian(fwhm=1.0) * bulge_SED).shear(shear)

    a.drawImage(bandpass, image=im1, method='no_pixel')
    b.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(
            im1.array, im2.array, 5,
            "ChromaticObject.shear not equal to ChromaticTransformation.shear")

    check_pickle(a)
    check_pickle(b)

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (b * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="rotated ChromaticObject * 2 resulted in wrong flux.")

    # Test the wavelength-dependent shear option
    pixel_scale = 0.3
    sigma = 1.7
    gal = galsim.Gaussian(sigma=sigma, flux=132)
    cgal = galsim.ChromaticObject(gal) * galsim.SED('1', 'nm', 'fphotons')
    im0 = gal.drawImage(scale=pixel_scale, dtype=float, method='no_pixel')
    # Initial distortion should be (0,0).
    mom = galsim.utilities.unweighted_moments(im0)
    print('mxx+myy = ',mom['Mxx']+mom['Myy'])
    print('e1 = ',(mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']))
    print('e2 = ',(2*mom['Mxy'])/(mom['Mxx']+mom['Myy']))
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), 0.0, decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), 0.0, decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, 1.0, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, 1.0, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0.0, decimal=4)

    # First a very simple case with no actual wavelength dependence, but using the
    # functional syntax.
    gal1 = cgal.shear(lambda w: galsim.Shear(e1=0.23, e2=0.13))
    bp = galsim.Bandpass(lambda w: 1. - 0.12*(w-600)**2/100**2, 'nm', 500, 700)
    im1 = gal1.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im1)
    print('mxx+myy = ',mom['Mxx']+mom['Myy'])
    print('simple e1 = ',(mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']))
    print('simple e2 = ',(2.*mom['Mxy'])/(mom['Mxx']+mom['Myy']))
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), 0.23, decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), 0.13, decimal=4)
    fact = np.sqrt(1. - 0.23**2 - 0.13**2)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, (1+0.23)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, (1-0.23)/fact, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, 0.13/fact, decimal=4)

    # Use a shear that varies quadratically within the range 500-700
    shear = lambda w: galsim.Shear(e1=0.23 + 0.11*(w-600)/100 - 0.36*(w-600)**2/100**2,
                                   e2=0.13 + 0.19*(w-600)/100 - 0.09*(w-600)**2/100**2)
    gal2 = cgal.shear(shear)
    im2 = gal2.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im2)
    # The analytic integrals:
    # mxx = int(dw (1 + e1(w))/sqrt(1-e(w)^2) * F(w))  / int(dw F(w))
    # myy = int(dw (1 - e1(w))/sqrt(1-e(w)^2) * F(w))  / int(dw F(w))
    # mxy = int(dw e2(w)/sqrt(1-e(w)^2) * F(w))  / int(dw F(w))
    sh_mxx = galsim.integ.int1d(lambda w: (1.+shear(w).e1)/np.sqrt(1.-shear(w).e**2)*bp(w),500,700)
    sh_mxx /= galsim.integ.int1d(lambda w: bp(w),500,700)
    sh_myy = galsim.integ.int1d(lambda w: (1.-shear(w).e1)/np.sqrt(1.-shear(w).e**2)*bp(w),500,700)
    sh_myy /= galsim.integ.int1d(lambda w: bp(w),500,700)
    sh_mxy = galsim.integ.int1d(lambda w: shear(w).e2/np.sqrt(1.-shear(w).e**2)*bp(w),500,700)
    sh_mxy /= galsim.integ.int1d(lambda w: bp(w),500,700)
    print('mxx+myy = ',mom['Mxx']+mom['Myy'])
    print('shear e1 = ',(mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), (sh_mxx-sh_myy)/(sh_mxx+sh_myy))
    print('shear e2 = ',(2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), (2*sh_mxy)/(sh_mxx+sh_myy))
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), (sh_mxx-sh_myy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), (2*sh_mxy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, sh_mxx, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, sh_myy, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, sh_mxy, decimal=4)

    # Repeat using transform rather than shear
    gal3 = cgal.transform(
                lambda w: (1.+shear(w).g1)/np.sqrt(1.-shear(w).g**2),
                lambda w: shear(w).g2/np.sqrt(1.-shear(w).g**2),
                lambda w: shear(w).g2/np.sqrt(1.-shear(w).g**2),
                lambda w: (1.-shear(w).g1)/np.sqrt(1.-shear(w).g**2) )
    im3 = gal3.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im3)
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), (sh_mxx-sh_myy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), (2*sh_mxy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, sh_mxx, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, sh_myy, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, sh_mxy, decimal=4)

    # Repeat using lens
    gal4 = cgal.lens(lambda w: shear(w).g1, lambda w: shear(w).g2, 1.)
    im4 = gal4.drawImage(bp, scale=pixel_scale, dtype=float, method='no_pixel')
    mom = galsim.utilities.unweighted_moments(im4)
    np.testing.assert_almost_equal((mom['Mxx']-mom['Myy'])/(mom['Mxx']+mom['Myy']), (sh_mxx-sh_myy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal((2*mom['Mxy'])/(mom['Mxx']+mom['Myy']), (2*sh_mxy)/(sh_mxx+sh_myy), decimal=4)
    np.testing.assert_almost_equal(mom['Mxx'] / (sigma/pixel_scale)**2, sh_mxx, decimal=4)
    np.testing.assert_almost_equal(mom['Myy'] / (sigma/pixel_scale)**2, sh_myy, decimal=4)
    np.testing.assert_almost_equal(mom['Mxy'] / (sigma/pixel_scale)**2, sh_mxy, decimal=4)

    assert_raises(TypeError, cgal.shear, 0.1, 0.3)
    assert_raises(TypeError, cgal.shear, 0.1)
    assert_raises(TypeError, cgal.shear, shear, g1=0.1, g2=0.2)
    assert_raises(TypeError, cgal.shear, shear=shear, g1=0.1, g2=0.2)


@timer
def test_ChromaticObject_shift():
    im1 = galsim.ImageD(32, 32, scale=0.2)
    im2 = galsim.ImageD(32, 32, scale=0.2)
    shift = (0.1, 0.3)
    a = galsim.Gaussian(fwhm=1.0).shift(shift) * bulge_SED
    b = (galsim.Gaussian(fwhm=1.0) * bulge_SED).shift(shift)

    a.drawImage(bandpass, image=im1, method='no_pixel')
    b.drawImage(bandpass, image=im2, method='no_pixel')
    printval(im1, im2)
    np.testing.assert_array_almost_equal(
            im1.array, im2.array, 5,
            "ChromaticObject.shift not equal to ChromaticTransformation.shift")

    check_pickle(a)
    check_pickle(b)

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (b * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="rotated ChromaticObject * 2 resulted in wrong flux.")

    cgal = galsim.Gaussian(fwhm=1.0) * bulge_SED
    assert_raises(TypeError, cgal.shift)
    assert_raises(TypeError, cgal.shift, 0.1)
    assert_raises(TypeError, cgal.shift, shift, 0.1, 0.2)
    assert_raises(TypeError, cgal.shift, shift, dx=0.1, dy=0.2)
    assert_raises(TypeError, cgal.shift, shift=shift)

@timer
def test_ChromaticObject_compound_affine_transformation():
    """ Check that making a (separable) object chromatic before a bunch of transformations is
    equivalent to making it chromatic after a bunch of transformations.
    """
    shear = galsim.Shear(eta=1.0, beta=0.3*galsim.radians)
    scale = 1.1
    theta = 0.1 * galsim.radians
    shift = (0.1, 0.3)
    sed = galsim.SED('wave**0.3', 'nm', 'fphotons')
    bandpass = galsim.Bandpass('1', 'nm', blue_limit=400, red_limit=550)

    a = galsim.Gaussian(fwhm=1.0)
    a = a.shear(shear).shift(shift).rotate(theta).dilate(scale)
    a = a._shear(shear)._shift(*shift).rotate(theta).expand(scale)
    a = a._lens(g1=0.1, g2=0.1, mu=1.1).shift(shift).rotate(theta).magnify(scale)
    a = a * sed

    b = galsim.Gaussian(fwhm=1.0) * sed
    b = b.shear(shear).shift(shift).rotate(theta).dilate(scale)
    b = b._shear(shear)._shift(*shift).rotate(theta).expand(scale)
    b = b._lens(g1=0.1, g2=0.1, mu=1.1).shift(shift).rotate(theta).magnify(scale)

    # Include a few gratuitous combinations of functional and static values.
    pshift = galsim.PositionD(*shift)
    c = galsim.Gaussian(fwhm=1.0) * sed
    c = galsim.Transform(c, flux_ratio=lambda w:1, jac=shear.getMatrix())
    c = c.shift(lambda w: shift).rotate(theta).dilate(lambda w: scale)
    c = c.shear(shear).shift(lambda w: pshift).rotate(theta).expand(scale)
    c = c.lens(g1=lambda w:0.1, g2=0.1, mu=lambda w:1.1).shift(shift).rotate(theta).magnify(scale)

    d = galsim.Gaussian(fwhm=1.0)
    d = galsim.Transform(d, flux_ratio=sed, jac=shear.getMatrix())
    d = d.shift(pshift).rotate(lambda w: theta).dilate(scale)
    d = d.shear(lambda w: shear).shift(shift).rotate(theta)
    d = d.transform(scale, lambda w:0, lambda w:0, scale)
    d = d.lens(g1=0.1, g2=lambda w:0.1, mu=1.1).shift(shift).rotate(theta).magnify(scale)

    im1 = galsim.ImageD(32, 32, scale=0.2)
    im1 = a.drawImage(bandpass, image=im1.copy(), method='no_pixel')
    im2 = b.drawImage(bandpass, image=im1.copy(), method='no_pixel')
    im3 = c.drawImage(bandpass, image=im1.copy(), method='no_pixel')
    im4 = d.drawImage(bandpass, image=im1.copy(), method='no_pixel')
    printval(im1, im2)
    for test_im in (im2, im3, im4):
        np.testing.assert_array_almost_equal(test_im.array, im1.array, 5,
                                            "ChromaticObject affine transformation not equal to "
                                            "GSObject affine transformation")

    # Check flux scaling
    flux = im2.array.sum()
    im2 = (b * 2.).drawImage(bandpass, image=im2, method='no_pixel')
    flux2 = im2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux, 5,
        err_msg="transformed ChromaticObject * 2 resulted in wrong flux.")

    # Just check that the cache resizing routines are what the docs say they are.
    galsim.ChromaticObject.resize_multiplier_cache(100)
    galsim.ChromaticConvolution.resize_effective_prof_cache(100)

    # Check some branches in repr that we wouldn't hit otherwise.
    repr(a); repr(b); repr(c); repr(d)


@timer
def test_analytic_integrator():
    """Test that the analytic (i.e., not sampled) versions of SEDs and Bandpasses produce the
    same results as the sampled versions.
    """
    # Make sure to use something non-separable so the ImageIntegrators actually get triggered.
    psf = galsim.ChromaticObject(galsim.Moffat(fwhm=1.0, beta=2.7)).dilate(lambda w:(w/500)**1.1)

    # pure analytic
    band1 = galsim.Bandpass('1', 'nm', blue_limit=500, red_limit=750)
    sed1 = galsim.SED('wave**1.1', wave_type='nm', flux_type='fphotons').withFluxDensity(1.0, 500)
    gal1 = galsim.Gaussian(fwhm=1.0) * sed1
    final1 = galsim.Convolve(gal1, psf)
    image1 = galsim.ImageD(32, 32, scale=0.2)
    kimage1 = galsim.ImageCD(32, 32, scale=0.2)
    kimage1.setCenter(0,0)  # Necessary when adding to image
    assert len(band1.wave_list) == 0
    assert len(sed1.wave_list) == 0
    # Gratuitously use add_to_image=True so we test that branch.  Initial image is 0.
    final1.drawImage(band1, image=image1, add_to_image=True)
    final1.drawKImage(band1, image=kimage1, add_to_image=True)

    # try making the SED sampled
    band2 = band1
    N = 250 # default N for ContinuousIntegrator
    h = (band2.red_limit*1.0 - band2.blue_limit)/N
    x = [band2.blue_limit + h * i for i in range(N+1)]
    # make a sampled SED
    sed2 = galsim.SED(galsim.LookupTable(x, sed1(x), interpolant='linear'),
                      wave_type='nm', flux_type='fphotons')
    gal2 = galsim.Gaussian(fwhm=1.0) * sed2
    final2 = galsim.Convolve(gal2, psf)
    image2 = galsim.ImageD(32, 32, scale=0.2)
    kimage2 = galsim.ImageCD(32, 32, scale=0.2)
    assert len(band2.wave_list) == 0
    assert len(sed2.wave_list) != 0
    final2.drawImage(band1, image=image2)
    final2.drawKImage(band1, image=kimage2)

    # try making the Bandpass sampled
    sed3 = sed1
    band3 = galsim.Bandpass(galsim.LookupTable(x, band1(x), interpolant='linear'), 'nm')
    gal3 = galsim.Gaussian(fwhm=1.0) * sed3
    final3 = galsim.Convolve(gal3, psf)
    image3 = galsim.ImageD(32, 32, scale=0.2)
    kimage3 = galsim.ImageCD(32, 32, scale=0.2)
    assert len(band3.wave_list) != 0
    assert len(sed3.wave_list) == 0
    final3.drawImage(band3, image=image3)
    final3.drawKImage(band3, image=kimage3)

    printval(image1, image2)
    np.testing.assert_array_almost_equal(image1.array, image2.array, 5,
                                         "Analytic integrator doesn't match sample integrator")
    printval(image1, image3)
    np.testing.assert_array_almost_equal(image1.array, image3.array, 5,
                                         "Analytic integrator doesn't match sample integrator")

    np.testing.assert_array_almost_equal(kimage1.array, kimage2.array, 5,
                                         "Analytic integrator doesn't match sample integrator")
    np.testing.assert_array_almost_equal(kimage1.array, kimage3.array, 5,
                                         "Analytic integrator doesn't match sample integrator")

    # Test that attempting to use SampleIntegrator with analytic sed, bandpass raises an Error:
    with assert_raises(ValueError):
        final1.drawImage(band1,
                         integrator=galsim.integ.SampleIntegrator(rule=galsim.integ.trapzRule))


@timer
def test_gsparams():
    """Check that gsparams actually gets processed by ChromaticObjects.
    """
    # Setting maximum_fft_size this low causes a warning to be emitted for GSObjects, so
    # make sure it does for ChromaticObjects too, thereby assuring that gsparams is really
    # getting properly forwarded through the internals of ChromaticObjects.
    gsparams = galsim.GSParams(maximum_fft_size=16)
    gal = galsim.Gaussian(fwhm=1, gsparams=gsparams) * bulge_SED
    with assert_warns(galsim.GalSimFFTSizeWarning):
        gal.drawImage(bandpass)
    assert (galsim.Gaussian(fwhm=1) * bulge_SED) != gal
    assert (galsim.Gaussian(fwhm=1) * bulge_SED).withGSParams(gsparams) == gal
    assert (galsim.Gaussian(fwhm=1) * bulge_SED).withGSParams(maximum_fft_size=16) == gal

    # Should raise an exception if raise_fft_size_error is True
    galsim.errors.raise_fft_size_error = True
    with assert_raises(galsim.GalSimFFTSizeError):
        gal.drawImage(bandpass)
    galsim.errors.raise_fft_size_error = False

    # Repeat, putting the gsparams argument in after the ChromaticObject constructor.
    gal = galsim.Gaussian(fwhm=1) * bulge_SED
    psf = galsim.Gaussian(sigma=0.4)
    final = galsim.Convolve([gal, psf], gsparams=gsparams)
    with assert_warns(galsim.GalSimFFTSizeWarning):
        final.drawImage(bandpass)

    galsim.errors.raise_fft_size_error = True
    with assert_raises(galsim.GalSimFFTSizeError):
        final.drawImage(bandpass)
    galsim.errors.raise_fft_size_error = False

    # Use a restrictive one this time, so we test the "most restrictive gsparams" feature
    gsp2 = galsim.GSParams(folding_threshold=1.e-4, maxk_threshold=1.e-4, maximum_fft_size=1.e4)
    final = galsim.Convolve([gal, psf])
    final1 = galsim.Convolve([gal, psf], gsparams=gsp2)
    final2 = galsim.Convolve([gal, psf]).withGSParams(gsp2)
    final3 = galsim.Convolve([gal.withGSParams(gsp2), psf])
    final4 = galsim.Convolve([gal, psf.withGSParams(gsp2)])
    assert final1 != final
    assert final2 == final1
    assert final3 == final1
    assert final4 == final1
    assert final1.gsparams == gsp2
    assert final1.obj_list[0].gsparams == gsp2
    assert final1.obj_list[1].gsparams == gsp2

    final5 = galsim.Convolve(gal, psf, gsparams=gsp2, propagate_gsparams=False)
    assert final5 != final1
    assert final5.gsparams == gsp2
    assert final5.obj_list[0].gsparams == galsim.GSParams()
    assert final5.obj_list[1].gsparams == galsim.GSParams()

    final6 = final5.withGSParams(gsparams)
    assert final6 != final5
    assert final6.gsparams == gsparams
    assert final6.obj_list[0].gsparams == galsim.GSParams()
    assert final6.obj_list[1].gsparams == galsim.GSParams()

    check_pickle(final1)

    # Repeat similar tests for ChromaticSum
    bulge = galsim.Gaussian(half_light_radius=1) * bulge_SED
    disk = galsim.Exponential(half_light_radius=2) * disk_SED
    final = galsim.Add([bulge, disk])
    final1 = galsim.Add([bulge, disk], gsparams=gsp2)
    final2 = galsim.Add([bulge, disk]).withGSParams(gsp2)
    final3 = galsim.Add([bulge.withGSParams(gsp2), disk])
    final4 = galsim.Add([bulge, disk.withGSParams(gsp2)])
    assert final1 != final
    assert final2 == final1
    assert final3 == final1
    assert final4 == final1
    assert final1.gsparams == gsp2
    assert final1.obj_list[0].gsparams == gsp2
    assert final1.obj_list[1].gsparams == gsp2

    final5 = galsim.Add(bulge, disk, gsparams=gsp2, propagate_gsparams=False)
    assert final5 != final1
    assert final5.gsparams == gsp2
    assert final5.obj_list[0].gsparams == galsim.GSParams()
    assert final5.obj_list[1].gsparams == galsim.GSParams()

    final6 = final5.withGSParams(gsparams)
    assert final6 != final5
    assert final6.gsparams == gsparams
    assert final6.obj_list[0].gsparams == galsim.GSParams()
    assert final6.obj_list[1].gsparams == galsim.GSParams()


@timer
def test_separable_ChromaticSum():
    """ Test ChromaticSum separable profile grouping.
    """
    psf = galsim.Gaussian(fwhm=1)
    mono_gal1 = galsim.Gaussian(fwhm=1)
    mono_gal2 = galsim.Gaussian(fwhm=1.1)
    mono_gal3 = galsim.Gaussian(fwhm=1.2)

    # 1) check 2 summands with same SED.  (Used to be considered separable.)
    cgal1 = mono_gal1 * bulge_SED + mono_gal2 * bulge_SED
    img1 = galsim.ImageD(32, 32, scale=0.2)
    kimg1 = galsim.ImageCD(32, 32, scale=0.5)
    assert not cgal1.separable

    # check that drawing the profile works as expected
    final1 = galsim.Convolve(cgal1, psf)
    final1.drawImage(bandpass, image=img1)
    final1.drawKImage(bandpass, image=kimg1)

    img1b = img1.copy()
    kimg1b = kimg1.copy()
    component1 = galsim.Convolve(mono_gal1*bulge_SED, psf)
    component1.drawImage(bandpass, image=img1b)
    component1.drawKImage(bandpass, image=kimg1b)

    component2 = galsim.Convolve(mono_gal2*bulge_SED, psf)
    component2.drawImage(bandpass, image=img1b, add_to_image=True)
    component2.drawKImage(bandpass, image=kimg1b, add_to_image=True)

    np.testing.assert_array_almost_equal(img1.array, img1b.array, 5,
                                         "separable ChromaticSum not correctly drawn")
    np.testing.assert_array_almost_equal(kimg1.array, kimg1b.array, 5,
                                         "separable ChromaticSum not correctly k-drawn")
    check_pickle(final1)

    # 2) Check flux scaling
    img2 = img1.copy()
    flux1 = img1.array.sum()
    (final1 * 2).drawImage(bandpass, image=img2)
    flux2 = img2.array.sum()
    np.testing.assert_array_almost_equal(
        flux2, 2.*flux1, 5,
        err_msg="ChromaticConvolution containing separable ChromaticSum * 2 resulted in wrong flux.")

    final2 = galsim.Convolve(cgal1 * 2, psf)
    img2b = img1.copy()
    final2.drawImage(bandpass, image=img2b)
    flux2b = img2b.array.sum()
    np.testing.assert_array_almost_equal(
        flux2b, 2.*flux1, 5,
        err_msg="separable ChromaticSum * 2 resulted in wrong flux.")

    check_pickle(final2)

    # 3) check that 3 summands, 2 with the same SED, 1 with a different SED, make an
    # inseparable sum.
    cgal3 = galsim.Add(mono_gal1 * bulge_SED, mono_gal2 * bulge_SED, mono_gal3 * disk_SED)
    assert not cgal3.separable
    assert len(cgal3.obj_list) == 3
    assert all([isinstance(obj, galsim.ChromaticTransformation) for obj in cgal3.obj_list])

    # check that drawing the profile works as expected
    final3 = galsim.Convolve(cgal3, psf)
    img3 = img1.copy()
    final3.drawImage(bandpass, image=img3)

    check_pickle(final3)


    component3 = galsim.Convolve(mono_gal3*disk_SED, psf)
    img3b = img1.copy()
    component1.drawImage(bandpass, image=img3b)
    component2.drawImage(bandpass, image=img3b, add_to_image=True)
    component3.drawImage(bandpass, image=img3b, add_to_image=True)

    np.testing.assert_array_almost_equal(img3.array, img3b.array, 5,
                                         "inseparable ChromaticSum not correctly drawn")

    check_pickle(component3)

    # 4) What about if we add them first before applying the SED?
    cgal4 = (mono_gal1 + mono_gal2) * bulge_SED
    img4 = img1.copy()
    assert cgal4.separable

    final4 = galsim.Convolve(cgal4, psf)
    final4.drawImage(bandpass, image=img4)
    np.testing.assert_array_almost_equal(img1.array, img4.array, 5,
                                         "separable ChromaticSum not correctly drawn")

    assert_raises(TypeError, galsim.ChromaticSum,
                  [mono_gal1 * bulge_SED, mono_gal2 * bulge_SED], invalid=3)
    assert_raises(TypeError, galsim.ChromaticSum)
    assert_raises(TypeError, galsim.ChromaticSum, bulge_SED)
    with assert_raises(galsim.GalSimIncompatibleValuesError):
        sum = mono_gal1 * bulge_SED + mono_gal2


@timer
def test_centroid():
    """Test the ChromaticObject.calculateCentroid function."""
    sed = galsim.SED('wave', wave_type='nm', flux_type='fphotons')
    bp = galsim.Bandpass('wave', 'nm', blue_limit=0, red_limit=1)
    shift_fn = lambda w: (w, 0)
    gal = sed * galsim.Gaussian(fwhm=1)
    gal = gal.shift(shift_fn)
    # The sed and bandpass each contribute a factor of wavelength to the flux integrand of the
    # galaxy.  The shift function contributes an additional factor of wavelength to the x-centroid
    # integrand.  The end result is that the x-centroid should be:
    # int(w^3, 0, 1) / int(w^2, 0, 1) = (1/4)/(1/3) = 3/4.
    cen = gal.calculateCentroid(bp)
    np.testing.assert_almost_equal(cen.x, 0.75, 5, "ChromaticObject.calculateCentroid() failed")
    np.testing.assert_almost_equal(cen.y, 0.0, 5, "ChromaticObject.calculateCentroid() failed")

    # Now check the centroid sampling integrator...
    # (Writing to the dict is a quick and dirty way to overwrite an attribute that is a property.)
    gal.__dict__['wave_list'] = np.linspace(0.0, 1.0, 500)
    cen = gal.calculateCentroid(bp)
    np.testing.assert_almost_equal(cen.x, 0.75, 5, "ChromaticObject.calculateCentroid() failed")
    np.testing.assert_almost_equal(cen.y, 0.0, 5, "ChromaticObject.calculateCentroid() failed")


@timer
def test_interpolated_ChromaticObject():
    """Test the ChromaticObject interpolation functionality."""
    # Define a ChromaticObject subclass for which interpolation could be helpful, but calculations
    # can be done quickly:
    # It's a Gaussian with a wavelength-dependent sigma [sigma = sigma_0 (wave / 500 nm)]
    # and shear g1 [g1 = 0.1 ((wave)/(500 nm) - 1)]
    class ChromaticGaussian(galsim.ChromaticObject):
        def __init__(self, sigma):
            # First, take the basic info.
            self.sigma = sigma
            self.separable = False
            self.interpolated = False
            self.deinterpolated = self

        @property
        def sed(self):
            return galsim.SED(1, 'nm', '1')

        @property
        def wave_list(self):
            return np.array([], dtype=float)

        def evaluateAtWavelength(self, wave):
            this_sigma = self.sigma * (wave / 500.)
            this_shear = 0.1 * ((wave/500.)-1.)
            ret = galsim.Gaussian(sigma = this_sigma)
            ret = ret.shear(g1 = this_shear)
            return ret

        @property
        def gsparams(self): return galsim.GSParams.default

        def withGSParams(self, gsp): return self

        def __eq__(self, other):
            return (isinstance(other, ChromaticGaussian) and
                    self.sigma == other.sigma)

        def __hash__(self):
            return hash(("ChromaticGaussian", self.sigma))

        def __repr__(self):
            return 'galsim.ChromaticGaussian(%r)'%self.sigma

    # For this test, we're going to use the ChromaticGaussian defined above.  This class is simple
    # enough that evaluation of both the interpolated and exact versions is very fast, so it won't
    # slow down the tests too much to do both ways.  Note that for initial tests (fair comparison
    # between exact and interpolated version) we're going to use a lot of images for the
    # interpolation, n_interp=100.  For more slowly-evaluated profiles, this is not going to be
    # feasible, and tests need to be done to check that a given n_interp is enough.  But that's an
    # accuracy test, not a unit test.
    sigma_0 = 0.06
    scale = 0.02
    n_interp = 100
    oversample_fac = 1.6
    blue_limit = min(bandpass.blue_limit, bandpass_g.blue_limit)
    red_limit = max(bandpass.red_limit, bandpass_g.red_limit)
    waves = np.linspace(blue_limit, red_limit, n_interp)
    # Make a version that tests the ability of the routine to handle non-sorted / non-array input.
    tricky_waves = list(waves[::-1])

    # Make a star.
    star = galsim.Gaussian(fwhm=1.e-8) * bulge_SED

    # First, compare images that are drawn with exact and interpolated ChromaticGaussian.
    exact_psf = ChromaticGaussian(sigma_0)
    interp_psf = exact_psf.interpolate(waves, oversample_fac=oversample_fac)
    exact_obj = galsim.Convolve(star, exact_psf)
    interp_obj = galsim.Convolve(star, interp_psf)
    im_exact = exact_obj.drawImage(bandpass, scale=scale, nx=40, ny=40)
    im_interp = im_exact.copy()
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, scale=scale)
    # Note: peak value of array is around 0.3, so going to 4 decimal places is a reasonably
    # stringent test considering how different the exact vs. interpolated rendering process is.
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-4,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated')

    # And test with midpoint rule (non-default).
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, integrator='midpoint', scale=scale)
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=2.e-4,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated (midpoint)')

    # Check kimage with midpoint rule too, also non-default.
    kscale = 2*np.pi/(scale*40)
    kim_exact = exact_obj.drawKImage(bandpass, scale=kscale, nx=40, ny=40, integrator='midpoint')
    kim_interp = interp_obj.drawKImage(bandpass, scale=kscale, nx=40, ny=40, integrator='midpoint')
    np.testing.assert_allclose(
        kim_interp.array, kim_exact.array, rtol=0, atol=1e-4*kim_exact.array.real.max(),
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated (midpoint)')

    # Also trapezoidal
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, integrator='trapezoidal',
                                     scale=scale)
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=2.e-4,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated (trapz)')

    kscale = 2*np.pi/(scale*40)
    kim_exact = exact_obj.drawKImage(bandpass, scale=kscale, nx=40, ny=40, integrator='trapezoidal')
    kim_interp = interp_obj.drawKImage(bandpass, scale=kscale, nx=40, ny=40,
                                       integrator='trapezoidal')
    np.testing.assert_allclose(
        kim_interp.array, kim_exact.array, rtol=0, atol=1e-4*kim_exact.array.real.max(),
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated (trapz)')

    # Check that we can turn interpolation off and on at will.
    other_psf = interp_psf.deinterpolated
    other_obj = galsim.Convolve(star, other_psf)
    im_other = im_exact.copy()
    im_other = other_obj.drawImage(bandpass, image=im_other, scale=scale)
    # Can test to very high accuracy.
    np.testing.assert_allclose(
        im_other.array, im_exact.array, atol=1.e-8,
        err_msg='Failure to turn off interpolation in ChromaticObject')

    # Check that when an interpolated ChromaticObject is convolved with a ChromaticObject that has a
    # non-trivial surface brightness profile (i.e., a galaxy rather than a star), the image that we
    # get and the total flux normalization is correct. We are forcing it to go through the
    # evaluation of ChromaticConvolution.drawImage().
    gal = galsim.Exponential(half_light_radius = 2.*scale)
    gal = gal.shear(g2 = 0.3)
    gal = disk_SED*gal
    obj_exact = galsim.Convolve(exact_psf, gal)
    obj_interp = galsim.Convolve(interp_psf, gal)
    im_exact = obj_exact.drawImage(bandpass_g, scale=scale)
    im_interp = im_exact.copy()
    im_interp = obj_interp.drawImage(bandpass_g, image=im_interp, scale=scale)
    expected_flux = disk_SED.calculateFlux(bandpass_g)
    np.testing.assert_allclose(
        im_exact.array.sum(), expected_flux, rtol=1.e-3,
        err_msg='ChromaticObject flux is wrong when convolved with ChromaticObject '
        ' (exact calculation)')
    np.testing.assert_allclose(
        im_interp.array.sum(), expected_flux, rtol=1.e-3,
        err_msg='ChromaticObject (with interpolation) flux is wrong when convolved with '
        ' ChromaticObject')
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-4,
        err_msg='ChromaticObject results differ for interpolated vs. exact'
        ' when convolving with ChromaticObject')

    # Check that when an ChromaticObject with interpolation is convolved with a ChromaticSum, the
    # image and flux normalization is correct.
    bulge = galsim.DeVaucouleurs(half_light_radius = 1.5*scale)
    bulge = bulge_SED*bulge
    tot_gal = gal + bulge
    obj_exact = galsim.Convolve(exact_psf, tot_gal)
    obj_interp = galsim.Convolve(interp_psf, tot_gal)
    im_exact = obj_exact.drawImage(bandpass, scale=scale)
    im_interp = im_exact.copy()
    im_interp = obj_interp.drawImage(bandpass, image=im_interp, scale=scale)
    expected_flux = disk_SED.calculateFlux(bandpass) + bulge_SED.calculateFlux(bandpass)
    # Check to 2%
    np.testing.assert_allclose(
        im_exact.array.sum(), expected_flux, rtol=0.02,
        err_msg='ChromaticObject flux is wrong when convolved with ChromaticSum'
        ' (exact calculation)')
    np.testing.assert_allclose(
        im_interp.array.sum(), expected_flux, rtol=0.02,
        err_msg='ChromaticObject flux is wrong when convolved with ChromaticSum'
        ' (interpolated calculation)')
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-3,
        err_msg='ChromaticObject results differ for interpolated vs. exact'
        ' when convolving with ChromaticSum')

    # Check that InterpolatedChromaticObject.drawImage works.  Need an exact object with
    # an SED for this to work.
    exact_psf *= bulge_SED
    interp_psf = exact_psf.interpolate(waves, oversample_fac=oversample_fac)
    im_exact = exact_psf.drawImage(bandpass, scale=0.2, nx=32, ny=32)
    im_interp = interp_psf.drawImage(bandpass, scale=0.2, nx=32, ny=32)
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-3,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated')

    # Check that we can render an image with chromatic transformations directly, and with
    # interpolation.  Use a ChromaticAtmosphere just because that's easily transformed.
    atm_fwhm = 0.7
    atm_scale = 0.2
    exact_psf = galsim.ChromaticAtmosphere(
        galsim.Kolmogorov(atm_fwhm), 500., zenith_angle=0.*galsim.degrees,
        parallactic_angle=0.*galsim.degrees)
    chrom_shear = lambda w: galsim.Shear(g1=0.2+0.2*(w-500.)/500.,g2=0.) if w<1000. else \
        galsim.Shear(g1=0.4, g2=0.)
    chrom_shift_y = lambda w: scale*(w-500.)
    chrom_dilate = lambda w: 1.0+0.1*(w-500.)/500.
    exact_psf = (exact_psf
                 .shear(shear=chrom_shear)
                 .shift(dx=0.,dy=chrom_shift_y)
                 .dilate(chrom_dilate))
    # Note here we are checking the use of more difficult input wavelengths.
    interp_psf = exact_psf.interpolate(tricky_waves, oversample_fac=oversample_fac)
    exact_obj = galsim.Convolve(star, exact_psf)
    interp_obj = galsim.Convolve(star, interp_psf)
    im_exact = exact_obj.drawImage(bandpass, scale=atm_scale)
    im_interp = im_exact.copy()
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, scale=atm_scale)
    # Note: peak value of array is around 4, so going to 3 decimal places is a reasonably
    # stringent test considering how different the exact vs. interpolated rendering process is.
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-3,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated'+
        ' when including chromatic transformations')

    # Check that we can render an image with achromatic transformations applied after
    # interpolation.
    exact_psf = galsim.ChromaticAtmosphere(
        galsim.Kolmogorov(atm_fwhm), 500., zenith_angle=0.*galsim.degrees,
        parallactic_angle=0.*galsim.degrees)
    interp_psf = exact_psf.interpolate(waves, oversample_fac=oversample_fac)

    achrom_shear = galsim.Shear(g1=0.05, g2=-0.1)
    exact_psf = exact_psf.shear(shear=achrom_shear)
    exact_obj = galsim.Convolve(star, exact_psf)
    interp_psf = interp_psf.shear(shear=achrom_shear)
    interp_obj = galsim.Convolve(star, interp_psf)

    im_exact = exact_obj.drawImage(bandpass, scale=atm_scale)
    im_interp = im_exact.copy()
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, scale=atm_scale)
    # Note: peak value of array is around 4, so going to 3 decimal places is a reasonably
    # stringent test considering how different the exact vs. interpolated rendering process is.
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-3,
        err_msg='Interpolated ChromaticObject results differ for exact vs. interpolated'+
        ' when including achromatic transformations after precomputation')

    # Check that the routine does not interpolate outside of its original bounds.
    with assert_raises(galsim.GalSimError):
        obj_interp.drawImage(bandpass_z)

    # Make sure it behaves appropriately when asked to apply chromatic transformations after
    # interpolating: it should do the job properly after un-setting the interpolation.
    # As a really strict test, we will apply chromatic transformations twice (before and after
    # setting up the interpolation).  This will make sure that it's really tracking and combining
    # all transformations.
    exact_psf = galsim.ChromaticAtmosphere(
        galsim.Kolmogorov(atm_fwhm), 500., zenith_angle=0.*galsim.degrees,
        parallactic_angle=0.*galsim.degrees)
    exact_psf = \
        exact_psf.shear(shear=chrom_shear).shift(dx=0.,dy=chrom_shift_y).dilate(chrom_dilate)
    interp_psf = exact_psf.interpolate(waves, oversample_fac=oversample_fac)
    trans_exact_psf = \
        exact_psf.shear(shear=chrom_shear).shift(dx=0.,dy=chrom_shift_y).dilate(chrom_dilate)

    with assert_warns(galsim.GalSimWarning):
        trans_interp_psf = interp_psf.shear(
            shear=chrom_shear).shift(dx=0.,dy=chrom_shift_y).dilate(chrom_dilate)
    exact_obj = galsim.Convolve(star, trans_exact_psf)
    interp_obj = galsim.Convolve(star, trans_interp_psf)
    im_exact = exact_obj.drawImage(bandpass, scale=atm_scale)
    im_interp = im_exact.copy()
    im_interp = interp_obj.drawImage(bandpass, image=im_interp, scale=atm_scale)
    # Note: since the image rendering should have been done in exactly the same way (it should
    # have trashed the interpolation entirely), test to high precision.
    np.testing.assert_allclose(
        im_interp.array, im_exact.array, atol=1.e-9,
        err_msg='Did not do exact chromatic transformation by discarding interpolation')
    # Also make sure that it ditched the interpolation.
    assert not hasattr(trans_interp_psf, 'waves')

    # test alternate initialization method, "from_images()", of InterpolatedChromaticObject
    # that uses images at discrete wavelengths to initialize object.

    # check sorting is done correctly for unsorted wavelengths/images
    PSF_exact = ChromaticGaussian(sigma_0)
    PSF = PSF_exact.interpolate(tricky_waves, oversample_fac=oversample_fac)
    tricky_images = list(PSF.ims[::-1])
    int_psf = galsim.InterpolatedChromaticObject.from_images(tricky_images, tricky_waves,_force_stepk =PSF.stepk_vals,
                                                              _force_maxk = PSF.maxk_vals)
    np.testing.assert_allclose(int_psf.waves, waves, atol=0,
      err_msg='InterpolatedChromaticObject from_images initialization fails to sort wavelengths')

    for i in range(len(int_psf.ims)):
        np.testing.assert_allclose(int_psf.ims[i].array, PSF.ims[i].array, atol=1e-17,
      err_msg='InterpolatedChromaticObject from_images initialization fails to sort images correctly')

    # check error messages from images with a non-PixelScale wcs, inconsistent pixel scales and image dimensions
    incorrect_ims = PSF.ims.copy()
    incorrect_ims[0].wcs = None
    # wcs = None
    assert_raises(galsim.GalSimValueError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)
    affine_wcs = galsim.AffineTransform(1.0, 0.1, 0.1, 1.0)
    incorrect_ims[0].wcs = affine_wcs
    # non-PixelScale wcs
    assert_raises(galsim.GalSimValueError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)
    incorrect_ims[0].wcs = incorrect_ims[1].wcs
    incorrect_ims[0].scale += 0.01
    # incosnistent pixel scales
    assert_raises(galsim.GalSimValueError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)
    incorrect_ims[0].scale -= 0.01
    smaller_image = incorrect_ims[0].array[1:]
    incorrect_img = galsim.Image(smaller_image, scale=incorrect_ims[0].scale)
    incorrect_ims[0] = incorrect_img
    # incosnistent image dimensions
    assert_raises(galsim.GalSimValueError, galsim.InterpolatedChromaticObject.from_images, incorrect_ims, PSF.waves)

    # check input images are correctly intialized from underscored function directly
    PSF = PSF_exact.interpolate(waves, oversample_fac=oversample_fac)
    int_psf = galsim.InterpolatedChromaticObject._from_images(PSF.ims, list(PSF.waves),_force_stepk =PSF.stepk_vals,
                                                              _force_maxk = PSF.maxk_vals)
    assert isinstance(int_psf.waves, np.ndarray)
    for i in range(len(int_psf.ims)):
        np.testing.assert_allclose(int_psf.ims[i].array, PSF.ims[i].array, atol=1e-17,
      err_msg='InterpolatedChromaticObject from_images initialization fails to initialize correct images')
    test_obj = galsim.Convolve(int_psf, star)
    true_obj = galsim.Convolve(PSF, star)
    true_img = true_obj.drawImage(bandpass, scale=scale, method = 'auto')
    im_interp = true_img.copy()
    test_img = test_obj.drawImage(bandpass, image=im_interp, scale=scale)
    # check drawing objects is identical when using same stepk and maxk as images
    np.testing.assert_allclose(test_img.array, true_img.array, atol=1e-17,
            err_msg='InterpolatedChromaticObject from_images initialization fails to reproduce default init. images')

    # without specifying the same stepk and maxk for each image, stepk and maxk are caclulated
    # based on the input image pixel scale and dimensions and have no wavelength dependance.
    int_psf = galsim.InterpolatedChromaticObject.from_images(PSF.ims, PSF.waves)
    test_obj = galsim.Convolve(int_psf, star)
    test_img = test_obj.drawImage(bandpass, image=im_interp, scale=scale, method = 'auto')
    # check images are within 0.01% of the total flux, the flux and maximum within 0.1%
    np.testing.assert_allclose(test_img.array, true_img.array, atol=1e-4*np.sum(true_img.array),
        err_msg='InterpolatedChromaticObject images differ when using from_images initialization with default stepk, maxk')
    np.testing.assert_allclose(test_img.array.sum(), true_img.array.sum(), rtol=1e-3,
        err_msg='InterpolatedChromaticObject images flux differ when using from_images initialization with default stepk, maxk')
    np.testing.assert_allclose(test_img.array.max(), true_img.array.max(), rtol=1e-3,
        err_msg='InterpolatedChromaticObject images maximum differ when using from_images initialization with default stepk, maxk')
    # check moments
    truth_mom = galsim.hsm.FindAdaptiveMom(true_img)
    test_mom = galsim.hsm.FindAdaptiveMom(test_img)

    np.testing.assert_allclose(test_mom.moments_amp,
                               truth_mom.moments_amp,
                               rtol=1e-3, atol=0,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments amplitude ')
    np.testing.assert_allclose(test_mom.moments_centroid.x,
                               truth_mom.moments_centroid.x,
                               rtol=0., atol=1e-2,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.x ')
    np.testing.assert_allclose(test_mom.moments_centroid.y,
                               truth_mom.moments_centroid.y,
                               rtol=0., atol=1e-2,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.y ')
    np.testing.assert_allclose(test_mom.moments_sigma,
                               truth_mom.moments_sigma,
                               rtol=1e-3, atol=0,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments sigma ')
    np.testing.assert_allclose(test_mom.observed_shape.g1,
                               truth_mom.observed_shape.g1,
                               rtol=0, atol=1e-4,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments g1 ')
    np.testing.assert_allclose(test_mom.observed_shape.g2,
                               truth_mom.observed_shape.g2,
                               rtol=0, atol=1e-4,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments g2 ')

    # check that flux macthes when convolved with non-trivial light profile
    gal = galsim.Exponential(half_light_radius = 2.*scale)
    gal = gal.shear(g2 = 0.3)
    gal = disk_SED*gal
    obj_exact = galsim.Convolve(PSF, gal)
    obj_interp = galsim.Convolve(int_psf, gal)
    true_img = obj_exact.drawImage(bandpass_g, scale=scale)
    im_interp = true_img.copy()
    test_img = obj_interp.drawImage(bandpass_g, image=im_interp, scale=scale)
    np.testing.assert_allclose(
        test_img.array, true_img.array, atol=1.e-4,
        err_msg='Interpolated ChromaticObject from_images initialization differs from default when drawing galaxy')
    # check moments
    truth_mom = galsim.hsm.FindAdaptiveMom(true_img)
    test_mom = galsim.hsm.FindAdaptiveMom(test_img)

    np.testing.assert_allclose(test_mom.moments_amp,
                               truth_mom.moments_amp,
                               rtol=1e-3, atol=0,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments amplitude ')
    np.testing.assert_allclose(test_mom.moments_centroid.x,
                               truth_mom.moments_centroid.x,
                               rtol=0., atol=1e-2,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.x ')
    np.testing.assert_allclose(test_mom.moments_centroid.y,
                               truth_mom.moments_centroid.y,
                               rtol=0., atol=1e-2,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments centroid.y ')
    np.testing.assert_allclose(test_mom.moments_sigma,
                               truth_mom.moments_sigma,
                               rtol=1e-3, atol=0,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments sigma ')
    np.testing.assert_allclose(test_mom.observed_shape.g1,
                               truth_mom.observed_shape.g1,
                               rtol=0, atol=1e-4,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments g1 ')
    np.testing.assert_allclose(test_mom.observed_shape.g2,
                               truth_mom.observed_shape.g2,
                               rtol=0, atol=1e-4,
                err_msg='InterpolatedChromaticObject from_images init. differs in moments g2 ')


@timer
def test_ChromaticOpticalPSF(run_slow):
    """Test the ChromaticOpticalPSF functionality."""
    import time

    # For ChromaticOpticalPSF, exact evaluation is too slow for routine unit tests.  So, for
    # this unit test, we use an interpolated version only.  The tests of
    # ChromaticObject in the previous unit test should be enough to ensure that exact
    # and interpolated evaluation match in general (given reasonable settings).

    # First, compare the interpolated result with saved, exact results.
    # Note that if the saved results file cannot be found, then this function will (slowly)
    # generate the exact image.  Note that this may need to be done outside of pytest to avoid
    # the default pytest time limit.

    # Note that exact results will have to be regenerated if any of the bandpasses or other
    # parameters defined here are changed.  Because of the parameters chosen here, there is a lot of
    # non-trivially complex structure in the PSFs, so this is a stringent test.
    aberrations = np.zeros(12)
    aberrations[4] = 40. # nm
    aberrations[7] = 20. # nm
    lam = 600. # nm
    aberrations /= lam
    diam = 2.4 # meters
    obscuration = 0.18
    nstruts = 4
    scale = 0.02
    n_interp = 15
    oversample_fac = 5.0
    waves = np.linspace(bandpass.blue_limit, bandpass.red_limit, n_interp)

    # from pylab import *
    # wexact = np.union1d(bandpass.wave_list, disk_SED.wave_list)
    # wexact = wexact[(wexact>bandpass.blue_limit) & (wexact<bandpass.red_limit)]
    # plot(wexact, bandpass(wexact)*disk_SED(wexact), marker='o', c='b')
    # plot(waves, bandpass(waves)*disk_SED(waves), ls=' ', marker='o', c='r', markersize=10)
    # show()

    star = galsim.Gaussian(fwhm=1.e-8) * disk_SED
    psf = galsim.ChromaticOpticalPSF(lam=lam, diam=diam, aberrations=aberrations,
                                     obscuration=obscuration, nstruts=nstruts)
    check_pickle(psf)

    # Test some other valid initialization options
    pupil_plane_im = os.path.join('Optics_comparison_images', 'sample_pupil_rolled.fits')
    psf2 = galsim.ChromaticOpticalPSF(lam=lam, diam=diam,
                                      astig1=0.2, astig2=-0.1, coma1=0.03, coma2=0.04,
                                      fft_sign='-', geometric_shooting=True,
                                      pupil_plane_im=pupil_plane_im)
    check_pickle(psf2)

    # And some invalid options.
    with assert_raises(galsim.GalSimIncompatibleValuesError):
        galsim.ChromaticOpticalPSF(lam=lam, diam=diam, aberrations=aberrations, lam_over_diam=0.02)
    with assert_raises(galsim.GalSimValueError):
        galsim.ChromaticOpticalPSF(lam=lam, diam=diam, fft_sign=0)
    with assert_raises(TypeError):
        galsim.ChromaticOpticalPSF(lam=lam, diam=diam, aper=psf2._aper, obscuration=0.3)

    # We will compare to this reference image below.  Make it if it doesn't exist yet.
    if not os.path.isfile(os.path.join(refdir, 'r_exact.fits')):
        import warnings
        warnings.warn("Could not find file r_exact.fits, so generating it from scratch.  This "
                      "should only happen if you intentionally deleted the file in order to "
                      "regenerate it!")

        # Generate exact results inside this if-block.
        # This block took ~30 seconds to run on a ~2016 Macbook Pro, nearly all in the image
        # rendering process.  In contrast, the ChromaticOpticalPSF with interpolation that is used
        # for this unit test takes about ~9 seconds to initialize, and ~1s for the image rendering
        # process.  Obviously, if many images are to be rendered after incurring the overhead of
        # initializing this object, the interpolated calculation leads to a huge savings compared to
        # doing the exact calculation each time.
        t2 = time.time()
        obj = galsim.Convolve(star, psf)
        im_r_ref = obj.drawImage(bandpass, scale=scale)
        im_r_ref.write(os.path.join(refdir, 'r_exact.fits'))
        t3 = time.time()
        print("Time to draw ChromaticOpticalPSF: {0}s".format(t3-t2))

    t4 = time.time()
    psf = psf.interpolate(waves, oversample_fac=oversample_fac)
    t5 = time.time()
    print("Time to initialize InterpolatedChromaticObject: {0}s".format(t5-t4))
    obj = galsim.Convolve(star, psf)

    if run_slow:
        # This is slow, but it is worth testing the pickling of InterpolatedChromaticObjects.
        check_pickle(psf)
    else:
        repr(psf)

    im_r_ref = galsim.fits.read(os.path.join(refdir, 'r_exact.fits'))
    im_r = im_r_ref.copy()
    t6 = time.time()
    obj.drawImage(bandpass, image=im_r, scale=scale)
    t7 = time.time()
    print("Time to draw InterpolatedChromaticObject: {0}s".format(t7-t6))
    printval(im_r, im_r_ref)
    # Check that arrays agree to within 1e-4 of integrated flux
    np.testing.assert_allclose(im_r.array, im_r_ref.array, atol=1e-4*im_r_ref.array.sum())
    # Check sums to 1e-3 of sum
    np.testing.assert_allclose(im_r.array.sum(), im_r_ref.array.sum(), atol=1e-3*im_r.array.sum())
    # Check peak to 1e-3 of peak
    np.testing.assert_allclose(im_r.array.max(), im_r_ref.array.max(), atol=1e-3*im_r.array.max())

    im_r_ref /= im_r_ref.array.max()
    im_r /= im_r.array.max()
    # Test nearly passes at decimal=4, but 0.08% of pixels disagree.  However, decimal=3 after
    # normalization with peak flux is still very good.
    np.testing.assert_array_almost_equal(
        im_r.array, im_r_ref.array, decimal=3,
        err_msg='Interpolated ChromaticOpticalPSF results disagree with reference in r band')

    # Finally, check that flux normalization is preserved when we convolve with a chromatic object.
    gal = galsim.Exponential(half_light_radius = 2.*scale)
    gal = gal.shear(g2 = 0.3)
    gal = disk_SED*gal
    obj_conv = galsim.Convolve(psf, gal)
    im = obj_conv.drawImage(bandpass, scale=scale)
    expected_flux = disk_SED.calculateFlux(bandpass)

    assert np.isclose(im.array.sum(), expected_flux, rtol=0.02), \
        "ChromaticObject flux is wrong when convolved with ChromaticOpticalPSF " \
        " (interpolated calculation)"


@timer
def test_ChromaticAiry():
    """Test the ChromaticAiry functionality."""
    # First, compare the interpolated result with saved, exact results.
    # Exact results were generated using the following code, sitting in this directory:
    #
    # import galsim
    # import os
    # import numpy as np
    #
    # path, filename = os.path.split(__file__)
    # datapath = os.path.abspath(os.path.join(path, "../examples/data/"))
    # bandpass = (galsim.Bandpass(os.path.join(datapath, 'LSST_r.dat'))
    #             .truncate(relative_throughput=1e-3)
    #             .thin(rel_err=1e-3))
    # disk_SED = (galsim.SED(os.path.join(datapath, 'CWW_Sbc_ext.sed'), wave_type='ang')
    #             .thin(rel_err=1e-3)
    #             .withFluxDensity(target_flux_density=0.3, wavelength=500.0))
    #
    # star = galsim.Gaussian(sigma=1.e-8)*disk_SED
    #
    # lam = 750. # nm
    # diam = 3.1 # meters
    # obscuration = 0.11
    # scale = 0.02
    #
    # psf = galsim.ChromaticAiry(lam=lam, diam=diam, obscuration=obscuration)
    # obj = galsim.Convolve(psf, star)
    # im_r = obj.drawImage(bandpass, scale=scale)
    # im_r.write('./chromatic_reference_images/r_exact_Airy.fits')
    #
    # Note that exact results will have to be regenerated if any of the bandpasses or other
    # parameters defined here are changed.  For example, had to regenerate on #590 because of new
    # version of r bandpass.

    # Define parameters:
    lam = 750. # nm
    diam = 3.1 # meters
    obscuration = 0.11
    scale = 0.02

    psf = galsim.ChromaticAiry(lam=lam, diam=diam, obscuration=obscuration)
    check_pickle(psf)

    # Generate a reference image
    star = galsim.Gaussian(fwhm=1.e-8) * disk_SED
    obj = galsim.Convolve(psf, star)
    im_r = obj.drawImage(bandpass, scale=scale)

    # Initialize object in different way, make sure results are identical:
    lam_over_diam = (1.e-9*lam/diam)*galsim.radians
    psf = galsim.ChromaticAiry(lam=lam, lam_over_diam=lam_over_diam/galsim.arcsec,
                               obscuration=obscuration)
    obj = galsim.Convolve(psf, star)
    im_r_2 = im_r.copy()
    obj.drawImage(bandpass, image=im_r_2, scale=scale)
    np.testing.assert_array_almost_equal(
        im_r_2.array, im_r.array, decimal=8,
        err_msg='Inconsistent ChromaticAiry image when initializing a different way')

    with assert_raises(galsim.GalSimIncompatibleValuesError):
        galsim.ChromaticAiry(lam=lam, diam=diam, lam_over_diam=lam_over_diam/galsim.arcsec)

    # Also check evaluation at a single wavelength.
    chromatic_psf_400 = psf.evaluateAtWavelength(400.)
    new_lam_over_diam = (1.e-9*400/diam)*galsim.radians
    exact_psf_400 = galsim.Airy(lam_over_diam=new_lam_over_diam/galsim.arcsec,
                                obscuration=obscuration)
    chr_im = chromatic_psf_400.drawImage(scale=scale)
    exact_im = chr_im.copy()
    exact_psf_400.drawImage(image=exact_im, scale=scale)
    np.testing.assert_array_almost_equal(
        chr_im.array, exact_im.array, decimal=8,
        err_msg='ChromaticAiry evaluated at a single wavelength is not as expected')

    # Finally, check that flux normalization is preserved when we convolve with a chromatic object.
    gal = galsim.Exponential(half_light_radius = 2.*scale)
    gal = gal.shear(g2 = 0.3)
    gal = disk_SED*gal
    obj_conv = galsim.Convolve(psf, gal)
    im = obj_conv.drawImage(bandpass, scale=scale)
    expected_flux = disk_SED.calculateFlux(bandpass)
    frac_diff_exact = abs(im.array.sum()/expected_flux-1.0)
    # Check to 1%
    np.testing.assert_almost_equal(
        frac_diff_exact, 0.0, decimal=2,
        err_msg='ChromaticObject flux is wrong when convolved with ChromaticAiry')

@timer
def test_phot(run_slow):
    """Test photon shooting with various chromatic PSFs.
    """
    import time

    # Use an emission line SED to exaggerate the effect of the wavelength dependence of the
    # Airy profile.  Otherwise we need a lot of photons to notice the effect.
    bandpass = galsim.Bandpass(galsim.LookupTable([500,1000], [1,1], 'linear'), wave_type='nm')
    sed = galsim.EmissionLine(510, fwhm=5) + galsim.EmissionLine(975, fwhm=5)
    flux = 1.e6
    gal_achrom = galsim.Sersic(n=2.8, half_light_radius=0.03, flux=flux)
    gal = (gal_achrom * sed).withFlux(flux, bandpass=bandpass)

    # 1. ChromaticAiry
    diam = 3.1 # meters
    obscuration = 0.11
    psf1 = galsim.ChromaticAiry(lam=bandpass.effective_wavelength, diam=diam,
                                obscuration=obscuration)
    # Do this first so we can see the real timing for photon shooting separate from the setup time.
    t0 = time.time()
    psf_achrom = galsim.Airy(lam=bandpass.effective_wavelength, diam=diam, obscuration=obscuration)
    psf_achrom.shoot(1)
    t1 = time.time()
    print('Airy setup time = ',t1-t0)
    print('psf hlr = ',psf_achrom.calculateHLR())
    # hlr = 0.04 arcsec for this combination.  So our galaxy needs to be ~this small for the
    # effects of changing the PSF to be important.  Even smaller is better.

    # 2. ChromaticAiry with wavelength-dependent transformation
    psf2 = galsim.Transform(psf1,
                            jac=lambda w: np.array([[1+0.1*(w/1000)**2, 0.1*(w-1000)/500],
                                                    [-0.15*(w-1000)/500, 1-0.08*(w/1000)**2]]),
                            offset=lambda w: np.array([1.e-4*(w-500), 1.e-4*(w-1000)]),
                            flux_ratio=lambda w: w/1000)

    # 3. Moffat with dimensionless SED
    psf3 = galsim.Moffat(beta=2.5, fwhm=0.12)
    psf3 = psf3 * galsim.SED('((wave-500)/500)**3', wave_type='nm', flux_type='1')

    # 4. Moffat just pretending to be chromatic
    psf4 = galsim.ChromaticObject(galsim.Moffat(beta=2.5, fwhm=0.12))

    # 5. ChromaticOpticalPSF with geometric shooting
    aberrations = np.array([0,0,0,0, 0.008, -0.005, -0.005, -0.002])
    psf5 = galsim.ChromaticOpticalPSF(lam=bandpass.effective_wavelength, diam=diam,
                                      aberrations=aberrations, obscuration=obscuration,
                                      geometric_shooting=True)

    # 6. ChromaticOpticalPSF without geometric shooting
    # Use struts this time.  (geometric_shooting=True doesn't handle struts well.)
    aberrations = np.array([0,0,0,0, 0.03, -0.05, -0.05, -0.02, 0.04, 0.02, -0.02, -0.03, 0.02])
    psf6 = galsim.ChromaticOpticalPSF(lam=bandpass.effective_wavelength, diam=diam,
                                      aberrations=aberrations, obscuration=obscuration,
                                      nstruts=6, geometric_shooting=False)

    # 7. InterpolatedChromaticObject
    psf7 = psf6.interpolate(waves=[500, 700, 1000])

    # 8. ChromaticSum
    # Most useful as a way to linearly interpolate between two (or more) values, so do that.
    ab1 = np.array([0,0,0,0, 0.008, -0.005, -0.005, -0.002])
    ab2 = np.array([0,0,0,0, 0.003, 0.003, -0.007, -0.001])
    psf8a = galsim.ChromaticOpticalPSF(lam=bandpass.effective_wavelength, diam=diam,
                                       aberrations=ab1, obscuration=obscuration,
                                       geometric_shooting=True)
    psf8b = galsim.ChromaticOpticalPSF(lam=bandpass.effective_wavelength, diam=diam,
                                       aberrations=ab2, obscuration=obscuration,
                                       nstruts=6, geometric_shooting=True)
    psf8 = 0.7 * psf8a + 0.3 * psf8b

    for psf in [psf1, psf2, psf3, psf4, psf5, psf6, psf7, psf8]:
        print('psf = ',psf)
        atol = 4.e-4
        if psf in [psf5, psf6, psf7, psf8]:
            atol = 6e-4   # OpticalPSF doesn't match quite as well as the others.
        rng = galsim.BaseDeviate(1234)

        # First draw with FFT
        obj = galsim.Convolve(gal, psf)
        pixel_scale = 0.01
        im1 = galsim.ImageD(100, 100, scale=pixel_scale)
        t0 = time.time()
        obj.drawImage(bandpass, image=im1)
        t1 = time.time()
        print('fft time = ',t1-t0)
        print('im1.max,sum = ', im1.array.max(), im1.array.sum())

        if run_slow:
            # Now the direct photon shooting method
            t0 = time.time()
            # This is the old way that photon shooting used to work.  The new way will be tested
            # below.  But since it now uses photon_ops, we'll wait to test that last.
            effective_psf = galsim.ChromaticConvolution._get_effective_prof(
                    psf*gal.sed, bandpass, integrator='trapezoidal', gsparams=psf.gsparams,
                    iimult=None)
            temp_obj = galsim.Convolve(gal_achrom/flux,effective_psf)
            im2 = temp_obj.drawImage(image=im1.copy(), method='phot', rng=rng)
            t1 = time.time()
            print('old method phot time = ',t1-t0)
            print('max diff/flux = ',np.max(np.abs(im1.array-im2.array)/flux))
            print('im2.max,sum = ', im2.array.max(), im2.array.sum())
            np.testing.assert_allclose(im2.array/flux, im1.array/flux, atol=atol)

            if psf is not psf4:
                # Make sure we would notice if psf wasn't applying the wavelength scaling properly.
                # Note: This test is why we're using such a crazy bandpass and sed here.
                # With a realistic SED and bandpass, the difference we're looking for is too subtle
                # to see with only 10^6 photons.
                psf_achrom = psf.evaluateAtWavelength(bandpass.effective_wavelength)
                achrom = galsim.Convolve(gal_achrom, psf_achrom)
                t0 = time.time()
                im2b = achrom.drawImage(image=im1.copy(), method='phot', rng=rng)
                t1 = time.time()
                print('achrom phot time = ',t1-t0)
                print('max diff/flux = ',np.max(np.abs(im1.array-im2b.array)/flux))
                print('im2b.max,sum = ', im2b.array.max(), im2b.array.sum())
                # This is about 1.5e-3.  So ~5x the tolerance we're using for the correct method.
                with assert_raises(AssertionError):
                    np.testing.assert_allclose(im2b.array/flux, im1.array/flux, atol=atol)

            # Now using photon_ops with both wave_sampler and psf.
            wave_sampler = galsim.WavelengthSampler(sed, bandpass)
            t0 = time.time()
            im3 = gal_achrom.drawImage(image=im1.copy(), method='phot', rng=rng,
                                       photon_ops=[wave_sampler, psf])
            t1 = time.time()
            print('wave_sampler time = ',t1-t0)
            print('max diff/flux = ',np.max(np.abs(im1.array-im3.array)/flux))
            print('im3.max,sum = ', im3.array.max(), im3.array.sum())
            np.testing.assert_allclose(im3.array/flux, im1.array/flux, atol=atol)

        # Error if wavelengths aren't set.
        with assert_raises(galsim.GalSimError):
            gal_achrom.drawImage(image=im1.copy(), method='phot', rng=rng, photon_ops=[psf])

        # Finally, the chromatic drawImage function should handle the wavelength sampling for us.
        # First do this with just the galaxy as the driver.
        # (Also add a gratuitous pupil_sampler, time_sampler to test the case where u,v,t are set.)
        pupil_sampler = galsim.PupilImageSampler(diam=diam, lam=bandpass.effective_wavelength)
        time_sampler = galsim.TimeSampler()
        t0 = time.time()
        im4 = gal.drawImage(bandpass, image=im1.copy(), method='phot', rng=rng,
                            photon_ops=[pupil_sampler, time_sampler, psf])
        t1 = time.time()
        print('auto wave time = ',t1-t0)
        print('max diff/flux = ',np.max(np.abs(im1.array-im4.array)/flux))
        print('im4.max,sum = ', im4.array.max(), im4.array.sum())
        np.testing.assert_allclose(im4.array/flux, im1.array/flux, atol=atol)

        # Now let the ChromaticConvolution reorganize this for us.
        t0 = time.time()
        calc_flux = obj.calculateFlux(bandpass)
        im5 = obj.drawImage(bandpass, image=im1.copy(), method='phot', rng=rng,
                            n_photons=flux, save_photons=True)
        t1 = time.time()
        print('regular phot time = ',t1-t0)
        print('max diff/flux = ',np.max(np.abs(im1.array-im5.array)/flux))
        print('im5.max,sum = ', im5.array.max(), im5.array.sum())
        np.testing.assert_allclose(im5.array/flux, im1.array/flux, atol=atol)

        assert hasattr(im5, 'photons')
        assert len(im5.photons) == flux
        assert np.allclose(np.sum(im5.photons.flux), calc_flux, rtol=1.e-3)

        # Check that n_photons=1 doesn't do something bad, like produce NaNs.
        (obj/flux).drawImage(bandpass, image=im5, method='phot', rng=rng, n_photons=1,
                             add_to_image=True)
        print('im5.max,sum => ', im5.array.max(), im5.array.sum())
        np.testing.assert_allclose(im5.array/flux, im1.array/flux, atol=atol)

        # Different code path when scale=1, since no jac (in some cases at least)
        im6 = obj.drawImage(bandpass, nx=50, ny=50, scale=1., method='phot', rng=rng)
        print('im6.max,sum = ', im6.array.max(), im6.array.sum())
        np.testing.assert_allclose(im6.array.sum(), im1.array.sum(), rtol=0.1)

    # Invalid to shoot photons outside of InterpolatedChromaticObject's range.
    # Note: It requires some care to avoid all the other range checks and actually hit the
    #       one in InterpolatedChromaticObject._shoot.
    rng = galsim.BaseDeviate(1234)
    bp2 = galsim.Bandpass(galsim.LookupTable([400,601], [1,1], 'linear'), wave_type='nm')
    sed2 = galsim.EmissionLine(415, fwhm=5) + 2*galsim.EmissionLine(515, fwhm=5)
    gal2 = (gal_achrom * sed2).withFlux(flux, bandpass=bp2)
    obj2 = galsim.Convolve(gal2, psf7)
    with assert_raises(galsim.GalSimRangeError):
        obj2.drawImage(bp2, image=im5, method='phot', rng=rng, n_photons=10)

@timer
def test_low_flux_phot():
    """ Check that objects with 0 realized flux don't have problems.
    """

    bandpass = galsim.Bandpass("LSST_r.dat", wave_type="nm")
    sed = galsim.SED('vega.txt', 'nm', 'flambda').thin(rel_err=1.e-2)
    sed = sed.withFlux(1.e-5, bandpass)

    base_psf = galsim.Gaussian(fwhm=0.7)
    psf = galsim.ChromaticAtmosphere(
        base_psf,
        700,
        alpha=-0.3,
        zenith_angle=0 * galsim.degrees,
        parallactic_angle=0 * galsim.degrees,
    )

    observed = galsim.Convolve(psf, sed * galsim.Exponential(half_light_radius=0.3))
    rng = galsim.BaseDeviate(1234)
    image = observed.drawImage(nx=53, ny=53, bandpass=bandpass, method="phot", rng=rng)
    np.testing.assert_array_equal(image.array, 0.)


@timer
def test_chromatic_fiducial_wavelength():
    """ Check that chromatic code can handle profiles with flux(effective_wavelength) = 0.
    """
    waves = np.arange(500., 600.1, 10.)
    blue_flux = waves < 550.0
    red_flux = waves > 550.0
    bp = galsim.Bandpass(galsim.LookupTable(waves, waves**0, interpolant='linear'), 'nm')
    blue_sed = galsim.SED(galsim.LookupTable(waves, blue_flux, interpolant='linear'),
                          'nm', 'flambda')
    red_sed = galsim.SED(galsim.LookupTable(waves, red_flux, interpolant='linear'),
                         'nm', 'flambda')

    gal1 = galsim.Gaussian(fwhm=1) * blue_sed
    gal2 = galsim.Gaussian(fwhm=1) * red_sed
    img1 = gal1.drawImage(bp)
    img2 = gal2.drawImage(bp)
    # Didn't see an obvious np.testing method for checking finiteness, so just use assert.
    assert np.isfinite(img1.array.sum()), "drawImage failed to identify fiducial wavelength"
    assert np.isfinite(img2.array.sum()), "drawImage failed to identify fiducial wavelength"

    # Pathalogical sed that is zero across the band.
    bad_sed = galsim.SED(galsim.LookupTable([300,498,499,601,602,800],
                                            [  1,  1,  0,  0,  1,  1], 'linear'), 'nm', 'flambda')
    gal3 = galsim.Gaussian(fwhm=1) * bad_sed
    with assert_raises(galsim.GalSimError):
        gal3.drawImage(bp)

@timer
def test_convolve_pixel():
    """Check that convolving a chromatic object with a Pixel works
    """
    # In response to issue #1302
    # The problem here had been that if some of the components of the convolution were
    # inseparable, then when separating out the separable and inseparable components,
    # the code called `obj._fiducual_profile(bandpass)` on the separable ones.
    # However, if any of these were simple GSObjects, this gave an error, since we didn't have
    # that method for GSObject.  The fix was to add a trivial version of that for GSObject.

    bandpass = galsim.Bandpass("LSST_r.dat", wave_type="nm").thin(rel_err=1.e-2)
    sed = galsim.SED('vega.txt', 'nm', 'flambda').thin(rel_err=1.e-2)
    insep_psf = galsim.ChromaticAiry(lam=550, diam=0.1)
    sep_psf = galsim.Kolmogorov(fwhm=0.65) * galsim.SED('(wave/500)**-0.2', 'nm', '1')
    scale = 0.2
    pixel = galsim.Pixel(scale)

    for psf in [sep_psf, insep_psf]:
        star = galsim.DeltaFunction() * sed
        star = star.withFlux(1000., bandpass)
        eff_psf = galsim.Convolve(psf, pixel)

        im1 = galsim.Convolve(star, psf).drawImage(bandpass, scale=scale, nx=32, ny=32)
        im2 = galsim.Convolve(star, eff_psf).drawImage(bandpass, scale=scale, nx=32, ny=32,
                                                       method='no_pixel')
        np.testing.assert_allclose(im1.array, im2.array, atol=1.e-5)

        galaxy = galsim.Exponential(half_light_radius=1.3) * sed
        galaxy = galaxy.withFlux(1000., bandpass)
        im3 = galsim.Convolve(galaxy, psf).drawImage(bandpass, scale=scale, nx=32, ny=32)
        im4 = galsim.Convolve(galaxy, eff_psf).drawImage(bandpass, scale=scale, nx=32, ny=32,
                                                         method='no_pixel')
        np.testing.assert_allclose(im3.array, im4.array, atol=1.e-5)

@timer
def test_chromatic_image_setup():
    """Test ability for chromatic drawImage to setup output image."""
    psf = galsim.ChromaticAiry(lam=550, diam=0.1)
    gal1 = galsim.Sersic(n=1, half_light_radius=1.0) * bulge_SED
    gal2 = galsim.Sersic(n=2, half_light_radius=0.5) * disk_SED

    # Just going to try drawing a few different combinations of profiles to make sure that
    # the automatic image construction logic doesn't crash.  Most possibilities effectively get
    # tested in other scripts above anyway, so just focus on possibilities related to known previous
    # failures here; specifically, drawing a convolution of two inseparable profiles while
    # specifying `nx`, `ny`, and `scale` as keywords.
    img = galsim.Convolve(gal1+gal2, psf).drawImage(bandpass, nx=32, ny=32, scale=0.2)
    bds = galsim.BoundsI(1, 32, 1, 32)
    img2 = galsim.Convolve(gal1+gal2, psf).drawImage(bandpass, bounds=bds, scale=0.2)
    np.testing.assert_array_equal(img.array.shape, (32, 32),
                                  "Got wrong size output image using nx=, ny= keywords.")
    np.testing.assert_array_equal(img2.array.shape, (32, 32),
                                  "Got wrong size output image using bounds= keyword.")
    np.testing.assert_almost_equal(img.scale, 0.2, 9,
                                   "Got wrong output image scale using nx=, ny= keywords.")
    np.testing.assert_almost_equal(img2.scale, 0.2, 9,
                                   "Got wrong output image scale using bounds= keyword.")


@timer
def test_convolution_of_spectral():
    """Test that convolution of >1 spectral ChromaticObjects raises a ValueError."""
    cgal1 = galsim.Gaussian(fwhm=1) * bulge_SED              # spectral
    cgal2 = galsim.Gaussian(fwhm=2)                          # dimensionless
    cgal3 = galsim.ChromaticObject(galsim.Gaussian(fwhm=3))  # dimensionless

    assert cgal1.spectral
    assert not cgal2.spectral
    assert not cgal3.spectral
    assert not cgal1.dimensionless
    assert cgal2.dimensionless
    assert cgal3.dimensionless

    # These should be okay
    galsim.Convolve(cgal1, cgal2)
    galsim.Convolve(cgal1, cgal3)
    galsim.Convolve(cgal1, cgal2, cgal3)

    # This should raise a ValueError
    assert_raises(ValueError, galsim.Convolve, cgal1, cgal1)
    assert_raises(ValueError, galsim.Convolve, cgal1, cgal1, cgal1)
    assert_raises(ValueError, galsim.Convolve, cgal1, cgal1, cgal2, cgal3)


def check_chromatic_invariant(obj, bps=None, waves=None):
    """ Helper function to check that ChromaticObjects satisfy intended invariants.
    """
    if bps is None:
        # load a filter
        bppath = os.path.join(galsim.meta_data.share_dir, 'bandpasses')
        bandpass = (galsim.Bandpass(os.path.join(bppath, 'LSST_r.dat'), 'nm')
                    .truncate(relative_throughput=1e-3)
                    .thin(rel_err=1e-3))
        bps = [bandpass]

    if waves is None:
        waves = [500.]

    assert isinstance(obj.wave_list, np.ndarray)
    assert isinstance(obj.separable, bool)
    assert isinstance(obj.interpolated, bool)
    assert isinstance(obj.deinterpolated, (galsim.ChromaticObject, galsim.GSObject))

    for wave in waves:
        desired = obj.sed(wave)
        # Since InterpolatedChromaticObject.evaluateAtWavelength involves actually drawing an
        # image, which implies flux can be lost off of the edges of the image, we don't expect
        # its accuracy to be nearly as good as for other objects.
        decimal = 2 if obj.interpolated else 7
        np.testing.assert_almost_equal(obj.evaluateAtWavelength(wave).flux, desired,
                                       decimal)
        # Don't bother trying to draw a deconvolution.
        if isinstance(obj, galsim.ChromaticDeconvolution):
            continue
        np.testing.assert_allclose(
                obj.evaluateAtWavelength(wave).drawImage().array.sum(dtype=float),
                desired,
                rtol=1e-2)

    if obj.sed.spectral:
        for bp in bps:
            calc_flux = obj.calculateFlux(bp)
            np.testing.assert_equal(obj.sed.calculateFlux(bp), calc_flux)
            np.testing.assert_allclose(calc_flux,
                                       obj.drawImage(bp).array.sum(dtype=float), rtol=1e-2)
            # Also try manipulating exptime and area.
            np.testing.assert_allclose(
                    calc_flux * 10,
                    obj.drawImage(bp, exptime=5, area=2).array.sum(dtype=float), rtol=1e-2)

            assert_raises(galsim.GalSimSEDError, galsim.Deconvolve, obj)
            assert_raises(galsim.GalSimSEDError, galsim.AutoConvolve, obj)
            assert_raises(galsim.GalSimSEDError, galsim.AutoCorrelate, obj)
            assert_raises(galsim.GalSimSEDError, galsim.FourierSqrt, obj)

        try:
            obj = copy.copy(obj)
            obj.sed = galsim.SED('1', 'nm', '1')
        except AttributeError:
            return
    if isinstance(obj, galsim.GSObject):
        with assert_raises(galsim.GalSimSEDError):
            (obj * galsim.SED(1, 'nm', '1')).drawImage(bps[0])
        return

    # When made with the same gsparams, it returns itself
    print('obj = ',obj)
    assert obj.withGSParams(obj.gsparams) is obj
    alt_gsp = galsim.GSParams(xvalue_accuracy=0.23, folding_threshold=7.3e-4)
    obj_alt = obj.withGSParams(alt_gsp)
    obj_alt2 = obj.withGSParams(xvalue_accuracy=0.23, folding_threshold=7.3e-4)
    assert isinstance(obj_alt, obj.__class__)
    assert obj_alt.gsparams == alt_gsp
    assert obj_alt2.gsparams == alt_gsp
    print('obj_alt = ',obj_alt)
    print('obj.gsp = ',obj.gsparams)
    print('obj_alt.gsp = ',obj_alt.gsparams)
    print('gsp == ? ',obj.gsparams == obj_alt.gsparams)
    print('obj == ? ',obj == obj_alt)
    assert obj_alt != obj  # Assuming none of our tests use this exact gsparams choice.
    # Back to the original, ==, but not is
    assert obj_alt.withGSParams(obj.gsparams) is not obj
    assert obj_alt.withGSParams(obj.gsparams) == obj

    # Test errors for dimensionless SEDs.
    with assert_raises(galsim.GalSimSEDError):
        obj.drawImage(bps[0])
    with assert_raises(galsim.GalSimSEDError):
        obj.drawKImage(bps[0])
    with assert_raises(galsim.GalSimSEDError):
        obj.withFluxDensity(100., 500)
    with assert_raises(galsim.GalSimSEDError):
        obj.withFluxDensity(100., 500)
    with assert_raises(galsim.GalSimSEDError):
        obj.calculateFlux(bps[0])
    with assert_raises(galsim.GalSimSEDError):
        obj.calculateMagnitude(bps[0])


@timer
def test_chromatic_invariant():
    # Test atomic and non-transformed objects first.

    # GSObject
    flux = 2.0
    gsobj = galsim.Kolmogorov(fwhm=0.6, flux=flux)
    check_chromatic_invariant(gsobj)
    rng = galsim.BaseDeviate(1234)
    waves = rng.np.uniform(500, 1000, 10)
    np.testing.assert_allclose(gsobj.sed(waves), flux)

    # ChromaticObject
    check_chromatic_invariant(galsim.ChromaticObject(gsobj))
    np.testing.assert_allclose(galsim.ChromaticObject(gsobj).sed(waves), flux)
    chrom1 = galsim.ChromaticObject(gsobj) * bulge_SED
    chrom2 = gsobj * bulge_SED
    chrom3 = galsim.ChromaticObject(gsobj * bulge_SED)
    check_pickle(chrom1)
    check_pickle(chrom2)
    check_pickle(chrom3)
    check_pickle(galsim.ChromaticObject(gsobj))
    p1 = galsim.PhotonArray(1)
    p2 = galsim.PhotonArray(1)
    p2.wavelength = 500

    with assert_raises(TypeError):
        galsim.ChromaticObject(bulge_SED)

    check_chromatic_invariant(chrom1)
    check_chromatic_invariant(chrom2)
    check_chromatic_invariant(chrom3)
    # also check that these end up with the same SED
    np.testing.assert_allclose(chrom1.sed(waves), flux * bulge_SED(waves))
    np.testing.assert_allclose(chrom2.sed(waves), flux * bulge_SED(waves))
    np.testing.assert_allclose(chrom3.sed(waves), flux * bulge_SED(waves))
    assert chrom1.sed == chrom2.sed == chrom3.sed

    # And that they make the same image through a given bandpass
    bp = galsim.Bandpass('0.8 + wave/800*0.1', 'nm', blue_limit=700, red_limit=900)
    img = galsim.ImageF(64, 64, scale=0.2)
    img1 = chrom1.drawImage(bp, img)
    img2 = chrom2.drawImage(bp, img)
    img3 = chrom3.drawImage(bp, img)
    print('im1.sum = ',img1.array.sum())
    print('im2.sum = ',img2.array.sum())
    print('im3.sum = ',img3.array.sum())
    np.testing.assert_almost_equal(img1.array, img2.array, decimal=5)
    np.testing.assert_almost_equal(img1.array, img3.array, decimal=5)

    # ChromaticAtmosphere
    chrom_atm = galsim.ChromaticAtmosphere(gsobj, 500.0, zenith_angle=20.0 * galsim.degrees,
                                           pressure=70., temperature=285., H2O_pressure=1.05)
    check_chromatic_invariant(chrom_atm)
    check_pickle(chrom_atm)

    assert_raises(TypeError, galsim.ChromaticAtmosphere, gsobj,
                  500.0, zenith_angle=20.0 * galsim.degrees, invalid=3)

    # SimpleChromaticTransformation
    chrom = gsobj * bulge_SED
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves))
    check_pickle(chrom)

    with assert_raises(galsim.GalSimError):
        chrom.noise

    # Still Simple
    chrom = (gsobj * bulge_SED).shear(g1=0.2).shift((0.1,0.2))
    assert isinstance(chrom, galsim.SimpleChromaticTransformation)
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves))
    check_pickle(chrom)

    chrom = (gsobj * bulge_SED) * (lambda w: w**1.03)
    assert isinstance(chrom, galsim.SimpleChromaticTransformation)
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves) * waves**1.03)
    # Not picklable, but run str, repr
    str(chrom)
    repr(chrom)

    # ChromaticTransformation
    chrom = galsim.Transform(gsobj, offset=(0.1, 0.3), flux_ratio=bulge_SED)
    assert isinstance(chrom, galsim.ChromaticTransformation)
    assert not isinstance(chrom, galsim.SimpleChromaticTransformation)
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves))
    check_pickle(chrom)

    with assert_raises(galsim.GalSimError):
        chrom.noise

    chrom = galsim.Transform(gsobj, jac=galsim.Shear(g1=0.1).getMatrix(), flux_ratio=bulge_SED)
    assert isinstance(chrom, galsim.ChromaticTransformation)
    assert not isinstance(chrom, galsim.SimpleChromaticTransformation)
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves))
    check_pickle(chrom)

    scaling = galsim.SED(lambda w: w**1.03, 'nm', '1')
    chrom = galsim.Transform(gsobj * bulge_SED, offset=(0.1, 0.3), flux_ratio=scaling)
    assert isinstance(chrom, galsim.ChromaticTransformation)
    assert not isinstance(chrom, galsim.SimpleChromaticTransformation)
    check_chromatic_invariant(chrom)
    np.testing.assert_allclose(chrom.sed(waves), flux * bulge_SED(waves) * waves**1.03, rtol=1.e-4)
    # Not picklable, but run str, repr
    str(chrom)
    repr(chrom)

    chrom = galsim.Transform(gsobj, jac = (0.5, 0.5, 1, -1), flux_ratio=bulge_SED)
    expected_flux = flux*bulge_SED.calculateFlux(bandpass)
    chromobj_flux = chrom.calculateFlux(bandpass)
    np.testing.assert_allclose(expected_flux, chromobj_flux)

    # ChromaticOpticalPSF
    chrom_opt = galsim.ChromaticOpticalPSF(lam=500.0, diam=2.0, tip=2.0, tilt=3.0, defocus=0.2,
                                           scale_unit='arcmin')
    check_chromatic_invariant(chrom_opt)
    np.testing.assert_allclose(chrom_opt.sed(waves), 1.0)
    check_pickle(chrom_opt)

    # ChromaticAiry
    chrom_airy = galsim.ChromaticAiry(lam=500.0, diam=3.0, scale_unit=galsim.arcmin)
    check_chromatic_invariant(chrom_airy)
    np.testing.assert_allclose(chrom_airy.sed(waves), 1.0)
    check_pickle(chrom_airy)

    # Now start testing compound objects...
    # ChromaticSum
    chrom_sum_noSED = chrom_airy + chrom_opt
    check_chromatic_invariant(chrom_sum_noSED)
    np.testing.assert_allclose(chrom_sum_noSED.sed(waves), 2.0)
    # TODO: Seems like this should be picklable. Probably anything that doesn't include
    #       unpicklable user input should be picklable.
    #       e.g. autoconv2 has no hope.  But there are a few check_pickle calls that are commented
    #       out that we should probably try to make work.  A job for another day, though...
    #check_pickle(chrom_sum_noSED)
    repr(chrom_sum_noSED)
    str(chrom_sum_noSED)
    assert_raises(galsim.GalSimError, chrom_sum_noSED.applyTo, p1)

    chrom = gsobj * bulge_SED
    chrom_sum_SED = chrom + chrom  # used to be considered separable, but not anymore.
    check_chromatic_invariant(chrom_sum_SED)
    np.testing.assert_allclose(chrom_sum_SED.sed(waves), 2*flux*bulge_SED(waves))
    check_pickle(chrom_sum_SED)
    assert not chrom_sum_SED.separable

    gsobj2 = galsim.Kolmogorov(fwhm=0.7)
    chrom2 = gsobj2 * disk_SED
    chrom_sum_SED2 = chrom + chrom2
    check_chromatic_invariant(chrom_sum_SED2)
    np.testing.assert_allclose(chrom_sum_SED2.sed(waves), flux*bulge_SED(waves) + disk_SED(waves))
    check_pickle(chrom_sum_SED2)
    assert not chrom_sum_SED2.separable
    assert_raises(galsim.GalSimError, chrom_sum_SED.applyTo, p1)

    # ChromaticConvolution
    conv1 = galsim.Convolve(chrom, chrom_airy)  # SEDed
    check_chromatic_invariant(conv1)
    np.testing.assert_allclose(conv1.sed(waves), flux*bulge_SED(waves))
    check_pickle(conv1)

    with assert_raises(galsim.GalSimError):
        conv1.noise

    conv2 = galsim.Convolve(chrom_airy, chrom_opt)  # Non-SEDed
    check_chromatic_invariant(conv2)
    np.testing.assert_allclose(conv2.sed(waves), 1.)
    check_pickle(conv2)
    assert_raises(galsim.GalSimError, conv1.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, conv1.applyTo, p2)
    assert_raises(galsim.GalSimError, conv2.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, conv2.applyTo, p2)

    # ChromaticDeconvolution
    deconv = galsim.Deconvolve(chrom_airy)
    check_chromatic_invariant(deconv)
    np.testing.assert_allclose(deconv.sed(waves), 1.0)
    #check_pickle(deconv)
    repr(deconv) # gratuitous coverage of repr until check_pickle works.
    str(deconv)
    assert_raises(galsim.GalSimError, deconv.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, deconv.applyTo, p2)
    # Repeat with non-unit SED; deconv SED should be inverse of base SED.
    deconv2 = galsim.Deconvolve(chrom_airy * (lambda w:(w/500)**3))
    np.testing.assert_allclose(deconv2.sed(waves), (waves/500)**-3)

    # ChromaticAutoConvolution
    autoconv1 = galsim.AutoConvolve(chrom_airy)
    check_chromatic_invariant(autoconv1)
    np.testing.assert_allclose(autoconv1.sed(waves), 1.0)
    autoconv2 = galsim.AutoConvolve(chrom_airy * (lambda w: (w/500.0)**0.1))
    check_chromatic_invariant(autoconv2)
    np.testing.assert_allclose(autoconv2.sed(waves), (waves/500)**0.2)
    check_pickle(autoconv1)
    assert_raises(galsim.GalSimError, autoconv1.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, autoconv1.applyTo, p2)

    # ChromaticAutoCorrelation
    autocorr1 = galsim.AutoCorrelate(chrom_airy)
    check_chromatic_invariant(autocorr1)
    np.testing.assert_allclose(autocorr1.sed(waves), 1.0)
    autocorr2 = galsim.AutoCorrelate(chrom_airy * (lambda w: (w/500.0)**0.1))
    check_chromatic_invariant(autocorr2)
    np.testing.assert_allclose(autocorr2.sed(waves), (waves/500)**0.2)
    check_pickle(autocorr1)
    assert_raises(galsim.GalSimError, autocorr1.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, autocorr1.applyTo, p2)

    # ChromaticFourierSqrt
    four1 = galsim.FourierSqrt(chrom_airy)
    check_chromatic_invariant(four1)
    np.testing.assert_allclose(four1.sed(waves), 1.0)
    four2 = galsim.FourierSqrt(chrom_airy * (lambda w: (w/500.0)**0.1))
    check_chromatic_invariant(four2)
    np.testing.assert_allclose(four2.sed(waves), (waves/500)**0.05)
    #check_pickle(four1)
    repr(four1) # gratuitous coverage of repr until check_pickle works.
    str(four1)
    assert_raises(galsim.GalSimError, four1.applyTo, p1)
    assert_raises(galsim.GalSimNotImplementedError, four1.applyTo, p2)

    # And a few transforms too...
    # ChromaticTransformation
    sheared_chrom = chrom.shear(g1=0.1)
    check_chromatic_invariant(sheared_chrom)
    np.testing.assert_allclose(sheared_chrom.sed(waves), chrom.sed(waves))
    check_pickle(sheared_chrom)

    scaled_chrom = 2 * chrom
    check_chromatic_invariant(scaled_chrom)
    np.testing.assert_allclose(scaled_chrom.sed(waves), 2*chrom.sed(waves))
    check_pickle(scaled_chrom)

    complex_scaled_chrom = chrom * (lambda w: (w/500.0)**0.1)
    np.testing.assert_allclose(complex_scaled_chrom.sed(waves), chrom.sed(waves) * (waves/500)**0.1)
    check_chromatic_invariant(complex_scaled_chrom)

    chrom_added_SED = chrom_airy * bulge_SED
    check_chromatic_invariant(chrom_added_SED)
    np.testing.assert_allclose(chrom_added_SED.sed(waves), bulge_SED(waves))
    check_pickle(chrom_added_SED)

    complex_expanded_chrom = chrom.expand(lambda w: (w/500.0)**0.1)
    check_chromatic_invariant(complex_expanded_chrom)
    np.testing.assert_allclose(complex_expanded_chrom.sed(waves),
                               chrom.sed(waves) * (waves/500)**0.2)

    complex_dilated_chrom = chrom.dilate(lambda w: (w/500.0)**0.1)
    check_chromatic_invariant(complex_dilated_chrom)
    np.testing.assert_allclose(complex_dilated_chrom.sed(waves), chrom.sed(waves))

    complex_transformed_chrom = chrom.transform(
            lambda w: 1. + 0.2*(w/500)**0.1,
            lambda w: -0.2*(w/500)**0.2,
            lambda w: -0.3*(w/500)**0.3,
            lambda w: 1. + 0.1*(w/500)**0.2
    )
    check_chromatic_invariant(complex_transformed_chrom)
    jac = ((1+0.2*(waves/500)**0.1) * (1+0.1*(waves/500)**0.2) -
           (-0.2*(waves/500)**0.2) * (-0.3*(waves/500)**0.3))
    np.testing.assert_allclose(complex_transformed_chrom.sed(waves), chrom.sed(waves) * jac)

    # ChromaticInterpolatedObject
    chrom_interp = chrom_airy.interpolate(waves=[400.0, 500.0, 600.0])
    check_chromatic_invariant(chrom_interp)
    np.testing.assert_allclose(chrom_interp.sed(waves), 1.0)
    check_pickle(chrom_interp)
    assert_raises(galsim.GalSimError, chrom_interp.applyTo, p1)


@timer
def test_ne():
    """Test chromatic.py objects for not-equals."""
    # Define some universal objects.
    gsp = galsim.GSParams(maxk_threshold=1.1e-3, folding_threshold=5.1e-3)
    gal1 = galsim.Gaussian(fwhm=1.0)
    gal2 = galsim.Gaussian(fwhm=1.1)
    cgal1 = galsim.ChromaticObject(gal1).dilate(lambda w:1)
    cgal2 = galsim.ChromaticObject(gal2).dilate(lambda w:1)
    cgal3 = cgal1.interpolate(np.arange(400, 550, 10))

    # ChromaticObject.  Only param is the GSObject to chromaticize.
    # The following should test unequal:
    gals = [cgal1, cgal2]
    check_all_diff(gals)

    # # Check that setifying doesn't remove any duplicate items.
    # assert len(gals) == len(set(gals))
    # hashes = [hash(item) for item in gals]
    # assert len(hashes) == len(set(hashes))

    # InterpolatedChromaticObject.  Params include object to be chromaticized, waves, and
    # oversample_fac.
    # Also get a copy of cgal1 and make it interpolatable, but with a different waves argument.
    gals = [cgal1,
            cgal2,
            cgal3,
            galsim.InterpolatedChromaticObject(cgal1, np.arange(500, 700, 50)),
            galsim.InterpolatedChromaticObject(cgal2, np.arange(500, 700, 50)),
            galsim.InterpolatedChromaticObject(cgal1, np.arange(500, 700, 25)),
            galsim.InterpolatedChromaticObject(cgal1, np.arange(500, 700, 50), oversample_fac=1.1),
            cgal1.interpolate(np.arange(500, 700, 10)),
            cgal1.interpolate(np.arange(500, 700, 10), oversample_fac=1.2)]
    check_all_diff(gals)

    # ChromaticAtmosphere.  Params include base_obj, base_wavelength, scale_unit, alpha,
    # zenith_angle, parallactic_angle, obj_coord, zenith_coord, HA, latitude, pressure, temperature,
    # and H2O_pressure.

    # Some params to test the celestial coords inputs.
    ra = galsim.Angle.from_hms("14:03:13") # hours : minutes : seconds
    dec = galsim.Angle.from_dms("54:20:57") # degrees : minutes : seconds
    m101 = galsim.CelestialCoord(ra, dec)
    latitude = 19.8207 * galsim.degrees # latitude of Mauna Kea
    HA = -1.0 * galsim.hours # Hour angle = one hour before transit
    zenith_coord = galsim.CelestialCoord(0.*galsim.degrees, 40*galsim.degrees) # Make something up
    # The following should test unequal.
    gals = [galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees),
            galsim.ChromaticAtmosphere(gal2, 500.0, zenith_angle=30*galsim.degrees),
            galsim.ChromaticAtmosphere(gal1, 600.0, zenith_angle=30*galsim.degrees),
            galsim.ChromaticAtmosphere(gal1, 600.0, zenith_angle=30*galsim.degrees,
                                       scale_unit=galsim.degrees),
            galsim.ChromaticAtmosphere(gal1, 600.0, zenith_angle=30*galsim.degrees,
                                       scale_unit='arcmin'),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees, alpha=-0.1),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees,
                                       parallactic_angle=45*galsim.degrees),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees,
                                       temperature=290.0),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees,
                                       pressure=75.0),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees,
                                       H2O_pressure=2.0),
            galsim.ChromaticAtmosphere(gal1, 500.0, obj_coord=m101, HA=HA, latitude=latitude),
            galsim.ChromaticAtmosphere(gal1, 500.0, obj_coord=m101, zenith_coord=zenith_coord)]
    check_all_diff(gals)

    # Chromatic.  Params are a gsobject and an SED.
    sed1 = galsim.SED(lambda w: w, 'nm', 'flambda')
    sed2 = galsim.SED(lambda w: 2*w, 'nm', 'flambda')
    # The following should test unequal.
    gals = [gal1 * sed1,
            gal1 * sed2,
            gal2 * sed1,
            gal2 * sed2]
    check_all_diff(gals)

    # ChromaticTransformation.  Params are an object (possibly chromatic), a jacobian jac, an
    # offset, a flux_ratio, and gsparams.  For coverage, test jac, offset, and flux_ratio as
    # consts and functions.
    jac1 = lambda w: [[w, 0], [0, 1]]
    jac2 = lambda w: [[w, 0], [0, w]]
    offset1 = lambda w: (0, w)
    offset2 = lambda w: (w, 0)
    flux_ratio1 = lambda w: w
    flux_ratio2 = lambda w: w**2
    # The following should test unequal.
    with assert_warns(galsim.GalSimWarning):
        trans_cgal3 = galsim.ChromaticTransformation(
                cgal3, jac=jac1, offset=offset1, flux_ratio=flux_ratio1),
    gals = [galsim.ChromaticTransformation(cgal1),
            galsim.ChromaticTransformation(cgal3),
            galsim.ChromaticTransformation(gal1, jac=[[1, 1.1], [0.1, 1]]),
            galsim.ChromaticTransformation(gal1, jac=[[1, 0.1], [0.1, 1]]),
            galsim.ChromaticTransformation(gal1, jac=jac1),
            galsim.ChromaticTransformation(gal1, jac=jac2),
            galsim.ChromaticTransformation(gal1, offset=(0.1, 0.0)),
            galsim.ChromaticTransformation(gal1, offset=(0.0, 0.1)),
            galsim.ChromaticTransformation(gal1, offset=offset1),
            galsim.ChromaticTransformation(gal1, offset=offset2),
            galsim.ChromaticTransformation(gal1, flux_ratio=1.1),
            galsim.ChromaticTransformation(gal1, flux_ratio=1.4),
            galsim.ChromaticTransformation(gal1, flux_ratio=flux_ratio1),
            galsim.ChromaticTransformation(gal1, flux_ratio=flux_ratio2),
            galsim.ChromaticTransformation(cgal1, jac=jac1, offset=offset1, flux_ratio=flux_ratio1),
            trans_cgal3,
            galsim.ChromaticTransformation(cgal1, gsparams=gsp),
            galsim.ChromaticTransformation(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticTransformation(cgal3).withGSParams(gsp),
            galsim.ChromaticTransformation(cgal3, propagate_gsparams=False).withGSParams(gsp),
            galsim.SimpleChromaticTransformation(gal1),
            galsim.SimpleChromaticTransformation(gal1, gsparams=gsp),
            galsim.SimpleChromaticTransformation(gal1, gsparams=gsp, propagate_gsparams=False)]
    check_all_diff(gals)

    # ChromaticSum.  Params are objs to add and potentially gsparams.
    # The following should test unequal.
    gals = [galsim.ChromaticSum(cgal1),
            galsim.ChromaticSum(cgal1, cgal2),
            galsim.ChromaticSum(cgal3, cgal2),
            galsim.ChromaticSum(cgal2, cgal1),  # Not! commutative.
            galsim.ChromaticSum(galsim.ChromaticSum(cgal1, cgal2), cgal2),
            galsim.ChromaticSum(cgal1, galsim.ChromaticSum(cgal2, cgal2)),  # Not! associative.
            galsim.ChromaticSum(cgal1, gsparams=gsp),
            galsim.ChromaticSum(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticSum(cgal3).withGSParams(gsp),
            galsim.ChromaticSum(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticConvolution.  Params are objs to convolve and potentially gsparams.
    # The following should test unequal
    with assert_warns(galsim.GalSimWarning):
        conv_32 = galsim.ChromaticConvolution(cgal3, cgal2),
    gals = [galsim.ChromaticConvolution(cgal1),
            galsim.ChromaticConvolution(cgal1, cgal2),
            conv_32,
            galsim.ChromaticConvolution(cgal2, cgal1),  # Not! commutative.
            galsim.ChromaticConvolution(galsim.ChromaticConvolution(cgal1, cgal2), cgal2),
            # ChromaticConvolution is associative! (unlike galsim.Convolution)
            # galsim.ChromaticConvolution(cgal1, galsim.ChromaticConvolution(cgal2, cgal2)),
            galsim.ChromaticConvolution(cgal1, gsparams=gsp),
            galsim.ChromaticConvolution(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticConvolution(cgal3).withGSParams(gsp),
            galsim.ChromaticConvolution(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticDeconvolution.  Only params here are obj to deconvolve and gsparams.
    gals = [galsim.ChromaticDeconvolution(cgal1),
            galsim.ChromaticDeconvolution(cgal2),
            galsim.ChromaticDeconvolution(cgal3),
            galsim.ChromaticDeconvolution(cgal1, gsparams=gsp),
            galsim.ChromaticDeconvolution(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticDeconvolution(cgal3).withGSParams(gsp),
            galsim.ChromaticDeconvolution(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticAutoConvolution.
    gals = [galsim.ChromaticAutoConvolution(cgal1),
            galsim.ChromaticAutoConvolution(cgal2),
            galsim.ChromaticAutoConvolution(cgal3),
            galsim.ChromaticAutoConvolution(cgal1, gsparams=gsp),
            galsim.ChromaticAutoConvolution(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticAutoConvolution(cgal3).withGSParams(gsp),
            galsim.ChromaticAutoConvolution(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticAutoCorrelation.
    gals = [galsim.ChromaticAutoCorrelation(cgal1),
            galsim.ChromaticAutoCorrelation(cgal2),
            galsim.ChromaticAutoCorrelation(cgal3),
            galsim.ChromaticAutoCorrelation(cgal1, gsparams=gsp),
            galsim.ChromaticAutoCorrelation(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticAutoCorrelation(cgal3).withGSParams(gsp),
            galsim.ChromaticAutoCorrelation(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticFourierSqrt.
    gals = [galsim.ChromaticFourierSqrtProfile(cgal1),
            galsim.ChromaticFourierSqrtProfile(cgal2),
            galsim.ChromaticFourierSqrtProfile(cgal3),
            galsim.ChromaticFourierSqrtProfile(cgal1, gsparams=gsp),
            galsim.ChromaticFourierSqrtProfile(cgal1, gsparams=gsp, propagate_gsparams=False),
            galsim.ChromaticFourierSqrtProfile(cgal3).withGSParams(gsp),
            galsim.ChromaticFourierSqrtProfile(cgal3, propagate_gsparams=False).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticOpticalPSF.  Params include: lam, (diam or lam_over_diam), aberrations, nstruts,
    # strut_thick, strut_angle, obscuration, oversampling, pad_factor, flux, gsparams, ...
    # Most of these get tested in the same way, (via a kwargs dict comparison), so only test a few
    # here.
    gals = [galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0),
            galsim.ChromaticOpticalPSF(lam=1.0, diam=1.1),
            galsim.ChromaticOpticalPSF(lam=1.0, diam=1.1, fft_sign='-'),
            galsim.ChromaticOpticalPSF(lam=1.0, diam=1.1, scale_unit=galsim.arcmin),
            galsim.ChromaticOpticalPSF(lam=1.0, diam=1.1, scale_unit='radians'),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, obscuration=0.5),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, aberrations=[0, 0, 0, 0, 0.1]),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, defocus=0.2),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, flux=0.2),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, gsparams=gsp),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0, flux=0.2).withGSParams(gsp)]
    check_all_diff(gals)

    # ChromaticAiry.  Params include: lam, diam, lam_over_diam, scale_unit, flux, obscuration,
    # gsparams.
    gals = [galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0),
            galsim.ChromaticAiry(lam=1.0, diam=1.0),
            galsim.ChromaticAiry(lam=1.0, diam=1.0, scale_unit=galsim.arcmin),
            galsim.ChromaticAiry(lam=1.0, diam=1.0, scale_unit='deg'),
            galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0, obscuration=0.5),
            galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0, flux=1.1),
            galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0, gsparams=gsp),
            galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0, flux=1.1).withGSParams(gsp)]
    check_all_diff(gals)

    # Check that all the various combinations are properly unequal
    gals = [cgal1,
            galsim.ChromaticObject(gal1),
            galsim.InterpolatedChromaticObject(cgal1, np.arange(500, 700, 50)),
            galsim.ChromaticAtmosphere(gal1, 500.0, zenith_angle=30*galsim.degrees),
            gal1 * sed1,
            galsim.ChromaticTransformation(cgal1),
            galsim.ChromaticSum(cgal1),
            galsim.ChromaticConvolution(cgal1),
            galsim.ChromaticDeconvolution(cgal1),
            galsim.ChromaticAutoConvolution(cgal1),
            galsim.ChromaticAutoCorrelation(cgal1),
            galsim.ChromaticOpticalPSF(lam=1.0, lam_over_diam=1.0),
            galsim.ChromaticAiry(lam=1.0, lam_over_diam=1.0)]
    check_all_diff(gals)

@timer
def test_shoot_transformation():
    """Check that transformed chromatic objects can be photon shot.
    """
    # In response to #1229
    bandpass = galsim.Bandpass('LSST_r.dat', 'nm')
    rng = galsim.BaseDeviate(1234)
    flux = 1000

    # Dilate
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=0.5)
    ).dilate(lambda w: (w/500)**-0.3)
    star = galsim.DeltaFunction()*galsim.SED('vega.txt', 'nm', 'flambda')
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    # The real test is that this didn't error.  So this assert isn't very exciting.
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux)

    # We used to do the wrong thing with poisson_flux=True.  Check that the flux isn't
    # just 1000 in that case.
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Rotate
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).rotate(lambda w: (w/500)*8. * galsim.degrees)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Expand
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).expand(lambda w: (w/500)**-0.3)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    # Transformations with a non-unit det get the flux right in an expecations sense, but
    # because the flux_ratios vary with wavelength, there is some noise on this match.
    # With only 1000 photons, it only matches to better than 2.e-3.
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Shift
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).shift(lambda w: ((w/500)*0.3, (w/500)*-0.4))
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Shear
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).shear(lambda w: galsim.Shear(g1=(w/500)*0.03, g2=(w/500)*0.05))
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Magnify
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).magnify(lambda w: (w/500)**-0.3)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Lens
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).lens(g1=lambda w: (w/500)*0.03, g2=lambda w: (w/500)*0.05, mu=lambda w: (w/500)**-0.3)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Transform
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).transform(dudx=lambda w: (w/500)**0.02, dudy=lambda w: (w/500)*0.11,
                dvdx=lambda w: (w/500)*0.08, dvdy=lambda w: (w/500)**0.03)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

    # Flux_scale
    psf = galsim.ChromaticObject(
        galsim.Gaussian(fwhm=1)
    ).withScaledFlux(lambda w: (w/500)**1.3)
    obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                        poisson_flux=False)
    print(img.added_flux)
    np.testing.assert_allclose(img.added_flux, flux, rtol=1.e-6)
    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
    print(img.added_flux)
    assert abs(img.added_flux - flux) > 0.1

@timer
def test_save_photons():
    # Test save_photons=True option for various Chromatic types
    # Note: there are additional tests of save_photons in test_phot, but that doesn't cover all
    # possible types and all possible paths through the various drawImage and shoot functions.
    # This test tries to be more comprehensive on that front.

    rng = galsim.BaseDeviate(1234)
    star_sed = galsim.SED('vega.txt', wave_type="nm", flux_type="fphotons")
    bandpass = galsim.Bandpass("LSST_r.dat", wave_type="nm")

    airy = galsim.ChromaticAiry(lam=500, diam=8)
    optical = galsim.ChromaticOpticalPSF(lam=500, diam=8, defocus=0.2, obscuration=0.3)
    disk = galsim.Exponential(half_light_radius=0.5).shear(g1=0.4, g2=0.2)
    bulge = galsim.Sersic(n=3, half_light_radius=0.3)
    atm = galsim.ChromaticAtmosphere(galsim.Kolmogorov(fwhm=0.9), base_wavelength=500.0,
                                     zenith_angle=20 * galsim.degrees)

    objs = [
        airy * star_sed,
        optical * star_sed,
        (airy * star_sed).expand(lambda w: (w/500)**0.5),
        disk * disk_SED.atRedshift(1.1),
        (bulge * bulge_SED + disk * disk_SED),
        galsim.Convolve(disk * disk_SED, optical, atm),
        galsim.Convolve(disk * disk_SED, atm.interpolate(np.linspace(500,900,5))),
        (atm * star_sed).interpolate(np.linspace(500,900,5)),
    ]

    flux = 1000
    for obj in objs:
        print('obj = ',obj)
        obj = obj.withFlux(flux, bandpass)
        image = obj.drawImage(bandpass=bandpass, method="phot",
                              n_photons=flux, save_photons=True,
                              scale=0.05, nx=32, ny=32, rng=rng)
        assert hasattr(image, 'photons')
        assert len(image.photons) == flux
        print(np.sum(image.photons.flux))
        # Note: tolerance is quite loose, since profiles that use InterpolatedImage can have
        # negative flux photons, which then don't necessarily sum to the right value.
        # Only the expectation value is right, and we're not shooting many photons here.
        assert np.allclose(np.sum(image.photons.flux), flux, rtol=0.1)

        # Sometimes there is a different path when n_photons is not given, so check that too.
        image = obj.drawImage(bandpass=bandpass, method="phot",
                              save_photons=True,
                              scale=0.05, nx=32, ny=32, rng=rng)
        assert hasattr(image, 'photons')
        print(np.sum(image.photons.flux))
        assert np.allclose(np.sum(image.photons.flux), flux, rtol=0.1)


if __name__ == "__main__":
    runtests(__file__)