Normalise the attenuation curves so that their E(B-V) corresponds to the...

Normalise the attenuation curves so that their E(B-V) corresponds to the E(B-V) requested by the user.

CHANGELOG.md

@@ -8,6 +8,7 @@
 - The estimates of the physical parameters from the analysis of the PDF and from the best fit were recorded in separate files. This can be bothersome when trying to compare quantities from different files. Rather, we generate a single file containing all quantities. The ones estimated from the analysis of the PDF are prefixed with "bayes" and the ones from the best fit with "best". (Médéric Boquien)
 - To homogenize input and output files, the "observation_id" has been changed to "id" in the output files. (Médéric Boquien)
 - The output files providing estimates of the physical properties are now generated both in the form of text and FITS files. (Médéric Boquien)
+- When using the dustatt_calzleit module, choosing ẟ≠0 leads to an effective E(B-V) different from the one set by the user. Now the E(B-V) will always correspond to the one specified by the user. This means that at fixed E(B-V), A(V) depends on ẟ. (Médéric Boquien)
 
 ### Fixed
 - To estimate parameters in log, pcigale determines which variables end with the "_log" string and removed it to compute the models. However in some circumstances, it was overzealous. This has been fixed. (Médéric Boquien)

pcigale/creation_modules/dustatt_calzleit.py

@@ -166,6 +166,16 @@ def a_vs_ebv(wavelength, bump_wave, bump_width, bump_ampl, power_slope):
     # UV bump
     attenuation += uv_bump(wavelength, bump_wave, bump_width, bump_ampl)
 
+    # As the powerlaw slope changes E(B-V), we correct this so that the curve
+    # always has the same E(B-V) as the starburst curve. This ensures that the
+    # E(B-V) requested by the user is the actual E(B-V) of the curve.
+    wl_BV = np.array([440., 550.])
+    EBV_calz = ((k_calzetti2000(wl_BV) * power_law(wl_BV, 0.)) +
+                uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
+    EBV = ((k_calzetti2000(wl_BV) * power_law(wl_BV, power_slope)) +
+           uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
+    attenuation *= (EBV_calz[1]-EBV_calz[0]) / (EBV[1]-EBV[0])
+
     return attenuation