Index: trunk/doc/pslib/psLibADD.tex
===================================================================
--- trunk/doc/pslib/psLibADD.tex	(revision 2779)
+++ trunk/doc/pslib/psLibADD.tex	(revision 3070)
@@ -1,3 +1,3 @@
-%%% $Id: psLibADD.tex,v 1.56 2004-12-21 21:37:08 price Exp $
+%%% $Id: psLibADD.tex,v 1.57 2005-01-22 01:57:42 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 
@@ -1646,5 +1646,5 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\subsubsection{Missing and Todo}
+\subsection{Missing and Todo}
 
 \tbd{define SINC, LAGRANGE interpolation}
@@ -1659,4 +1659,263 @@
 
 \tbd{define Brent's method \& minimization bracketing}
+
+\section{Pan-STARRS Modules}
+
+\subsection{Object Models}
+
+\subsubsection{Real 2D Gaussian}
+
+This function is a two-dimensional Gaussian with an elliptical
+cross-section and a constant local background:
+\[
+f(x,y) = Z_o e^{-z} + S_o
+\]
+where
+\[
+z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
+\]
+
+Below is the relationship between the \code{psModel} parameters and
+the function parameters, sample C-code implementing the function
+efficiently, and the value of the derivatives:
+
+\begin{verbatim}
+  param[0] = So;
+  param[1] = Zo;
+  param[2] = Xo;
+  param[3] = Yo;
+  param[4] = sqrt(2) / SigmaX;
+  param[5] = sqrt(2) / SigmaY;
+  param[6] = Sxy;
+
+  X = x[0] - param[2];
+  Y = x[1] - param[3];
+  
+  px = param[4]*X;
+  py = param[5]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
+  r = exp(-z);
+  f = param[1]*r + param[0];
+  /* f is the function value */
+
+  q = param[1]*r;
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q*(2*px*param[4] + param[6]*Y);
+  deriv[3] = q*(2*py*param[5] + param[6]*X);
+  deriv[4] = -2*q*px*X;
+  deriv[5] = -2*q*py*Y;
+  deriv[6] = -q*X*Y;
+\end{verbatim}
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\subsubsection{Pseudo-Gaussian}
+
+This function is a polynomial approximation of a 2D Gaussian.  The
+function is very similar to the real Gaussian:
+\[
+f(x,y) = Z_o (1 + z + z^2/2 + z^3/6)^{-1} + S_o
+\]
+where
+\[
+z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
+\]
+
+Below is the relationship between the \code{psModel} parameters and
+the function parameters, sample C-code implementing the function
+efficiently, and the value of the derivatives:
+
+\begin{verbatim}
+  param[0] = So;
+  param[1] = Zo;
+  param[2] = Xo;
+  param[3] = Yo;
+  param[4] = sqrt(2) / SigmaX;
+  param[5] = sqrt(2) / SigmaY;
+  param[6] = Sxy;
+
+  X = x[0] - param[2];
+  Y = x[1] - param[3];
+  
+  px = param[4]*X;
+  py = param[5]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
+  t = 1 + z + 0.5*z*z;
+  r = 1.0 / (t*(1 + z/3)); /* ~ exp (-Z) */
+  f = param[1]*r + param[0];
+  /* f is the function value */
+
+  /* note difference from a pure gaussian: q = param[1]*r */
+  q = param[1]*r*r*t;
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q*(2*px*param[4] + param[6]*Y);
+  deriv[3] = q*(2*py*param[5] + param[6]*X);
+  deriv[4] = -2*q*px*X;
+  deriv[5] = -2*q*py*Y;
+  deriv[6] = -q*X*Y;
+\end{verbatim}
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\subsubsection{Waussian}
+
+The Waussian is a modified polynomial approximation of a 2D Gaussian,
+with non-linear polynomial terms having variable coefficients, rather
+than the Taylor series values of 1/2 and 1/6.  The
+function is very similar to the pseudo-Gaussian:
+\[
+f(x,y) = Z_o (1 + z + B_2 (z^2/2 + B_3 z^3/6))^{-1} + S_o
+\]
+where
+\[
+z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
+\]
+
+Below is the relationship between the \code{psModel} parameters and
+the function parameters, sample C-code implementing the function
+efficiently, and the value of the derivatives.  Note the fudge factors
+of 100 in the derivatives of $B_2$ and $B_3$: these are included to
+slow the variation of these parameters, which are otherwise very
+sensitive to small errors.
+
+\begin{verbatim}
+  param[0] = So;
+  param[1] = Zo;
+  param[2] = Xo;
+  param[3] = Yo;
+  param[4] = Sx;
+  param[5] = Sy;
+  param[6] = Sxy;
+  param[7] = B2;
+  param[8] = B3;
+
+  X = x - param[2];
+  Y = y - param[2];
+  
+  px = param[4]*X;
+  py = param[5]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
+  t = 0.5*z*z*(1 + param[8]*z/3);
+  r = 1.0 / (1 + z + param[7]*t); /* ~ exp (-Z) */
+  f = param[1]*r + param[0];
+
+  /* note difference from gaussian: q = param[1]*r */
+  q = param[1]*r*r*(1 + param[7]*z*(1 + param[8]*z/2));
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q*(2*px*param[4] + param[6]*Y);
+  deriv[3] = q*(2*py*param[5] + param[6]*X);
+  deriv[4] = -2*q*px*X;
+  deriv[5] = -2*q*py*Y;
+  deriv[6] = -q*X*Y;
+  deriv[7] = -100*param[1]*r*r*t;
+  deriv[8] = -100*param[1]*r*r*param[7]*(z*z*z)/6;
+  /* the values of 100 dampen the swing of param[7,8] */
+\end{verbatim}
+
+\subsubsection{Twisted Gaussian}
+
+This function describes an object with power-law wings and a flattened
+core, where the core has a different contour from the wings.  
+
+\[
+f(x,y) = Z_{\rm pk} (1 + z_1 + z_2^M)^{-1} + Sky
+\]
+where
+\[
+z_1 = \frac{x^2}{2\sigma_{x,in}^2} + \frac{y^2}{2\sigma_{y,in}^2} + x y \sigma_{xy,in}
+z_2 = \frac{x^2}{2\sigma_{x,out}^2} + \frac{y^2}{2\sigma_{y,out}^2} + x y \sigma_{xy,out}
+\]
+
+\begin{verbatim}
+  param[0]  = So;
+  param[1]  = Zo;
+  param[2]  = Xo;
+  param[3]  = Yo;
+  param[4]  = SxInner;
+  param[5]  = SyInner;
+  param[6]  = SxyInner;
+  param[7]  = SxOuter;
+  param[8]  = SyOuter;
+  param[9]  = SxyOuter;
+  param[10] = N;
+
+  X = x - param[2];
+  Y = y - param[3];
+  
+  px1 = param[4]*X;
+  py1 = param[5]*Y;
+  px2 = param[7]*X;
+  py2 = param[8]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + param[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + param[9]*X*Y;
+
+  r  = 1.0 / (1 + z1 + pow(z2,param[10]));
+  f  = param[5]*r + param[6];
+
+  q1 = param[5]*SQ(r);
+  q2 = param[5]*SQ(r)*param[10]*pow(z2,(param[10]-1));
+
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q1*(2*px1*param[4] + param[6]*Y) + q2*(2*px2*param[7] + param[9]*Y);
+  deriv[3] = q1*(2*py1*param[5] + param[6]*X) + q2*(2*py2*param[8] + param[9]*X);
+
+  /* these fudge factors impede the growth of param[4] beyond param[7] */
+  f1 = fabs(param[7]) / fabs(param[4]);
+  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+  deriv[4] = -2*q1*px1*X*f2;
+
+  /* these fudge factors impede the growth of param[5] beyond param[8] */
+  f1 = fabs(param[8]) / fabs(param[5]);
+  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+  deriv[5] = -2*q1*py1*Y*f2;
+
+  deriv[6] = -q1*X*Y;
+
+  deriv[7] = -2*q2*px2*X;
+  deriv[8] = -2*q2*py2*Y;
+  deriv[9] = -q2*X*Y;
+  deriv[10] = -q1*ln(z2);
+\end{verbatim}
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\tbd{future galaxy models to be implemented}
+
+\begin{verbatim}
+float Sersic()
+  param[0] = So;
+  param[1] = Zo;
+  param[2] = Xo;
+  param[3] = Yo;
+  param[4] = Sx;
+  param[5] = Sy;
+  param[6] = Sxy;
+  param[7] = Nexp;
+
+float SersicBulge()
+  param[0]  So;
+  param[1]  Zo;
+  param[2]  Xo;
+  param[3]  Yo;
+  param[4]  SxInner;
+  param[5]  SyInner;
+  param[6]  SxyInner;
+  param[7]  Zd;
+  param[8]  SxOuter;
+  param[9]  SyOuter;
+  param[10] = SxyOuter;
+  param[11] = Nexp;
+\end{verbatim}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%