diff options
Diffstat (limited to 'math/slalib/dcmpf.f')
| -rw-r--r-- | math/slalib/dcmpf.f | 160 |
1 files changed, 160 insertions, 0 deletions
diff --git a/math/slalib/dcmpf.f b/math/slalib/dcmpf.f new file mode 100644 index 00000000..93d0d0c2 --- /dev/null +++ b/math/slalib/dcmpf.f @@ -0,0 +1,160 @@ + SUBROUTINE slDCMF (COEFFS,XZ,YZ,XS,YS,PERP,ORIENT) +*+ +* - - - - - - +* D C M F +* - - - - - - +* +* Decompose an [X,Y] linear fit into its constituent parameters: +* zero points, scales, nonperpendicularity and orientation. +* +* Given: +* COEFFS d(6) transformation coefficients (see note) +* +* Returned: +* XZ d x zero point +* YZ d y zero point +* XS d x scale +* YS d y scale +* PERP d nonperpendicularity (radians) +* ORIENT d orientation (radians) +* +* Called: slDA1P +* +* The model relates two sets of [X,Y] coordinates as follows. +* Naming the elements of COEFFS: +* +* COEFFS(1) = A +* COEFFS(2) = B +* COEFFS(3) = C +* COEFFS(4) = D +* COEFFS(5) = E +* COEFFS(6) = F +* +* the model transforms coordinates [X1,Y1] into coordinates +* [X2,Y2] as follows: +* +* X2 = A + B*X1 + C*Y1 +* Y2 = D + E*X1 + F*Y1 +* +* The transformation can be decomposed into four steps: +* +* 1) Zero points: +* +* x' = XZ + X1 +* y' = YZ + Y1 +* +* 2) Scales: +* +* x'' = XS*x' +* y'' = YS*y' +* +* 3) Nonperpendicularity: +* +* x''' = cos(PERP/2)*x'' + sin(PERP/2)*y'' +* y''' = sin(PERP/2)*x'' + cos(PERP/2)*y'' +* +* 4) Orientation: +* +* X2 = cos(ORIENT)*x''' + sin(ORIENT)*y''' +* Y2 =-sin(ORIENT)*y''' + cos(ORIENT)*y''' +* +* See also slFTXY, slPXY, slINVF, slXYXY +* +* P.T.Wallace Starlink 19 December 2001 +* +* Copyright (C) 2001 Rutherford Appleton Laboratory +* +* 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 + + DOUBLE PRECISION COEFFS(6),XZ,YZ,XS,YS,PERP,ORIENT + + DOUBLE PRECISION A,B,C,D,E,F,RB2E2,RC2F2,XSC,YSC,P1,P2,P,WS,WC, + : OR,HP,SHP,CHP,SOR,COR,DET,X0,Y0,slDA1P + + + +* Copy the six coefficients. + A = COEFFS(1) + B = COEFFS(2) + C = COEFFS(3) + D = COEFFS(4) + E = COEFFS(5) + F = COEFFS(6) + +* Scales. + RB2E2 = SQRT(B*B+E*E) + RC2F2 = SQRT(C*C+F*F) + IF (B*F-C*E.GE.0D0) THEN + XSC = RB2E2 + ELSE + B = -B + E = -E + XSC = -RB2E2 + END IF + YSC = RC2F2 + +* Non-perpendicularity. + IF (C.NE.0D0.OR.F.NE.0D0) THEN + P1 = ATAN2(C,F) + ELSE + P1 = 0D0 + END IF + IF (E.NE.0D0.OR.B.NE.0D0) THEN + P2 = ATAN2(E,B) + ELSE + P2 = 0D0 + END IF + P = slDA1P(P1+P2) + +* Orientation. + WS = C*RB2E2-E*RC2F2 + WC = B*RC2F2+F*RB2E2 + IF (WS.NE.0D0.OR.WC.NE.0D0) THEN + OR = ATAN2(WS,WC) + ELSE + OR = 0D0 + END IF + +* Zero points. + HP = P/2D0 + SHP = SIN(HP) + CHP = COS(HP) + SOR = SIN(OR) + COR = COS(OR) + DET = XSC*YSC*(CHP+SHP)*(CHP-SHP) + IF (ABS(DET).GT.0D0) THEN + X0 = YSC*(A*(CHP*COR-SHP*SOR)-D*(CHP*SOR+SHP*COR))/DET + Y0 = XSC*(A*(CHP*SOR-SHP*COR)+D*(CHP*COR+SHP*SOR))/DET + ELSE + X0 = 0D0 + Y0 = 0D0 + END IF + +* Results. + XZ = X0 + YZ = Y0 + XS = XSC + YS = YSC + PERP = P + ORIENT = OR + + END |