1 /*
2 *+
3 * Name:
4 * palDcmpf
5
6 * Purpose:
7 * Decompose an [x,y] linear fit into its constituent parameters:
8 * zero points, scales, nonperpendicularity and orientation.
9
10 * Language:
11 * Starlink ANSI C
12
13 * Type of Module:
14 * Library routine
15
16 * Invocation:
17 * palDcmpf ( double coeffs[6], double *xz, double *yz, double *xs,
18 * double *ys, double *perp, double *orient )
19
20 * Arguments:
21 * coeffs = double[6] (Given)
22 * transformation coefficients (see note)
23 * xz = double (Returned)
24 * x zero point
25 * yz = double (Returned)
26 * y zero point
27 * xs = double (Returned)
28 * x scale
29 * ys = double (Returned)
30 * y scale
31 * perp = double (Returned)
32 * nonperpendicularity (radians)
33 * orient = double (Returned)
34 * orientation (radians)
35
36 * Description:
37 * The model relates two sets of [x,y] coordinates as follows.
38 * Naming the elements of coeffs:
39 * ---
40 * coeffs[0] = A
41 * coeffs[1] = B
42 * coeffs[2] = C
43 * coeffs[3] = D
44 * coeffs[4] = E
45 * coeffs[5] = F
46 * ---
47 * the model transforms coordinates [x1,y1] into coordinates
48 * [x2,y2] as follows:
49 * ---
50 * x2 = A + B * x1 + C * y1
51 * y2 = D + E * x1 + F * y1
52 * ---
53 * The transformation can be decomposed into four steps:
54 * ---
55 * 1) Zero points:
56 * x' = xz + x1
57 * y' = yz + y1
58 *
59 * 2) Scales:
60 * x'' = xs * x'
61 * y'' = ys * y'
62 *
63 * 3) Nonperpendicularity:
64 * x''' = cos(perp / 2) * x'' + sin(perp / 2) * y''
65 * y''' = sin(perp / 2) * x'' + cos(perp / 2) * y''
66 *
67 * 4) Orientation:
68 * x2 = cos(orient) * x''' + sin(orient) * y'''
69 * y2 = -sin(orient) * y''' + cos(orient) * y'''
70 * ---
71
72 * See also:
73 * palFitxy, palPxy, palInvf and palXy2xy
74
75 * Authors:
76 * PTW: Pat Wallace (STFC)
77 * GSB: Graham Bell (EAO)
78
79 * History:
80 * 2001-12-19 (PTW):
81 * SLALIB implementation.
82 * 2018-09-13 (GSB):
83 * Initial version in C.
84
85 * Copyright:
86 * Copyright (C) 2001 Rutherford Appleton Laboratory.
87 * Copyright (C) 2018 East Asian Observatory.
88
89 * Licence:
90 * This program is free software; you can redistribute it and/or modify
91 * it under the terms of the GNU General Public License as published by
92 * the Free Software Foundation; either version 2 of the License, or
93 * (at your option) any later version.
94 *
95 * This program is distributed in the hope that it will be useful,
96 * but WITHOUT ANY WARRANTY; without even the implied warranty of
97 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
98 * GNU General Public License for more details.
99 *
100 * You should have received a copy of the GNU General Public License
101 * along with this program (see SLA_CONDITIONS); if not, write to the
102 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
103 * Boston, MA 02110-1301 USA
104 *-
105 */
106
107 #include <math.h>
108
109 #include "pal.h"
110 #include "palmac.h"
111
palDcmpf(double coeffs[6],double * xz,double * yz,double * xs,double * ys,double * perp,double * orient)112 void palDcmpf ( double coeffs[6], double *xz, double *yz, double *xs,
113 double *ys, double *perp, double *orient ) {
114
115 double a, b, c, d, e, f, rb2e2, rc2f2, xsc, ysc, p1, p2, p, ws, wc,
116 or, hp, shp, chp, sor, cor, det, x0, y0;
117
118 /* Copy the six coefficients. */
119 a = coeffs[0];
120 b = coeffs[1];
121 c = coeffs[2];
122 d = coeffs[3];
123 e = coeffs[4];
124 f = coeffs[5];
125
126 /* Scales. */
127 rb2e2 = sqrt(b*b + e*e);
128 rc2f2 = sqrt(c*c + f*f);
129 if (b*f - c*e >= 0.0) {
130 xsc = rb2e2;
131 } else {
132 b = -b;
133 e = -e;
134 xsc = -rb2e2;
135 }
136 ysc = rc2f2;
137
138 /* Non-perpendicularity. */
139 if (c != 0.0 || f != 0.0) {
140 p1 = atan2(c, f);
141 } else {
142 p1 = 0.0;
143 }
144 if (e != 0.0 || b != 0.0) {
145 p2 = atan2(e, b);
146 } else {
147 p2 = 0.0;
148 }
149 p = palDrange(p1 + p2);
150
151 /* Orientation. */
152 ws = c*rb2e2 - e*rc2f2;
153 wc = b*rc2f2 + f*rb2e2;
154 if (ws != 0.0 || wc != 0.0) {
155 or = atan2(ws, wc);
156 } else {
157 or = 0.0;
158 }
159
160 /* Zero points. */
161 hp = p / 2.0;
162 shp = sin(hp);
163 chp = cos(hp);
164 sor = sin(or);
165 cor = cos(or);
166 det = xsc * ysc * (chp + shp) * (chp - shp);
167 if (fabs(det) > 0.0) {
168 x0 = ysc * (a * (chp*cor - shp*sor) - d * (chp*sor + shp*cor)) / det;
169 y0 = xsc * (a * (chp*sor - shp*cor) + d * (chp*cor + shp*sor)) / det;
170 } else {
171 x0 = 0.0;
172 y0 = 0.0;
173 }
174
175 /* Results. */
176 *xz = x0;
177 *yz = y0;
178 *xs = xsc;
179 *ys = ysc;
180 *perp = p;
181 *orient = or;
182 }
183