D-term transfer

If an interferometric array contains antennas with known D-term (D-ref antennas), it is possible to determine D-term of other antennas through observations toward a polarization calibrator. The Stokes parameters of the calibrator can be assumed to be known because they can be measured using a sub array that consists of D-ref antennas.

# Formalism
Let the Stokes parameters and the correlator outputs (with Gain correction) as :

```math
\boldsymbol{S} = \left(
\begin{array}{c}
I \\
Q \\
U \\
V \\
\end{array}
\right), \ 
\boldsymbol{X} = \left(
\begin{array}{c}
\left< XX^* \right> / ( G^j_X G^{i*}_X ) \\
\left< XY^* \right> / ( G^j_X G^{i*}_Y ) \\
\left< YX^* \right> / ( G^j_Y G^{i*}_X ) \\
\left< YY^* \right> / ( G^j_Y G^{i*}_Y ) \\
\end{array}
\right).
```

and let the D-term matrix and parallactic-angle matrix as :

```math
D = 
\left(
\begin{array}{cccc}
1              & D^{i*}_X        & D^j_X           & D^j_X D^{i*}_X \\
D^{i*}_Y       & 1               & D^j_X D^{i*}_Y  & D^j_X          \\
D^j_Y          & D^j_Y D^{i*}_X  & 1               & D^{i*}_X       \\
D^j_Y D^{i*}_Y & D^j_Y           & D^{i*}_Y        & 1              \\
\end{array}
\right) \tag{1.1} \\
P = 
\left(
\begin{array}{cccc}
1 &  \cos 2\psi  &  \sin 2\psi & 0 \\
0 & -\sin 2\psi  &  \cos 2\psi & i \\
0 & -\sin 2\psi  &  \cos 2\psi & -i \\
1 & -\cos 2\psi  & -\sin 2\psi & 0 \\
\end{array}
\right)
```

then we have $\boldsymbol{X} = D P \boldsymbol{S}$. In (1.1), let $D^j$ known. The least-squared fit for unknown $D_i$ will be obtained to minimize the norm of the residual, $\boldsymbol{X} - D_{\rm trial}P\boldsymbol{S}$. The correction vector is given by $P \cdot$ residual where $P$ is the partial derivative matrix given as:

```math
P = \left(
\begin{array}{cccc}
\frac{\partial \left< XX^* \right>_{\rm real}}{\partial D_{\rm x.real}} &
\frac{\partial \left< XX^* \right>_{\rm real}}{\partial D_{\rm x.imag}} &
\frac{\partial \left< XX^* \right>_{\rm real}}{\partial D_{\rm y.real}} &
\frac{\partial \left< XX^* \right>_{\rm real}}{\partial D_{\rm y.imag}} \\
\frac{\partial \left< XX^* \right>_{\rm imag}}{\partial D_{\rm x.real}} &
\frac{\partial \left< XX^* \right>_{\rm imag}}{\partial D_{\rm x.imag}} & \cdots  \\
\frac{\partial \left< XY^* \right>_{\rm real}}{\partial D_{\rm x.real}} & \cdots  \\
\frac{\partial \left< XY^* \right>_{\rm imag}}{\partial D_{\rm x.real}} & \cdots  \\
\frac{\partial \left< YX^* \right>_{\rm real}}{\partial D_{\rm x.real}} & \cdots  \\
\frac{\partial \left< YX^* \right>_{\rm imag}}{\partial D_{\rm x.real}} & \cdots  \\
\frac{\partial \left< YY^* \right>_{\rm real}}{\partial D_{\rm x.real}} & \cdots  \\
\frac{\partial \left< YY^* \right>_{\rm imag}}{\partial D_{\rm x.real}} & \cdots  \\
\end{array} \right),
``` 

The matrix $P$ is divided into sub-matrices as 

```math
P = \left(
\begin{array}{cc}
A & O \\
O & B \\
\end{array} \right),
```

where 

