#include <cmath>
#include <vector>
#include <algorithm>

double jinc(double r) {
    // Returns the jinc function [J1(2*pi*r)/(2*pi*r)] evaluated at r.
    // Per this definition, the first zero of jinc function is at 0.61.
    if (r == 0.0) {
        return 0.5;
    }
    return std::real(std::cyl_bessel_j(1, r) / r);
}

std::vector<double> jinc(const std::vector<double>& r) {
    std::vector<double> y(r.size());
    std::transform(r.begin(), r.end(), y.begin(), [](double x) {
        return jinc(x);
    });
    return y;
}