The optional globular cluster is generated using a King (1962)
  star density law
  \begin{equation}
    D(r)=k((1+(r/r_{c})^{2})^{-1/2}-(1+(r_{t}/r_{c})^{2})^{-1/2})^{2}
  \end{equation}
  where $D(r)$ is the star surface density a projected distance $r$ from the
  centre of the cluster.
  $k$, $r_{c}$ and $r_{t}$ are constants, $k$ being a scale factor and $r_{c}$
  and $r_{t}$ the core and tidal radii respectively.
  
  Moffat’s formula gives the intensity I(r) at a radial distance r from the
  centre of the star image as
  \begin{equation}
    I(r)=I_{0}/(1+(\frac{r}{R})^{2})^{\beta}
  \end{equation}
  where $I_{0}$, R and $\beta$ are all constants.
  The total luminosity, Lt, of such a profile is
  \begin{equation}
    Lt=\frac{\pi R^{2}I_{0}}{\beta -1}
  \end{equation}
  so
  \begin{equation}
    I(r)=\frac{Lt(\beta -1)}{\pi^{2}R^{2}}/(1+(\frac{r}{R})^{2})^{\beta}
  \end{equation}
  
  If the intensity threshold beyond which the profile is truncated is $I_{th}$,
  then the corresponding radius $r_{th}$ is
  \begin{equation}
    r_{th}=((\frac{Lt(\beta -1)}{I_{th}\pi R^{2}})^{1/\beta}-1)^{1/2}R
  \end{equation}
  from this, the fraction of the total light emitted beyond this boundary, f, may
  be calculated
  \begin{equation}
    f=(1+(\frac{r_{th}}{R})^{2})^{1-\beta}
  \end{equation}
  \begin{displaymath}
    B=\frac {\pi( (X-(\frac{IXEXT}{2}-1))^{2} +
    (Y-\frac{IYEXT}{2})^{2} )^{\frac{1}{2}}} {A}
  \end{displaymath}

  An AR model is of the form
  \begin{equation}
    X_{i}=\sum_{j=1}^{M} A_{j}X_{i-j}+E_{i}
  \end{equation}
  for equally spaced observations $X_{i}$ and for a set of constants $A_{j}$.
  $E_{i}$ is the error in using this model.
  The method involves choosing the $A_{j}$ to minimize the $E_{i}$.
  
  The Q, U and E frames can now be calculated as follows:
  \[Q_{ij}=\frac{A_{ij}-B_{ij}}{A_{ij}+B_{ij}}\]
  \[U_{ij}=\frac{C_{ij}-D_{ij}}{C_{ij}+D_{ij}}\]
  \[E_{ij}=\frac{2}{A_{ij}+B_{ij}+C_{ij}+D_{ij}}\]
  The total polarization frame, P, and polarization angle frame, T, are generated
  from Q, U and (if required) E as follows:
  \[P_{ij}=\sqrt{{Q_{ij}}^{2}+{U_{ij}}^{2}}-E_{ij}\]
  \[T_{ij}=0.5\arctan \frac{U_{ij}}{Q_{ij}}\]
  The values of Tij are restricted to the range 0 to +$\pi$ radians.
  
  Every pixel (X,Y) in the image frame is multiplied by SIN(B)/B where
  \begin{displaymath}
    B=\frac {\pi( (X-(\frac{IXEXT}{2}-1))^{2} +
    (Y-\frac{IYEXT}{2})^{2} )^{\frac{1}{2}}} {A}
  \end{displaymath}