```math
A = \left(
\begin{array}{cc}
UCmQS + (I - QCpUS) D_{\rm x, real} & (I - QCpUS) D_{\rm x, imag} \\
(QCpUS - I) D_{\rm x, imag} & UCmQS + (I-QCpUS) D_{\rm x, imag} \\
I - QCpUS + UCmQS \ D_{\rm y, real} & UCmQS \ D_{\rm y, imag} \\
-UCmQS \ D_{\rm y, imag} & I - QCpUS + UCmQS \ D_{\rm x, real} \\
\end{array} \right), \ B = \left(
\begin{array}{cc}
I + QCpUS + UCmQS \ D_{\rm x, real} & UCmQS \ D_{\rm x, imag} \\
-UCmQS \ D_{\rm x, imag} & I + QCpUS + UCmQS \ D_{\rm x, real} \\
UCmQS + (I + QCpUS) D_{\rm x, real} & (I + QCpUS) D_{\rm y, imag} \\
-(I + QCpUS) D_{\rm y, imag} & UCmQS + (I + QCpUS) D_{\rm y, real} \\
\end{array} \right).
```

Then, we have

```math
P^TP = \left(
\begin{array}{cc}
A^TA & O \\
O & B^TB \\
\end{array} \right)
```

Here is a Python code to calculate $P^TP$ and solution for the correction vector.

```python
def TransferD(Vis, DtX, DtY, PS):
    refAntNum, PAnum = len(DtX), PS.shape[1]
    #
    A0 = np.repeat(PS[2], refAntNum)  +  np.outer(PS[3], DtX.real).reshape(PAnum* refAntNum)
    A1 = np.repeat(PS[3], refAntNum)  +  np.outer(PS[2], DtY.real).reshape(PAnum* refAntNum)
    A2 = -np.outer(PS[3], DtX.imag).reshape(PAnum* refAntNum)
    A3 = -np.outer(PS[2], DtY.imag).reshape(PAnum* refAntNum)
    #
    B0 = np.repeat(PS[0], refAntNum)  +  np.outer(PS[1], DtX.real).reshape(PAnum* refAntNum)
    B1 = np.repeat(PS[1], refAntNum)  +  np.outer(PS[0], DtY.real).reshape(PAnum* refAntNum)
    B2 = -np.outer(PS[1], DtX.imag).reshape(PAnum* refAntNum)
    B3 = -np.outer(PS[0], DtY.imag).reshape(PAnum* refAntNum)
    #
    resid = Vis.transpose(0,2,1).reshape(4, PAnum* refAntNum)
    resid[0] -= (np.repeat(PS[0], refAntNum) + np.outer(PS[1], DtX.conjugate()).reshape(PAnum* refAntNum))
    resid[1] -= (np.repeat(PS[1], refAntNum) + np.outer(PS[0], DtY.conjugate()).reshape(PAnum* refAntNum))
    resid[2] -= (np.repeat(PS[2], refAntNum) + np.outer(PS[3], DtX.conjugate()).reshape(PAnum* refAntNum))
    resid[3] -= (np.repeat(PS[3], refAntNum) + np.outer(PS[2], DtY.conjugate()).reshape(PAnum* refAntNum))
    #
    PTP_diag = np.array([
        A0.dot(A0) + A1.dot(A1) + A2.dot(A2) + A3.dot(A3), 
        A0.dot(A0) + A1.dot(A1) + A2.dot(A2) + A3.dot(A3), 
        B0.dot(B0) + B1.dot(B1) + B2.dot(B2) + B3.dot(B3),
        B0.dot(B0) + B1.dot(B1) + B2.dot(B2) + B3.dot(B3)])
    #
    PTdotR = np.array([
        A0.dot(resid[0].real) + A1.dot(resid[1].real) + A2.dot(resid[0].imag) + A3.dot(resid[1].imag),
       -A2.dot(resid[0].real) - A3.dot(resid[1].real) + A0.dot(resid[0].imag) + A1.dot(resid[1].imag),
        B0.dot(resid[2].real) + B1.dot(resid[3].real) + B2.dot(resid[2].imag) + B3.dot(resid[3].imag),
       -B2.dot(resid[2].real) - B3.dot(resid[3].real) + B0.dot(resid[2].imag) + B1.dot(resid[3].imag)])
    #
    Solution = PTdotR / PTP_diag
    return Solution[0] + 1.0j* Solution[1], Solution[2] + 1.0j* Solution[3]
#
```