## Bug #726

open### Incorrect management of internal dust extinction in disks?

70%

**Description**

(!) This is an issue previously reported by Damien Leborgne, which I erased by mistake...

In stuff version 1.26

---------------------

Unwanted or strange behaviour:

if a population is defined with `B/T (BULGE_FRACTION) < 1`

and `DISK_EXTINCT = 0`

, the absolute (and therefore apparent) magnitudes are too bright by a factor $-2.5 \log_{10} \left(1 + (1-B/T)*10^{0.12041\,\rm{extinct}}\right)$ with respect to the LF defined by the Schechter parameters.

In other words, I would have thought that the 2 following configs should lead to the same LF (the one defined by LF_PHISTAR, LF_MSTAR, and LF_ALPHA), or equivalently to the same apparent magnitude distribution, but it is not the case:

**Files**

#### Updated by Damien Le Borgne about 5 years ago

stuff replace $10^{0.12041\,\rm{extinct}}$ by $1.0-10^{-0.12041\, \rm{extinct}}$ seems to work.

#### Updated by Emmanuel Bertin about 5 years ago

**Status**changed from*New*to*In Progress***% Done**changed from*0*to*30*

The whole correction for disk extinction computation definitely needs to be revisited, based on new studies (e.g., SDSS). However for now, let's just keep the De Vaucouleurs et al. (1991) relation between the extinction in the blue $A_B$ and the disk inclination with respect to the line-of-sight $i$:

$$A_B = - \alpha(T)\log \cos i,$$

where $\alpha(T)$ is a type-dependent factor which corresponds to the `DISK_EXTINCT`

configuration parameter (see e.g., the last page of Bertin&Arnouts 1996). Averaging over all inclinations we get (see this link)

$$\langle A_B\rangle_i = +\alpha(T) \log 2$$.

Stuff must generate galaxies using a "face-on" luminosity function with a Schechter $M_0^*$, which should be related to the "observed" luminosity function $M^*$ using

$$M_0^* = M^* -2.5\log_{10} \left(B/T + \frac{(1-B/T)}{\rho}\right),$$

where $B/T$ is the average observed bulge-to-total flux ratio for this category of galaxies, and $\rho$ the average ratio between the observed disk flux and the face-on disk flux in the reference passband (which is assumed to be blue):

$$\rho= 10^{-0.4\,A_B} = 10^{-0.12041\,\alpha(T)}.$$

Hence, one should have, I think:

$$M_0^* = M^* -2.5\log_{10} \left(B/T + (1-B/T)10^{0.12041\,\alpha(T)}\right),$$

which does not seem to match the actual formula in the code. Note also that formally what should be averaged over all angles is $M_i^*$, not $A_B(i)$.

#### Updated by Emmanuel Bertin about 5 years ago

**% Done**changed from*30*to*40*

(!) Edited the formulas above to take into account averaging over all inclination angles.

#### Updated by Damien Le Borgne about 5 years ago

1) in the default stuff config, DISK_EXTINCT 0.0,0.75,1.23,1.47,1.47,1.23

so alpha is positive. Does it mean $A_B$ is negative ????

2) I didn't understand the code this way :

This idea, I think, as you explained to me, is that the mstar should first be shifted to brighter values in order to account for the fact that disks will be obscured later in the code.

Indeed,the bulge flux is unobscured and the disk obscured, so that the total flux variation that account for the extinction is

$\Delta F_{tot}$ = $F_{disk}^{obscured} - F_{disk}^{unobscured} = F_{disk}^{unobscured} ( \exp^{-\tau}-1)$

so that $\Delta F_{tot}/F_{tot} = (D/T)\, ( \exp^{-\tau}-1) = (1-B/T) ( \exp^{-\tau}-1)$, which is $\le 0$.

and $\Delta m = - 2.5 \log_{10} (1 + \Delta F_{tot}/F_{tot}) = -2.5 \log_{10} (1 + (1-B/T) ( \exp^{ - \tau}-1))$, which is $\ge 0$.

So the formula should be

```
/* Correct M* for disk internal extinction */
galtype->lf->mstar += 2.5*log10(1.0+(1.0-bt)*(DEXP(-0.12041*extinct)-1));
```

?

#### Updated by Damien Le Borgne about 5 years ago

From you paper, $A_B = - \alpha(T) \log \cos i$.

You are missing a $\log$ in your formula above.

#### Updated by Emmanuel Bertin about 5 years ago

**% Done**changed from*40*to*50*

Ah yes, thanks! Now at least I get the right 0.12041 factor. In my case I would only need to replace 1.0 with $B/T$.

#### Updated by Damien Le Borgne about 5 years ago

**% Done**changed from*50*to*60*

... and add the "minus" sign in the exponential (your last formula misses it).

Still, I'm not sure to agree with

$$M_0^* = M^* -2.5\log_{10} \left(B/T + (1-B/T)/\rho\right),$$

Why does the bulge ($B/T$, first term) affect M*, since it is not affected by dust?

