plan9port

[fork] Plan 9 from user space
git clone git://src.adamsgaard.dk/plan9port # fast
git clone https://src.adamsgaard.dk/plan9port.git # slow
Log | Files | Refs | README | LICENSE Back to index

posn.c (5532B)


      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 }