deredden V2.10.3 : weird behaviour at ~2000-3000A ?

hachi1 wrote on Jan 12, 2011

Dear all,

I programmed an own de-reddening routine in FORTRAN (CCM extinction law - just as in deredden) and compared the output in the end (flat spectrum with flux=1.0, deredden by E(B-V)=0.5 w/ R_V=3.1).

Plot:

http://www.mpa-garching.mpg.de/~shaching/extinction_problem.png

Now the weird thing is that IRAF and my routine are somewhat different at ~2000...3000A, and the IRAF routine in fact produces a little jump in extinction at ~3050A.

Now... is that a bug of my routine or of IRAF? Or is it a feature ... maybe it has to do with air-vacuum wl conversion?

Thanks in advance for any suggestions!

S. Hachinger

---

Francisco Valdes wrote on Jan 12, 2011

The deredden code is not very complex. There is no vacuum correction applied though you are right that one must be careful to use wavelengths as intended by the authors in the reference paper. I attach the IRAF code that computes extinction law. While this is in SPP you should have no problem understanding the language. Maybe you will spot something that you are doing differently.

Frank Valdes

# CCM -- Compute CCM Extinction Law

real procedure ccm (wavelength, rv)

real    wavelength              # Wavelength in Angstroms
real    rv                      # A(V) / E(B-V)

real    x, y, a, b

begin
        # Convert to inverse microns
        x = 10000. / wavelength

        # Compute a(x) and b(x)
        if (x < 0.3) {
            call error (1, "Wavelength out of range of extinction function")

        } else if (x < 1.1) {
            y = x ** 1.61
            a = 0.574 * y
            b = -0.527 * y

        } else if (x < 3.3) {
            y = x - 1.82
            a = 1 + y * (0.17699 + y * (-0.50447 + y * (-0.02427 +
                y * (0.72085 + y * (0.01979 + y * (-0.77530 + y * 0.32999))))))
            b = y * (1.41338 + y * (2.28305 + y * (1.07233 + y * (-5.38434 +
                y * (-0.62251 + y * (5.30260 + y * (-2.09002)))))))

        } else if (x < 5.9) {
            y = (x - 4.67) ** 2
            a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
            b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)

        } else if (x < 8.0) {
            y = (x - 4.67) ** 2
            a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
            b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)

            y = x - 5.9
            a = a - 0.04473 * y**2 - 0.009779 * y**3
            b = b + 0.2130 * y**2 + 0.1207 * y**3

        } else if (x <= 10.0) {
            y = x - 8
            a = -1.072 - 0.628 * y + 0.137 * y**2 - 0.070 * y**3
            b = 13.670 + 4.257 * y - 0.420 * y**2 + 0.374 * y**3

        } else {
            call error (1, "Wavelength out of range of extinction function")

        }

        # Compute A(lambda)/A(V)
        y = a + b / rv
        return (y)
end

---

hachi1 wrote on Jan 12, 2011

Dear Frank,

sorry for my late reply. I have figured out the reason for the difference: in fact there seems to be a bug in the IRAF deredden code in the section:

} else if (x < 8.0) {
y = (x - 4.67) ** 2
a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)

y = x - 5.9
a = a - 0.04473 * y**2 - 0.009779 * y**3
b = b + 0.2130 * y**2 + 0.1207 * y**3

}

According to Cardelli, Clayton & Mathis 1989 (p. 35), the code would have to be sth like:

} else if (x < 8.0) {
y = (x - 4.67) ** 2
a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
y = (x - 4.62) ** 2
b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)

y = x - 5.9
a = a - 0.04473 * y**2 - 0.009779 * y**3
b = b + 0.2130 * y**2 + 0.1207 * y**3

}

Please check this -- if I'm right: who takes care that such things are corrected, i.e. where can I submit a bug report?

Cheers

Stephan

---

Francisco Valdes wrote on Jan 12, 2011

Stephan,

You are correct that the IRAF routine did not correctly code the CCM formulae. Your post is sufficient to alert us to the error and I, the original author of the IRAF task, will fix it for the next release.

Thank you,
Frank Valdes

---

hachi1 wrote on Jan 12, 2011

Hi Frank,

ok, thanks very much!

Cheers

Stephan

---

Last post on Jan 12, 2011