Warning: Attempt to read property "date" on null in /usr/local/www/websvn.planix.org/blame.php on line 247

Warning: Attempt to read property "msg" on null in /usr/local/www/websvn.planix.org/blame.php on line 247
WebSVN – planix.SVN – Blame – /os/branches/feature_fixcpp/sys/src/cmd/scat/posn.c – Rev 2

Subversion Repositories planix.SVN

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include	<u.h>
2 
#include	<libc.h>
3 
#include	<bio.h>
4 
#include	"sky.h"
5 
 
6 
void
7 
amdinv(Header *h, Angle ra, Angle dec, float mag, float col)
8 
{
9 
	int i, max_iterations;
10 
	float tolerance;
11 
	float object_x, object_y, delta_x, delta_y, f, fx, fy, g, gx, gy;
12 
	float x1, x2, x3, x4;
13 
	float y1, y2, y3, y4;
14 
 
15 
	/*
16 
	 *  Initialize
17 
	 */
18 
	max_iterations = 50;
19 
	tolerance = 0.0000005;
20 
 
21 
	/*
22 
	 *  Convert RA and Dec to St.coords
23 
	 */
24 
	traneqstd(h, ra, dec);
25 
 
26 
	/*
27 
	 *  Set initial value for x,y
28 
	 */
29 
	object_x = h->xi/h->param[Ppltscale];
30 
	object_y = h->eta/h->param[Ppltscale];
31 
 
32 
	/*
33 
	 *  Iterate by Newtons method
34 
	 */
35 
	for(i = 0; i < max_iterations; i++) {
36 
		/*
37 
		 *  X plate model
38 
		 */
39 
		x1 = object_x;
40 
		x2 = x1 * object_x;
41 
		x3 = x1 * object_x;
42 
		x4 = x1 * object_x;
43 
 
44 
		y1 = object_y;
45 
		y2 = y1 * object_y;
46 
		y3 = y1 * object_y;
47 
		y4 = y1 * object_y;
48 
 
49 
		f =
50 
			h->param[Pamdx1] * x1 +
51 
			h->param[Pamdx2] * y1 +
52 
		 	h->param[Pamdx3] +
53 
			h->param[Pamdx4] * x2 +
54 
			h->param[Pamdx5] * x1*y1 +
55 
			h->param[Pamdx6] * y2 +
56 
		   	h->param[Pamdx7] * (x2+y2) +
57 
			h->param[Pamdx8] * x2*x1 +
58 
			h->param[Pamdx9] * x2*y1 +
59 
			h->param[Pamdx10] * x1*y2 +
60 
			h->param[Pamdx11] * y3 +
61 
			h->param[Pamdx12] * x1* (x2+y2) +
62 
			h->param[Pamdx13] * x1 * (x2+y2) * (x2+y2) +
63 
			h->param[Pamdx14] * mag +
64 
			h->param[Pamdx15] * mag*mag +
65 
			h->param[Pamdx16] * mag*mag*mag +
66 
			h->param[Pamdx17] * mag*x1 +
67 
			h->param[Pamdx18] * mag * (x2+y2) +
68 
			h->param[Pamdx19] * mag*x1 * (x2+y2) +
69 
			h->param[Pamdx20] * col;
70 
 
71 
		/*
72 
		 *  Derivative of X model wrt x
73 
		 */
74 
		fx =
75 
			h->param[Pamdx1] +
76 
			h->param[Pamdx4] * 2*x1 +
77 
			h->param[Pamdx5] * y1 +
78 
			h->param[Pamdx7] * 2*x1 +
79 
			h->param[Pamdx8] * 3*x2 +
80 
			h->param[Pamdx9] * 2*x1*y1 +
81 
			h->param[Pamdx10] * y2 +
82 
			h->param[Pamdx12] * (3*x2+y2) +
83 
			h->param[Pamdx13] * (5*x4 + 6*x2*y2 + y4) +
84 
			h->param[Pamdx17] * mag +
85 
			h->param[Pamdx18] * mag*2*x1 +
86 
			h->param[Pamdx19] * mag*(3*x2+y2);
87 
 
88 
		/*
89 
		 *  Derivative of X model wrt y
90 
		 */
91 
		fy =
92 
			h->param[Pamdx2] +
93 
			h->param[Pamdx5] * x1 +
94 
			h->param[Pamdx6] * 2*y1 +
95 
			h->param[Pamdx7] * 2*y1 +
96 
			h->param[Pamdx9] * x2 +
97 
			h->param[Pamdx10] * x1*2*y1 +
98 
			h->param[Pamdx11] * 3*y2 +
99 
			h->param[Pamdx12] * 2*x1*y1 +
100 
			h->param[Pamdx13] * 4*x1*y1* (x2+y2) +
101 
			h->param[Pamdx18] * mag*2*y1 +
102 
			h->param[Pamdx19] * mag*2*x1*y1;
103 
		/*
104 
		 *  Y plate model
105 
		 */
106 
		g =
107 
			h->param[Pamdy1] * y1 +
108 
			h->param[Pamdy2] * x1 +
109 
			h->param[Pamdy3] +
110 
			h->param[Pamdy4] * y2 +
111 
			h->param[Pamdy5] * y1*x1 +
112 
			h->param[Pamdy6] * x2 +
113 
			h->param[Pamdy7] * (x2+y2) +
114 
			h->param[Pamdy8] * y3 +
115 
			h->param[Pamdy9] * y2*x1 +
116 
			h->param[Pamdy10] * y1*x3 +
117 
			h->param[Pamdy11] * x3 +
118 
			h->param[Pamdy12] * y1 * (x2+y2) +
119 
			h->param[Pamdy13] * y1 * (x2+y2) * (x2+y2) +
120 
			h->param[Pamdy14] * mag +
121 
			h->param[Pamdy15] * mag*mag +
122 
			h->param[Pamdy16] * mag*mag*mag +
123 
			h->param[Pamdy17] * mag*y1 +
124 
			h->param[Pamdy18] * mag * (x2+y2) +
125 
			h->param[Pamdy19] * mag*y1 * (x2+y2) +
126 
			h->param[Pamdy20] * col;
127 
 
128 
		/*
129 
		 *  Derivative of Y model wrt x
130 
		 */
131 
		gx =
132 
			h->param[Pamdy2] +
133 
			h->param[Pamdy5] * y1 +
134 
			h->param[Pamdy6] * 2*x1 +
135 
			h->param[Pamdy7] * 2*x1 +
136 
			h->param[Pamdy9] * y2 +
137 
			h->param[Pamdy10] * y1*2*x1 +
138 
			h->param[Pamdy11] * 3*x2 +
139 
			h->param[Pamdy12] * 2*x1*y1 +
140 
			h->param[Pamdy13] * 4*x1*y1 * (x2+y2) +
141 
			h->param[Pamdy18] * mag*2*x1 +
142 
			h->param[Pamdy19] * mag*y1*2*x1;
143 
 
144 
		/*
145 
		 *  Derivative of Y model wrt y
146 
		 */
147 
		gy =
148 
			h->param[Pamdy1] +
149 
			h->param[Pamdy4] * 2*y1 +
150 
			h->param[Pamdy5] * x1 +
151 
			h->param[Pamdy7] * 2*y1 +
152 
			h->param[Pamdy8] * 3*y2 +
153 
			h->param[Pamdy9] * 2*y1*x1 +
154 
			h->param[Pamdy10] * x2 +
155 
			h->param[Pamdy12] * 3*y2 +
156 
			h->param[Pamdy13] * (5*y4 + 6*x2*y2 + x4) +
157 
			h->param[Pamdy17] * mag +
158 
			h->param[Pamdy18] * mag*2*y1 +
159 
			h->param[Pamdy19] * mag*(x2 + 3*y2);
160 
 
161 
		f = f - h->xi;
162 
		g = g - h->eta;
163 
		delta_x = (-f*gy+g*fy) / (fx*gy-fy*gx);
164 
		delta_y = (-g*fx+f*gx) / (fx*gy-fy*gx);
165 
		object_x = object_x + delta_x;
166 
		object_y = object_y + delta_y;
167 
		if((fabs(delta_x) < tolerance) && (fabs(delta_y) < tolerance))
168 
			break;
169 
	}
170 
 
171 
	/*
172 
	 *  Convert mm from plate center to pixels
173 
	 */
174 
	h->x = (h->param[Pppo3]-object_x*1000.0)/h->param[Pxpixelsz];
175 
	h->y = (h->param[Pppo6]+object_y*1000.0)/h->param[Pypixelsz];
176 
}
177 
 
