How should refpts work? about at HOT 5 CLOSED

def linopt(ring,
            dp=0.0,
            refpts=None,
            get_chrom=False,
            orbit=None,
            keep_lattice=False,
            ddp=DDP,
            coupled=False):

    uintrefs = uint32_refpts([] if refpts is None else refpts, len(ring))

    if orbit is None:
        orbit = find_orbit4(ring, dp, keep_lattice=keep_lattice)
        keep_lattice = True
    orbs = numpy.squeeze(lattice_pass(ring, orbit.copy(order='K'), refpts=refpts,
                                      keep_lattice=keep_lattice))
    m44, mstack = find_m44(ring, dp, uintrefs, orbit=orbit, keep_lattice=True)
    nrefs = uintrefs.size

    M = m44[:2, :2]
    N = m44[2:, 2:]
    m = m44[:2, 2:]
    n = m44[2:, :2]

    if coupled:
        # Calculate A, B, C, gamma at the first element
        H = m + md((_s2, n.T, _s2.T))
        t = numpy.trace(M - N)
        t2 = t * t
        t2h = t2 + 4.0 * numpy.linalg.det(H)

        g = sqrt(1.0 + sqrt(t2 / t2h)) / sqrt(2.0)
        G = numpy.diag((g, g))
        C = -H * numpy.sign(t) / (g * sqrt(t2h))
        A = md((G, G, M)) - numpy.dot(G, (md((m, _s2, C.T, _s2.T)) + md((C, n)))) + md((C, N, _s2, C.T, _s2.T))
        B = md((G, G, N)) + numpy.dot(G, (md((_s2, C.T, _s2.T, m)) + md((n, C)))) + md((_s2, C.T, _s2.T, M, C))
        # Get initial twiss parameters
        a0_a, b0_a, tune_a = _closure(A)
        a0_b, b0_b, tune_b = _closure(B)
    else:
        A = M
        B = N
    # Get initial twiss parameters
    a0_a, b0_a, tune_a = _closure(M)
    a0_b, b0_b, tune_b = _closure(N)
    tune = numpy.array([tune_a, tune_b])

    # Calculate chromaticity by calling this function again at a slightly
    # different momentum.
    if get_chrom:
        d0_up, tune_up, _, l_up = linopt0(ring, dp + 0.5 * ddp, uintrefs, keep_lattice=True, coupled=coupled)
        d0_down, tune_down, _, l_down = linopt0(ring, dp - 0.5 * ddp, uintrefs, keep_lattice=True, coupled=coupled)
        chrom = (tune_up - tune_down) / ddp
        dispersion = (l_up['closed_orbit'] - l_down['closed_orbit'])[:, :4] / ddp
        disp0 = (d0_up['closed_orbit'] - d0_down['closed_orbit'])[:4] / ddp
    else:
        chrom = None
        dispersion = numpy.NaN
        disp0 = numpy.NaN

    lindata0 = numpy.zeros((), dtype=LINDATA_DTYPE)
    lindata0['closed_orbit'] = orbit
    lindata0['dispersion'] = disp0
    lindata0['alpha'] = numpy.array((a0_a, a0_b))
    lindata0['beta'] = numpy.array((b0_a, b0_b))

    # Propagate to reference points
    if nrefs > 0:
        if coupled:
            MSA, MSB, gamma, CL, AL, BL = zip(*map(analyze1, numpy.split(mstack, mstack.shape[2], axis=2), repeat(A), repeat(B), repeat(C), repeat(G)))
            msa = numpy.stack(MSA, axis=2)
            msb = numpy.stack(MSB, axis=2)
        else:
            msa = mstack[:2, :2, :]
            msb = mstack[2:, 2:, :]
            CL = numpy.zeros((2, 2, mstack.shape[2]))
            AL = msa
            BL = msa
            gamma = 1
        alpha_a, beta_a, mu_a = _twiss22(msa, a0_a, b0_a)
        alpha_b, beta_b, mu_b = _twiss22(msb, a0_b, b0_b)

        lindata = numpy.zeros(nrefs, dtype=LINDATA_DTYPE)
        lindata['idx'] = uintrefs
        # Use rollaxis to get the arrays in the correct shape for the lindata
        # structured array - that is, with nrefs as the first dimension.
        lindata['s_pos'] = get_s_pos(ring, uintrefs)
        lindata['closed_orbit'] = numpy.rollaxis(orbs, -1)
        lindata['m44'] = numpy.rollaxis(mstack, -1)
        lindata['alpha'] = numpy.stack((alpha_a, alpha_b), axis=1)
        lindata['beta'] = numpy.stack((beta_a, beta_b), axis=1)
        lindata['mu'] = numpy.stack((mu_a, mu_b), axis=1)
        lindata['A'] = numpy.rollaxis(AL, -1)
        lindata['B'] = numpy.rollaxis(BL, -1)
        lindata['C'] = numpy.rollaxis(CL, -1)
        lindata['gamma'] = gamma
        lindata['dispersion'] = dispersion
    else:
        lindata = numpy.array([], dtype=LINDATA_DTYPE)

    return lindata0, tune, chrom, lindata