How to differentiate a function defined in polar coordinates over discrete cartesian coordinates? EDITED :
Let's say I have a grid of points in a polar coordinate system (i.e. a polar grid, evenly distributed). This is represented on the figure : polar coordinate grid.
I do want to differentiate a function $f(r,\theta)$ with respect to $\theta$. For this I simply have to calculate :
$\frac{\partial{f(r,\theta)}}{\partial \theta}=
\frac{f(r, \,\theta_{i} + \delta_{\theta})-f(r,\,\theta_i)}{\delta_{\theta}} \quad \text{Where } i$ is the index for $\theta$ of the points of the grid.
The solution is straightforward to compute (for instance with a python algorithm).
The Problem :
I am struggling on the following problem :
In a cartesian coordinate system $(\vec{x},\vec{y})$ I have a grid of $N* N$ cartesian coordinate points $a_{ij}$
such that :
$|\vec{a_{i,j}}\cdot \vec{x}-\vec{a_{i-1,j}} \cdot \vec{x}| = \delta \quad \forall\, i \in [[-N+1, N]] \quad \text{and} \; j \in [[-N,N]]$
$|\vec{a_{i,j}}\cdot \vec{y}-\vec{a_{i,j-1}} \cdot \vec{y}| = \delta \quad \forall\, i \in [[-N, N]] \quad \text{and} \; j \in [[ -N+1,N]]$
Where $\delta$ is a constant. See cartesian coordinate grid.
Simply said,
I have a cartesian grid, evenly distributed (contiguous points are equidistant with respect to the $x$ and $y$ axis).
I would like to differentiate (and enventually integrate) the function $f(r,\theta)$ over the set of discrete cartesian points with respect to $\theta$. However, these points are no more evenly distributed for a polar differentiation.
To differentiate $f$ with respect to $\theta$ other the imposed cartesian grid ($a_{i,j}$) I get lost.
$\frac{d f(r,\theta)}{d \theta}= \cdot\cdot\cdot  = \;?$
I have been searching a lot, but have not recognised the same problem anywhere. However I feel like the solution lies within the use of total derivative,
and/or the use of Jacobian.
Important note :
$f(r,\theta)$ is considered a black box (not known).
I thank you kindly for your answer,
 A: As I understand it, you have a grid of points $(\boldsymbol x_{i,j})$ and at each point you have a value for your function $f$, $f_{i,j}$. You want to approximate $\frac{\partial f}{\partial\theta}$ at the grid points. Here's how I would do it.
Recall that the position vector in polar coordinates is
$$\boldsymbol x=\begin{pmatrix}x^1 \\ x^2\end{pmatrix}=\begin{pmatrix}r\cos\theta \\ r\sin\theta\end{pmatrix}$$
The coordinates are related by
$$\begin{matrix}x^1 = r\cos\theta \\ x^2=r\sin\theta\end{matrix}~~~~,~~~~\begin{matrix}r = \sqrt{x^2+y^2} \\ \theta=\operatorname{atan2}(x,y)\end{matrix}$$
And the polar unit basis vectors are
$$\boldsymbol e_r=\left\Vert\frac{\partial\boldsymbol x}{\partial r}\right\Vert^{-1}\frac{\partial\boldsymbol x}{\partial r}=\begin{pmatrix}\cos\theta \\ \sin\theta\end{pmatrix} 
\\ \boldsymbol e_\theta=\left\Vert\frac{\partial\boldsymbol x}{\partial \theta}\right\Vert^{-1}\frac{\partial\boldsymbol x}{\partial \theta}=\begin{pmatrix}-r\sin\theta \\ r\cos\theta\end{pmatrix}$$
They are related to the Cartesian unit basis vectors by the expressions
$$\begin{bmatrix}\boldsymbol e_r \\ \boldsymbol e_\theta\end{bmatrix}=\begin{bmatrix}\cos\theta & \sin\theta \\ -\sin\theta 
 & \cos\theta\end{bmatrix}\begin{bmatrix}\boldsymbol e_1 \\ \boldsymbol e_2\end{bmatrix}=\frac{1}{\sqrt{(x^1)^2+(x^2)^2}}\begin{bmatrix}x^1 & x^2 \\ -x^2 & x^1\end{bmatrix}\begin{bmatrix}\boldsymbol e_1\\ \boldsymbol e_2\end{bmatrix}$$
Where does this get us? Well, since we have the unit vector for $\theta$, we can express
$$\frac{\partial f}{\partial\theta}=\boldsymbol e_\theta\cdot\nabla f$$
But of course we can represent both of these using our rectangular coordinates:
$$\boldsymbol e_\theta=\frac{1}{\sqrt{(x^1)^2+(x^2)^2}} \left(-x^2\boldsymbol e_1+x^1\boldsymbol e^2\right) \\ 
\nabla f=\begin{bmatrix}\frac{\partial f}{\partial x^1} \\ \frac{\partial f}{\partial x^2}\end{bmatrix}$$
Therefore, we can write the gradient at the grid point $(i,j)$ using e.g central differences. (I'm assuming, as you've mentioned, that you have a uniform grid, but it doesn't have to be that way.)
$$(\nabla f)_{i,j}=\frac{1}{2\delta}\begin{bmatrix}f_{i+1,j}-f_{i-1,j} \\ f_{i,j+1}-f_{i,j-1}\end{bmatrix}$$
Therefore
$$\left(\frac{\partial f}{\partial\theta}\right)_{i,j}=\frac{1}{\sqrt{(x^1_{i,j})^2+(x^2_{i,j})^2}}\left(-x^2_{i,j}\begin{bmatrix}1 \\ 0\end{bmatrix}+x^1_{i,j}\begin{bmatrix}0 \\ 1\end{bmatrix}\right)\frac{1}{2\delta}\begin{bmatrix}f_{i+1,j}-f_{i-1,j} \\ f_{i,j+1}-f_{i,j-1}\end{bmatrix}$$
Simplifying,
$$\boxed{\left(\frac{\partial f}{\partial\theta}\right)_{i,j} = \frac{1}{2\delta\sqrt{(x^1_{i,j})^2+(x^2_{i,j})^2}}\bigg(x^1_{i,j}(f_{i,j+1}-f_{i,j-1})+x^2_{i,j}(f_{i-1,j}-f_{i+1,j})\bigg)}$$
Easy enough to implement I hope.
