Create a function that maps (sub)stellar mass to radius, based on the
    famous models of Baraffe+ (1998A&A...337..403B).

    tablelines
      An iterable yielding lines from the table data file.
      I've named the file '1998A&A...337..403B_tbl1-3.dat'
      in some repositories (it's about 150K, not too bad).
    metallicity
      The metallicity of the model to select.
    heliumfrac
      The helium fraction of the model to select.
    age_gyr
      The age of the model to select, in Gyr.
    age_tol
      The tolerance on the matched age, in Gyr.

    Returns: a function mtor(mass_g), return a radius in cm as a function of a
    mass in grams. The mass must be between 0.05 and 0.7 Msun.

    The ages in the data table vary slightly at fixed metallicity and helium
    fraction. Therefore, there needs to be a tolerance parameter for matching
    the age.

    This function requires Scipy.

def mk_radius_from_mass_bcah98 (*args, **kwargs):
    """Create a function that maps (sub)stellar mass to radius, based on the
    famous models of Baraffe+ (1998A&A...337..403B).

    tablelines
      An iterable yielding lines from the table data file.
      I've named the file '1998A&A...337..403B_tbl1-3.dat'
      in some repositories (it's about 150K, not too bad).
    metallicity
      The metallicity of the model to select.
    heliumfrac
      The helium fraction of the model to select.
    age_gyr
      The age of the model to select, in Gyr.
    age_tol
      The tolerance on the matched age, in Gyr.

    Returns: a function mtor(mass_g), return a radius in cm as a function of a
    mass in grams. The mass must be between 0.05 and 0.7 Msun.

    The ages in the data table vary slightly at fixed metallicity and helium
    fraction. Therefore, there needs to be a tolerance parameter for matching
    the age.

    This function requires Scipy.

    """
    from scipy.interpolate import UnivariateSpline
    m, r = load_bcah98_mass_radius (*args, **kwargs)
    spl = UnivariateSpline (m, r, s=1)

    # This allows us to do range-checking with either scalars or vectors with
    # minimal gymnastics.
    @numutil.broadcastize (1)
    def interp (mass_g):
        if np.any (mass_g < 0.05 * cgs.msun) or np.any (mass_g > 0.7 * cgs.msun):
            raise ValueError ('mass_g must must be between 0.05 and 0.7 Msun')
        return spl (mass_g)

    return interp