What is this line doing?

[sturgman]

I have a question regarding the code in PyAstronomy/src/pyasl/asl/broad.py

Function broadGaussFast(x, y, sigma, edgeHandling=None, maxsig=None) has the following:

e = gf.evaluate(nx)
# This step ensured that the 
e /= np.sum(e)

My question is, why are we dividing by np.sum(e)?

[sczesla]

Hi,

Very good question. The comment has obviously been mutilated.

As far as I remember, this "normalization" is to ensure that, in the final result, the EWs of the original and broadened lines is the same, which appears to be the case; just re-checked that.

Stefan

[sturgman]

Sorry. I am not in this field and don't know what EW is. I was using the code as an example for my own routine that broadens using a Gaussian convolution. I am using a normalized Gaussian (area under curve equal 1) and then do the convolution like you do.

However, you set your Gaussian to have an area of 1 with the line: gf['A']=1 Correct? So the second line will divide by 1 as long as the x-axis interval is 1. If the x-axis interval is anything else, you will get poor normalization, no?

I realize this is an awkward conversation to have via this channel so just close the issue if it is not of interest.

[sczesla]

Sorry, the EW is the "equivalent width", which is the area under the curve (e.g., absorption line) normalized by the surrounding continuum (often used in astronomy). At any rate, the idea is exactly what you suggest: The area under the curve should not change as a result of the convolution/broadening, and I think the algorithm has this property. Below is a small "verification" of this:

from __future__ import print_function, division
from PyAstronomy import pyasl
import matplotlib.pylab as plt
import numpy as np
import scipy.integrate as scinteg
from PyAstronomy import funcFit as fuf

# Set up an input spectrum
x = np.linspace(5000.0,5100.0,2517)
y = np.zeros(x.size)

# Use a Gaussian as "signal"
gf = fuf.GaussFit1d()
gf["mu"] = 5050.0
gf["A"] = 1.0
gf["sig"] = 0.74

y += gf.evaluate(x)

# Apply Gaussian instrumental broadening, setting the resolution to 10000.
r, fwhm = pyasl.instrBroadGaussFast(x, y, 1000,
          edgeHandling="firstlast", fullout=True)

print("Area of input Gaussian: ", scinteg.trapz(y, x))
print("Area of Broadened Gaussian: ", scinteg.trapz(r, x))

# Apply Gaussian instrumental broadening, setting the resolution to 10000.
# Limit the extent of the Gaussian broadening kernel to five standard
# deviations.
r2, fwhm = pyasl.instrBroadGaussFast(x, y, 1000,
          edgeHandling="firstlast", fullout=True, maxsig=5.0)

print("Area of Gaussian (5sig kernel): ", scinteg.trapz(r2, x))

print("FWHM used for the Gaussian kernel: ", fwhm, " A")

# Plot the output
plt.plot(x,r, 'r--p', label="Broadened curve (full)")
plt.plot(x, r2, 'k:', label="Broadened curve (5 stds)")
plt.plot(x,y, 'b-', label="Input")
plt.legend(loc=4)
plt.show()

As far as I can see, changing the sampling (by adapting the linspace command) does not produce a change in the area prior or after the broadening. I think the reason for requiring the normalization by sum (and not area, which is, indeed, achieved by gf["A"]=1), is that the actual convolution a few lines below is carried out in bin/index space; it is not a proper integration. The algorithm is essentially ignorant of the spacing.

I think this is a perfectly suitable channel for having such discussions, and I would be most interested to correct any mistakes or logical errors in the algorithms. This is a non-trivial issue.

Stefan

[sturgman]

Thanks for the clarification. It makes perfect sense and I found no easy way of breaking your code to yield a different area. It is robust as far as I can tell as well.

[sczesla]

Great to hear!