Also, to be sure, if I measure blindly a LF from the stuff catalog, am I supposed to recover the LF defined in the stuff parameters by the Schechter function?

#### Updated by Damien Le Borgne about 5 years ago

Sorry, I guess you are right with your formula (including the minus sign of course which I overlooked).

#### Updated by Damien Le Borgne about 5 years ago

Still, I have a weird case: if $A_B \rightarrow +\infty$ (i.e. complete extinction of the disk), then $\rho \rightarrow 0$ and $M_0^* - M_0 \rightarrow -\infty$.

Is that expected?...

#### Updated by Emmanuel Bertin about 5 years ago

Still, I have a weird case: if $A_B \rightarrow +\infty$ (i.e. complete extinction of the disk), then $\rho \rightarrow 0$ and $M_0^* - M_0 \rightarrow -\infty$. Is that expected?...

Yes this is because in that case the formula with $\rho$ is no longer applicable, and you get instead

$$M_0^* = M^* -2.5\log_{10} B/T$$

#### Updated by Damien Le Borgne about 5 years ago

**File**counts.png counts.png added

With your formula above patched to stuff, I measure this in 3 stuff catalogs (details below):

with the following common parameters:

```
IMAGE_WIDTH 16384
IMAGE_HEIGHT 16384
PIXEL_SIZE 0.2
MAG_LIMITS 16.0,28.0
PASSBAND_REF couch/Bj
CALIBSED_REF AB
REFDETECT_TYPE PHOTONS
PASSBAND_OBS couch/Bj
CALIBSED_OBS AB
OBSDETECT_TYPE PHOTONS
HUBBLE_TYPE -6.0
SED_GALAXIES E
SEDINDEX_BULGE 1
SEDINDEX_DISK 1
LF_PHISTAR 4.95e-2
LF_MSTAR -19.58
LF_ALPHA -0.54
LF_MAGLIMITS -27.0,-13.0
LF_PHISTAREVOL 0.
LF_MSTAREVOL 0.
```

and the 3 sets of particular parameters

```
CATALOG_NAME bj.list
BULGE_FRACTION 1.
DISK_EXTINCT 0.0
```

```
CATALOG_NAME bj2.list
BULGE_FRACTION 0.2
DISK_EXTINCT 0.0
```

```
CATALOG_NAME bj3.list
BULGE_FRACTION 0.2
DISK_EXTINCT 10.0
```

Is it normal that bj3 is shifted bright-wards ?

#### Updated by Emmanuel Bertin about 5 years ago

This must be a side-effect. Did you turn off galaxy luminosity evolution?

#### Updated by Damien Le Borgne about 5 years ago

Yes, I set the LF evolution parameters to 0 as you can see, to make thinks easier to understand.

Note: this plot was made using $M_0^* = M^* -2.5\log_{10} \left(B/T + (1-B/T)10^{0.12041\,\alpha(T)}\right)$ only, not your simpler formula at large $A_B$.

Also, why you the formula not be valid for large $A_B$? What is the acceptable range of validity then?

Sorry to bother you... I don't understand everything.... but you can explain me next week if you prefer.

#### Updated by Emmanuel Bertin about 5 years ago

OK I have to check where side effects kick in (most probably boundaries that depend on $M_0^*$ somewhere).

#### Updated by Emmanuel Bertin about 5 years ago

**% Done**changed from*60*to*70*

Actually as this plot shows one cannot approximate $M_0^*-M^*$ satisfactorily for $\alpha(T)>>1$ using $\langle A_B \rangle_i$, especially for high $B/T$'s. So I guess I will have to numerically average $M_0^*(i)-M^*(i)$ over all inclination angles. *btw* note that $i$ never reaches $\frac{\pi}{2}$ in `stuff`

.

#### Updated by Damien Le Borgne about 5 years ago

I went through the code and did some tests.

If we stick to reasonnable alpha (`extinct`

) parameters , i.e. `extinct`

<3, I guess the current code (i.e. with $M_0^* = M^* -2.5\log_{10} \left(B/T + (1-B/T)10^{0.12041\,\alpha(T)}\right),$) is more or less good enough.

More precisely, consider this plot below:

- upper panel = histogram of mabs from the sampling of the face-on Schechter function in Stuff.
- lower panel = histogram of mabs + an inclination-dependent k-correction. (So it's almost apparent magnitudes except the distance effect is not applied)

- bj : no extinction, bulge only
- bj2 : 20% bulge, and no exinction for the disk
- bj3 : 20% bulge, and
`extinct=3`

for the disk

The bug is: **The red line (bj3) in the upper plot should not go above the other luminosity functions** : only M* changes (I checked in the code that Phistar and alpha are the same).

When extinction is applied to bj3 (lower panel), the magnitudes more or less fall in the expected range, but there are too many galaxies! (Remember the Phi* values are the same for bj, bj2, bj3, so the LFs should only be shifted along the magnitudes axis...)

In Stuff, I guess the problem comes from `lf_intschechter`

, or from the magmin, magmax variables which are used for the normalisation (?). I can't think of anything else...

Emmanuel can you please check?