def METHOD_NAME():
    # compare the differences in chromatic moment shifts between two SEDs as a function of
    # Bandpass thinning.  Goals should be to keep the error below:
    # sigma(dRbar) < 0.01 arcsec
    # sigma(dV) < 0.0001 arcsec^2
    # sigma(dseeing) < 0.0001
    import glob
    SED_files = glob.glob(os.path.join(sedpath, '*.sed'))
    bp_files = glob.glob(os.path.join(bppath, '*.dat'))
    SED1 = galsim.SED(os.path.join(sedpath, 'CWW_E_ext.sed'),
                      wave_type='nm', flux_type='flambda').withFluxDensity(0.01, 500.0)
    SEDs = dict([(os.path.basename(SED_file),
                  galsim.SED(SED_file, wave_type='nm', flux_type='flambda')) for SED_file in SED_files])
    del SEDs['CWW_E_ext.sed']
    bands = dict([(os.path.basename(bp_file),
                   galsim.Bandpass(bp_file, wave_type='nm')) for bp_file in bp_files])
    redshifts = [0.0, 0.5, 1.0]
    rel_errs = [1.e-4, 1.e-3]
    for SED_name, SED0 in SEDs.items():
        for redshift in redshifts:
            SED = SED0.atRedshift(redshift).withFluxDensity(0.01, 500.0)
            for bandname, band in bands.items():
                print('{0} SED at z={1} through {2} filter'.format(
                    SED_name, redshift, bandname))
                dDCR = dDCR_moments(SED1, SED, band)
                dseeing = dseeing_moments(SED1, SED, band)
                flux = SED.calculateFlux(band)
                hdr = '{0:8s} {1:>8s} {2:>8s} {3:>8s} {4:>8s} {5:>8s} {6:>8s} {7:>8s} {8:>8s}'
                print(hdr.format(
                    'rel_err', 'dRbar', 'dV', 'dseeing', 'flux',
                    'd(dRbar)', 'd(dV)', 'd(dseeing)', 'd(flux)/flux'))
                out = '{0:8} {1:8.5f} {2:8.5f} {3:8.5f} {4:8.5f}'
                print(out.format('full', dDCR[0], dDCR[1], dseeing, flux))
                for rel_err in rel_errs:
                    band1 = band.thin(rel_err=rel_err)
                    dDCR_thinned = dDCR_moments(SED1, SED, band1)
                    dseeing_thinned = dseeing_moments(SED1, SED, band1)
                    flux_thinned = SED.calculateFlux(band1)
                    out = ('{0:8s} {1:8.5f} {2:8.5f} {3:8.5f} {4:8.5f}'
                           +' {5:8.5f} {6:8.5f} {7:8.5f} {8:8.5f}')
                    print(out.format(
                        str(rel_err), dDCR_thinned[0], dDCR_thinned[1],
                        dseeing_thinned, flux_thinned,
                        dDCR_thinned[0] - dDCR[0], dDCR_thinned[1] - dDCR[1],
                        dseeing_thinned - dseeing, (flux_thinned - flux)/flux))
                print()
    print('{0:8s} {1:>8s} {2:>8s}'.format('rel_err', 'time', 'Neval'))
    t0 = time.time()
    for i in range(20):
        dDCR_thinned = dDCR_moments(SED1, SED, band)
        dseeing_thinned = dseeing_moments(SED1, SED, band)
        flux_thinned = SED.calculateFlux(band)
    t1 = time.time()
    print('{0:8s} {1:8.5f} {2:8d}'.format('full', t1-t0, dDCR_thinned[2]))
    for rel_err in rel_errs:
        band1 = band.thin(rel_err=rel_err)
        t0 = time.time()
        for i in range(20):
            dDCR_thinned = dDCR_moments(SED1, SED, band1)
            dseeing_thinned = dseeing_moments(SED1, SED, band1)
            flux_thinned = SED.calculateFlux(band1)
        t1 = time.time()
        print('{0:8s} {1:8.5f} {2:8d}'.format(str(rel_err), t1-t0, dDCR_thinned[2]))