Алгоритм

1. Как и с астигматизмом, сначала рассчитываем параметры радиально-тангенциального поля в точке (можно расширить функцию для астигаматизма
2. Для каждой точки свертки пересчитываем, куда она смещается

struct AstigCtx {
    float2 e_r;  // radial unit in uv metric
    float2 e_t;  // tangential unit (perp)
    float  s_r;  // radial scale
    float  s_t;  // tangential scale
    float  rhoN; // normalized field height (0 center .. 1 edge)
};

inline AstigCtx astigSetup(float2 uv, constant Parameters& P) {
    AstigCtx C;

    float ar = P.invSize.y / P.invSize.x; // width/height
    float2 oc = P.opticalCenter;
    float2 df = uv - oc;
    float2 dfIso = float2(df.x * ar, df.y);
    float  rhoN = clamp(length(dfIso) * 2.0f, 0.0f, 1.0f);

    float2 rIso = (length(dfIso) > 1e-6f) ? normalize(dfIso) : float2(1.0f, 0.0f);
    float2 e_r  = normalize(float2(rIso.x / ar, rIso.y));
    float2 e_t  = float2(-e_r.y, e_r.x);

    float mag   = fabs(P.astigK) * (rhoN * rhoN);
    float a     = 1.0f + mag;
    float s_r   = (P.astigK >= 0.0f) ? (1.0f / a) : a;
    float s_t   = (P.astigK >= 0.0f) ? a : (1.0f / a);

    C.e_r = e_r;
    C.e_t = e_t;
    C.s_r = s_r;
    C.s_t = s_t;
    C.rhoN = rhoN;

    return C;
}

inline float2 comaApply(float2 d, AstigCtx C, float comaK) {
    if (fabs(comaK) < 1e-6f) return d;
    // Project onto radial/tangential axes
    float dr = dot(d, C.e_r);
    float dt = dot(d, C.e_t);
    // Magnitude grows ~ rho^3 across the field (third-order behavior)
    float mag = comaK * (C.rhoN * C.rhoN * C.rhoN);
    // Quadratic skew in pupil coordinates → "comet" shape in image PSF
    float2 delta = ( (dr*dr - 0.5f*dt*dt) * C.e_r
                   + (dr*dt) * C.e_t );
    return d + mag * delta;
}