

def Broaden(self, data, model, oversampling=5, m=101, dimension=20, full_output=False):
    "\n        Fits the broadening profile using singular value decomposition. This function is\n        called by GenerateModel, and is not meant to be called directly!\n\n        -oversampling is the oversampling factor to use before doing the SVD\n        -m is the size of the broadening function, in oversampled units\n        -dimension is the number of eigenvalues to keep in the broadening function. (Keeping too many starts fitting noise)\n\n        -NOTE: This function works well when there are strong telluric lines and a flat continuum.\n               If there are weak telluric lines, it's hard to not fit noise.\n               If the continuum is not very flat (i.e. from the spectrum of the actual\n                 object you are trying to telluric correct), the broadening function\n                 can become multiply-peaked and oscillatory. Use with care!\n        "
    n = (data.x.size * oversampling)
    if ((n % 2) != 0):
        n += 1
    if ((m % 2) == 0):
        m += 1
    Spectrum = UnivariateSpline(data.x, (data.y / data.cont), s=0)
    Model = UnivariateSpline(model.x, model.y, s=0)
    xnew = np.linspace(data.x[0], data.x[(- 1)], n)
    ynew = Spectrum(xnew)
    model_new = FittingUtilities.RebinData(model, xnew).y
    design = np.zeros(((n - m), m))
    for j in range(m):
        for i in range((m // 2), ((n - (m // 2)) - 1)):
            design[((i - (m // 2)), j)] = model_new[((i - j) + (m // 2))]
    design = mat(design)
    try:
        (U, W, V_t) = svd(design, full_matrices=False)
    except np.linalg.linalg.LinAlgError:
        outfilename = 'SVD_Error.log'
        outfile = open(outfilename, 'a')
        np.savetxt(outfile, np.transpose((data.x, data.y, data.cont)))
        outfile.write('\n\n\n\n\n')
        np.savetxt(outfile, np.transpose((model.x, model.y, model.cont)))
        outfile.write('\n\n\n\n\n')
        outfile.close()
        sys.exit(('SVD did not converge! Outputting data to %s' % outfilename))
    W1 = (1.0 / W)
    U_t = np.transpose(U)
    V = np.transpose(V_t)
    W1[dimension:] = 0
    W2 = diagsvd(W1, m, m)
    spec = np.transpose(mat(ynew[(m / 2):((n - (m / 2)) - 1)]))
    temp = np.dot(U_t, spec)
    temp = np.dot(W2, temp)
    Broadening = np.dot(V, temp)
    spacing = (xnew[2] - xnew[1])
    xnew = np.arange(model.x[0], model.x[(- 1)], spacing)
    model_new = Model(xnew)
    Broadening = np.array(Broadening)[(..., 0)]
    maxindex = Broadening.argmax()
    badfit = False
    if ((maxindex > ((m / 2.0) + (m / 10.0))) or (maxindex < ((m / 2.0) - (m / 10.0)))):
        outfilename = 'SVD_Error2.log'
        np.savetxt(outfilename, np.transpose((Broadening,)))
        warnings.warn('Warning! SVD Broadening function peaked at the wrong location! See SVD_Error2.log for the broadening function')
        badfit = True
    elif (np.mean(Broadening[(maxindex - int((m / 10.0))):(maxindex + int((m / 10.0)))]) < (3 * np.mean(Broadening[int((m / 5.0)):]))):
        outfilename = 'SVD_Error2.log'
        np.savetxt(outfilename, np.transpose((Broadening,)))
        warnings.warn('Warning! SVD Broadening function is not strongly peaked! See SVD_Error2.log for the broadening function')
        badfit = True
    if badfit:
        resolution = self.GetValue('resolution')
        (model, resolution) = self.FitResolution(data, model, resolution)
        if full_output:
            centralwavelength = ((data.x[0] + data.x[(- 1)]) / 2.0)
            FWHM = (centralwavelength / resolution)
            sigma = (FWHM / (2.0 * np.sqrt((2.0 * np.log(2.0)))))
            left = 0
            right = np.searchsorted(xnew, (10 * sigma))
            x = np.arange(0, (10 * sigma), (xnew[1] - xnew[0]))
            gaussian = np.exp(((- ((x - (5 * sigma)) ** 2)) / (2 * (sigma ** 2))))
            return (model, [(gaussian / gaussian.sum()), xnew])
        else:
            return model
    model = DataStructures.xypoint(x=xnew)
    Broadened = UnivariateSpline(xnew, np.convolve(model_new, Broadening, mode='same'), s=0)
    model.y = Broadened(model.x)
    params = [0.0, (- Broadening[maxindex]), maxindex, 10.0]
    (params, success) = leastsq(self.GaussianErrorFunction, params, args=(np.arange(Broadening.size), Broadening))
    sigma = (params[3] * (xnew[1] - xnew[0]))
    FWHM = ((sigma * 2.0) * np.sqrt((2.0 * np.log(2.0))))
    resolution = (np.median(data.x) / FWHM)
    logging.info(('Approximate resolution = %g' % resolution))
    if full_output:
        return (FittingUtilities.RebinData(model, data.x), [Broadening, xnew])
    else:
        return FittingUtilities.RebinData(model, data.x)

dataset/ETHPy150Open kgullikson88/Telluric-Fitter/src/TelluricFitter.py TelluricFitter.Broaden/original