Initial commit

1 files changed, 329 insertions, 0 deletions

diff --git a/math/slalib/fitxy.f b/math/slalib/fitxy.f
new file mode 100644
index 00000000..5060510e
--- /dev/null
+++ b/math/slalib/fitxy.f
@@ -0,0 +1,329 @@
+      SUBROUTINE slFTXY (ITYPE,NP,XYE,XYM,COEFFS,J)
+*+
+*     - - - - - -
+*      F T X Y
+*     - - - - - -
+*
+*  Fit a linear model to relate two sets of [X,Y] coordinates.
+*
+*  Given:
+*     ITYPE    i        type of model: 4 or 6 (note 1)
+*     NP       i        number of samples (note 2)
+*     XYE     d(2,np)   expected [X,Y] for each sample
+*     XYM     d(2,np)   measured [X,Y] for each sample
+*
+*  Returned:
+*     COEFFS  d(6)      coefficients of model (note 3)
+*     J        i        status:  0 = OK
+*                               -1 = illegal ITYPE
+*                               -2 = insufficient data
+*                               -3 = no solution
+*
+*  Notes:
+*
+*  1)  ITYPE, which must be either 4 or 6, selects the type of model
+*      fitted.  Both allowed ITYPE values produce a model COEFFS which
+*      consists of six coefficients, namely the zero points and, for
+*      each of XE and YE, the coefficient of XM and YM.  For ITYPE=6,
+*      all six coefficients are independent, modelling squash and shear
+*      as well as origin, scale, and orientation.  However, ITYPE=4
+*      selects the "solid body rotation" option;  the model COEFFS
+*      still consists of the same six coefficients, but now two of
+*      them are used twice (appropriately signed).  Origin, scale
+*      and orientation are still modelled, but not squash or shear -
+*      the units of X and Y have to be the same.
+*
+*  2)  For NC=4, NP must be at least 2.  For NC=6, NP must be at
+*      least 3.
+*
+*  3)  The model is returned in the array COEFFS.  Naming the
+*      elements of COEFFS as follows:
+*
+*                  COEFFS(1) = A
+*                  COEFFS(2) = B
+*                  COEFFS(3) = C
+*                  COEFFS(4) = D
+*                  COEFFS(5) = E
+*                  COEFFS(6) = F
+*
+*      the model is:
+*
+*            XE = A + B*XM + C*YM
+*            YE = D + E*XM + F*YM
+*
+*      For the "solid body rotation" option (ITYPE=4), the
+*      magnitudes of B and F, and of C and E, are equal.  The
+*      signs of these coefficients depend on whether there is a
+*      sign reversal between XE,YE and XM,YM;  fits are performed
+*      with and without a sign reversal and the best one chosen.
+*
+*  4)  Error status values J=-1 and -2 leave COEFFS unchanged;
+*      if J=-3 COEFFS may have been changed.
+*
+*  See also slPXY, slINVF, slXYXY, slDCMF
+*
+*  Called:  slDMAT, slDMXV
+*
+*  Last revision:   8 September 2005
+*
+*  Copyright P.T.Wallace.  All rights reserved.
+*
+*  License:
+*    This program is free software; you can redistribute it and/or modify
+*    it under the terms of the GNU General Public License as published by
+*    the Free Software Foundation; either version 2 of the License, or
+*    (at your option) any later version.
+*
+*    This program is distributed in the hope that it will be useful,
+*    but WITHOUT ANY WARRANTY; without even the implied warranty of
+*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*    GNU General Public License for more details.
+*
+*    You should have received a copy of the GNU General Public License
+*    along with this program (see SLA_CONDITIONS); if not, write to the
+*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+*    Boston, MA  02110-1301  USA
+*
+*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
+*-
+
+      IMPLICIT NONE
+
+      INTEGER ITYPE,NP
+      DOUBLE PRECISION XYE(2,NP),XYM(2,NP),COEFFS(6)
+      INTEGER J
+
+      INTEGER I,JSTAT,IW(4),NSOL
+      DOUBLE PRECISION A,B,C,D,AOLD,BOLD,COLD,DOLD,SOLD,
+     :                 P,SXE,SXEXM,SXEYM,SYE,SYEYM,SYEXM,SXM,
+     :                 SYM,SXMXM,SXMYM,SYMYM,XE,YE,
+     :                 XM,YM,V(4),DM3(3,3),DM4(4,4),DET,
+     :                 SGN,SXXYY,SXYYX,SX2Y2,SDR2,XR,YR
+
+
+
+*  Preset the status
+      J=0
+
+*  Variable initializations to avoid compiler warnings
+      A = 0D0
+      B = 0D0
+      C = 0D0
+      D = 0D0
+      AOLD = 0D0
+      BOLD = 0D0
+      COLD = 0D0
+      DOLD = 0D0
+      SOLD = 0D0
+
+*  Float the number of samples
+      P=DBLE(NP)
+
+*  Check ITYPE
+      IF (ITYPE.EQ.6) THEN
+
+*
+*     Six-coefficient linear model
+*     ----------------------------
+
+*     Check enough samples
+         IF (NP.GE.3) THEN
+
+*     Form summations
+         SXE=0D0
+         SXEXM=0D0
+         SXEYM=0D0
+         SYE=0D0
+         SYEYM=0D0
+         SYEXM=0D0
+         SXM=0D0
+         SYM=0D0
+         SXMXM=0D0
+         SXMYM=0D0
+         SYMYM=0D0
+         DO I=1,NP
+            XE=XYE(1,I)
+            YE=XYE(2,I)
+            XM=XYM(1,I)
+            YM=XYM(2,I)
+            SXE=SXE+XE
+            SXEXM=SXEXM+XE*XM
+            SXEYM=SXEYM+XE*YM
+            SYE=SYE+YE
+            SYEYM=SYEYM+YE*YM
+            SYEXM=SYEXM+YE*XM
+            SXM=SXM+XM
+            SYM=SYM+YM
+            SXMXM=SXMXM+XM*XM
+            SXMYM=SXMYM+XM*YM
+            SYMYM=SYMYM+YM*YM
+         END DO
+
+*        Solve for A,B,C in  XE = A + B*XM + C*YM
+            V(1)=SXE
+            V(2)=SXEXM
+            V(3)=SXEYM
+            DM3(1,1)=P
+            DM3(1,2)=SXM
+            DM3(1,3)=SYM
+            DM3(2,1)=SXM
+            DM3(2,2)=SXMXM
+            DM3(2,3)=SXMYM
+            DM3(3,1)=SYM
+            DM3(3,2)=SXMYM
+            DM3(3,3)=SYMYM
+            CALL slDMAT(3,DM3,V,DET,JSTAT,IW)
+            IF (JSTAT.EQ.0) THEN
+               DO I=1,3
+                  COEFFS(I)=V(I)
+               END DO
+
+*           Solve for D,E,F in  YE = D + E*XM + F*YM
+               V(1)=SYE
+               V(2)=SYEXM
+               V(3)=SYEYM
+               CALL slDMXV(DM3,V,COEFFS(4))
+
+            ELSE
+
+*           No 6-coefficient solution possible
+               J=-3
+
+            END IF
+
+         ELSE
+
+*        Insufficient data for 6-coefficient fit
+            J=-2
+
+         END IF
+
+      ELSE IF (ITYPE.EQ.4) THEN
+
+*
+*     Four-coefficient solid body rotation model
+*     ------------------------------------------
+
+*     Check enough samples
+         IF (NP.GE.2) THEN
+
+*        Try two solutions, first without then with flip in X
+            DO NSOL=1,2
+               IF (NSOL.EQ.1) THEN
+                  SGN=1D0
+               ELSE
+                  SGN=-1D0
+               END IF
+
+*           Form summations
+               SXE=0D0
+               SXXYY=0D0
+               SXYYX=0D0
+               SYE=0D0
+               SXM=0D0
+               SYM=0D0
+               SX2Y2=0D0
+               DO I=1,NP
+                  XE=XYE(1,I)*SGN
+                  YE=XYE(2,I)
+                  XM=XYM(1,I)
+                  YM=XYM(2,I)
+                  SXE=SXE+XE
+                  SXXYY=SXXYY+XE*XM+YE*YM
+                  SXYYX=SXYYX+XE*YM-YE*XM
+                  SYE=SYE+YE
+                  SXM=SXM+XM
+                  SYM=SYM+YM
+                  SX2Y2=SX2Y2+XM*XM+YM*YM
+               END DO
+
+*
+*           Solve for A,B,C,D in:  +/- XE = A + B*XM - C*YM
+*                                    + YE = D + C*XM + B*YM
+               V(1)=SXE
+               V(2)=SXXYY
+               V(3)=SXYYX
+               V(4)=SYE
+               DM4(1,1)=P
+               DM4(1,2)=SXM
+               DM4(1,3)=-SYM
+               DM4(1,4)=0D0
+               DM4(2,1)=SXM
+               DM4(2,2)=SX2Y2
+               DM4(2,3)=0D0
+               DM4(2,4)=SYM
+               DM4(3,1)=SYM
+               DM4(3,2)=0D0
+               DM4(3,3)=-SX2Y2
+               DM4(3,4)=-SXM
+               DM4(4,1)=0D0
+               DM4(4,2)=SYM
+               DM4(4,3)=SXM
+               DM4(4,4)=P
+               CALL slDMAT(4,DM4,V,DET,JSTAT,IW)
+               IF (JSTAT.EQ.0) THEN
+                  A=V(1)
+                  B=V(2)
+                  C=V(3)
+                  D=V(4)
+
+*              Determine sum of radial errors squared
+                  SDR2=0D0
+                  DO I=1,NP
+                     XM=XYM(1,I)
+                     YM=XYM(2,I)
+                     XR=A+B*XM-C*YM-XYE(1,I)*SGN
+                     YR=D+C*XM+B*YM-XYE(2,I)
+                     SDR2=SDR2+XR*XR+YR*YR
+                  END DO
+
+               ELSE
+
+*              Singular: set flag
+                  SDR2=-1D0
+
+               END IF
+
+*           If first pass and non-singular, save variables
+               IF (NSOL.EQ.1.AND.JSTAT.EQ.0) THEN
+                  AOLD=A
+                  BOLD=B
+                  COLD=C
+                  DOLD=D
+                  SOLD=SDR2
+               END IF
+
+            END DO
+
+*        Pick the best of the two solutions
+            IF (SOLD.GE.0D0.AND.(SOLD.LE.SDR2.OR.NP.EQ.2)) THEN
+               COEFFS(1)=AOLD
+               COEFFS(2)=BOLD
+               COEFFS(3)=-COLD
+               COEFFS(4)=DOLD
+               COEFFS(5)=COLD
+               COEFFS(6)=BOLD
+            ELSE IF (JSTAT.EQ.0) THEN
+               COEFFS(1)=-A
+               COEFFS(2)=-B
+               COEFFS(3)=C
+               COEFFS(4)=D
+               COEFFS(5)=C
+               COEFFS(6)=B
+            ELSE
+
+*           No 4-coefficient fit possible
+               J=-3
+            END IF
+         ELSE
+
+*        Insufficient data for 4-coefficient fit
+            J=-2
+         END IF
+      ELSE
+
+*     Illegal ITYPE - not 4 or 6
+         J=-1
+      END IF
+
+      END