Avoiding convolution

[JohannesBuchner]

Hi Matthew,

This part of the code

        # now broaden bc to simulate effect of bulk-velocity shifts
        vsigma = vfwhm / 2.35
        wsigma = wavbe * vsigma*1e3 / _c_  # change vsigma from km/s to m/s
        winc = (self.wavlen[self.wav2num(wnorm)]
                - self.wavlen[self.wav2num(wnorm) - 1])
        psigma = wsigma / winc     # winc is wavelength increment at wnorm
        gauss = Gaussian1DKernel(stddev=psigma)
        flux_bc = convolve(flux_bc, gauss)
        # Performs a Gaussian smooth with dispersion psigma pixels

is probably not too fast because of the convolution.

On another project I spent some time looking into the BC. It turns out one can find an analytic approximating form, by considering that the BC can be approximated by a triangle (linear function with cut-off), which is convolved with a Gaussian.

            black_body_BC = (self.BC_wave**-5) / np.expm1(h_c_per_k_B / (self.BC_T * self.BC_wave))
            truncation = -np.expm1(-self.BC_tau * self.BC_wave_ratio**3)
            black_body_BC0 = (self.BE_wave**-5) / np.expm1(h_c_per_k_B / (self.BC_T * self.BE_wave))
            truncation0 = -np.expm1(-self.BC_tau)

            # the following is based on a derivation of 
            # convolving a gaussian with a linear approximation
            #    truncation ~= alpha * x + beta
            # with the values:
            alpha = 1.8
            beta = -0.8
            # this leads to a Gaussian CDF term (termB) and a
            #   more complicated result from the x * Gaussian(x) integration (termA1-3)
            sigma = (self.lines_width * 1000) / cst.c 
            x = self.BC_wave_ratio
            z = (x - 1) * 2**-0.5 / sigma
            termB = 0.5 * (1 - erf(z))
            termA1 = 0.5 * x
            termA2 = -0.5 * x * erf(z)
            termA3 = -sigma * (2 * np.pi)**-0.5 * np.exp(-z**2)
            convolved = (beta * termB + alpha * (termA1 + termA2 + termA3)) * (1 - np.exp(-1))
            # 250 nm is >3 sigma away from the balmer edge up to 30000km/s
            #   below we can use the original truncation formula
            truncation_convolved = np.where(self.BC_wave > 250, convolved, truncation)

            self.BC = black_body_BC / black_body_BC0 * truncation_convolved / truncation0

[JohannesBuchner]

You may also want to have a look at all calls to np.exp, and see if np.expm1 may not be preferrable for numerical stability.

For example np.exp(1.43877735e8 / (tbb*wav)) - 1.0 can be np.expm1(1.43877735e8 / (tbb*wav))