178 
void
179 
ppoinv(Header *h, Angle ra, Angle dec)
180 
{
181 
 
182 
	/*
183 
	 *  Convert RA and Dec to standard coords.
184 
	 */
185 
	traneqstd(h, ra, dec);
186 
 
187 
	/*
188 
	 *  Convert st.coords from arcsec to radians
189 
	 */
190 
	h->xi  /= ARCSECONDS_PER_RADIAN;
191 
	h->eta /= ARCSECONDS_PER_RADIAN;
192 
 
193 
	/*
194 
	 *  Compute PDS coordinates from solution
195 
	 */
196 
	h->x =
197 
		h->param[Pppo1]*h->xi +
198 
		h->param[Pppo2]*h->eta +
199 
		h->param[Pppo3];
200 
	h->y =
201 
		h->param[Pppo4]*h->xi +
202 
		h->param[Pppo5]*h->eta +
203 
		h->param[Pppo6];
204 
	/*
205 
	 * Convert x,y from microns to pixels
206 
	 */
207 
	h->x /= h->param[Pxpixelsz];
208 
	h->y /= h->param[Pypixelsz];
209 
 
210 
}
211 
 
212 
void
213 
traneqstd(Header *h, Angle object_ra, Angle object_dec)
214 
{
215 
	float div;
216 
 
217 
	/*
218 
	 *  Find divisor
219 
	 */
220 
	div =
221 
		(sin(object_dec)*sin(h->param[Ppltdec]) +
222 
		cos(object_dec)*cos(h->param[Ppltdec]) *
223 
		cos(object_ra - h->param[Ppltra]));
224 
 
225 
	/*
226 
	 *  Compute standard coords and convert to arcsec
227 
	 */
228 
	h->xi = cos(object_dec) *
229 
		sin(object_ra - h->param[Ppltra]) *
230 
		ARCSECONDS_PER_RADIAN/div;
231 
 
232 
	h->eta = (sin(object_dec)*cos(h->param[Ppltdec])-
233 
		cos(object_dec)*sin(h->param[Ppltdec])*
234 
		cos(object_ra - h->param[Ppltra]))*
235 
		ARCSECONDS_PER_RADIAN/div;
236 
 
237 
}
238 
 
239 
void
240 
xypos(Header *h, Angle ra, Angle dec, float mag, float col)
241 
{
242 
	if (h->amdflag) {
243 
		amdinv(h, ra, dec, mag, col);
244 
	} else {
245 
		ppoinv(h, ra, dec);
246 
	}
247 
}