SUBROUTINE sla_FITXY (ITYPE,NP,XYE,XYM,COEFFS,J) *+ * - - - - - - * F I 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 = singular 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 sla_PXY, sla_INVF, sla_XY2XY, sla_DCMPF * * Called: sla_DMAT, sla_DMXV * * P.T.Wallace Starlink 11 February 1991 * * Copyright (C) 1995 Rutherford Appleton Laboratory *- 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 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,A,B,C,D, : SDR2,XR,YR,AOLD,BOLD,COLD,DOLD,SOLD * Preset the status J=0 * 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 sla_DMAT(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 sla_DMXV(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 sla_DMAT(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