def dpdx(a, x, y, order=4):
    """Differential with respect to x
    The polynomial is defined as p(x,y) = a[i,j] * x**(i-j) * y**j summed over i and j
    Then dp/dx = (i-j) * a[i,j] * x**(i-j-1) * y**j


    Parameters
    ----------
    a      a linear array of polynomial coefficients in JWST order.
    x      an integer or float variable(or an array of same) representing pixel x positions
    y      a variable (or an array) representing  pixel y positions. x and y must be of the same shape.
    order  an integer, the polynomal order

    Returns
    -------
    dpdx    float value of dp/dx for the given (x,y) point(s)
    """

    dpdx = 0.0
    k = 1  # index for coefficients
    for i in range(1, order + 1):
        for j in range(i + 1):
            if i - j > 0:
                dpdx = dpdx + (i - j) * a[k] * x**(i - j - 1) * y**j
            k += 1
    return dpdx