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

    Parameters
    ----------
    a      an array of polynomial coefficients in JWST arrangement.
            The number of coefficients must be (order+1)(order+2)/2
    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
    order  an integer, the polynomal order

    Returns
    -------
    dpdy    float value of dp/dy for the given (x,y) point(s)
            where p(x,y) is the value of the polynomial
    """

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