/******************************************************************************* siderostat.c Source: ------- USNO/NRL Optical Interferometer 3450 Massachusetts Avenue NW Washington DC 20392-5420 Description: ------------ Implements functions related to siderostat pointing Model parameters are as described in Matt Kelly's thesis Modification history: --------------------- Written and debugged September 1, 1990. dm ported to C March 2, 1997 dm Some modifications by CAH, Mar 2001: Rename funcs to sid_funcs Aug 2001: disable "refrac", done in OYSTER *******************************************************************************/ #include #include #include #include #define LATITUDE 35.0967 #define PI 3.14159365358979323 #define RAD (PI/180.) #define NUM_SIDS 3 /* #define R0 0.0127 */ #define R0 0.01675 /* Forward function definitions */ int sid_funcs(double ha, double dec, double feedAz, double feedEl, double sidAz, double sidEl, double axisOff, double zeroAz, double mirrorOff, double zeroEl, double *az, double *el); int get_motor(double *thetaa, double *thetae, double *hastar, double *decstar, double *fazi, double *fele, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) ; int get_angles(double *azimuth, double *elevation, double *wp, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) ; int get_star(double *az, double *el, double *hastar, double *decstar, double *fazi, double *fele, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) ; int get_norm (double *az, double *el, double *wp, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0 ) ; int rotate (char *axis, double *theta, double *phi, double *in, double *out ) ; int angle_to_vector(double *azi, double *ele, double *vec) ; int vector_to_angle (double *azi, double *ele, double *vec) ; int refrac (double *s, double *sr) ; int un_refrac(double *s, double *sr) ; /* Table of constant values */ static double c_b123 = 180.; /*************************************************************** funcs For star positions, X, and model Q, calculate the motor positions, Y, and derivatives with respect to the model parameters, DYDA. There are NPAR model parameters. ***************************************************************/ int sid_funcs(double ha, double dec, double feedAz, double feedEl, double sidAz, double sidEl, double axisOff, double zeroAz, double mirrorOff, double zeroEl, double *az, double *el) { static double delta = .01; int ret_val, i__1, i__2; static double qtry[8]; static int i, j; static double y0[2], y1[2]; double q[8]; int npar = 8; q[0] = feedAz; q[1] = feedEl; q[2] = sidAz; q[3] = sidEl; q[4] = axisOff; q[5] = zeroAz; q[6] = mirrorOff; q[7] = zeroEl; /* Function Body */ if (npar != 8) { printf ( "funcs: incorrect initialization\n" ); exit ( 1 ) ; } i__1 = npar; for (i = 1; i <= i__1; ++i) { qtry[i - 1] = q[i-1]; } get_motor(az, el, &ha, &dec, qtry, &qtry[1], &qtry[2], &qtry[3] , &qtry[4], &qtry[5], &qtry[6], &qtry[7]); i__1 = npar; for (i = 1; i <= i__1; ++i) { i__2 = npar; for (j = 1; j <= i__2; ++j) { qtry[j - 1] = q[j-1]; } qtry[i - 1] = q[i-1] - delta; get_motor(y0, &y0[1], &ha, &dec, qtry, &qtry[1], &qtry[2], &qtry[3], &qtry[4], &qtry[5], &qtry[6], &qtry[7]); qtry[i - 1] = q[i-1] + delta; get_motor(y1, &y1[1], &ha, &dec, qtry, &qtry[1], &qtry[2], &qtry[3], &qtry[4], &qtry[5], &qtry[6], &qtry[7]); for (j = 1; j <= 2; ++j) { /*dyda[j + (i << 1)] = (y1[j - 1] - y0[j - 1]) / (delta * 2.);*/ } } ret_val = 0; return ret_val; } /*************************************************************** get_motor Calculates the rotation of the siderostat axes. HASTAR hour angle of star (in hours) DECSTAR declination of star THETAA angle of azimuth axis (output) THETAE angle of elevation axis (output) Fazi feed beam azimuth (0=north, 90 = west) Fele feed beam elevation (0=horizon, 90 = vertical) See comments in GET_ANGLES for definitions of other model parameters. ***************************************************************/ int get_motor(double *thetaa, double *thetae, double *hastar, double *decstar, double *fazi, double *fele, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) { double d__1; static double norm, f[3]; static int i; static double s[3], ap[3]; /* Calculate the vector toward the star, including atmospheric refraction. */ angle_to_vector(fazi, fele, f); d__1 = *hastar * -15.; angle_to_vector(&d__1, decstar, s); d__1 = 90. - LATITUDE; rotate ("Y ", &c_b123, &d__1, s, s); /* refrac(s, s); */ /* Calculate mirror normal, Ap, in local coordinates. */ norm = 0.; for (i = 1; i <= 3; ++i) { ap[i - 1] = s[i - 1] + f[i - 1]; norm += ap[i - 1] * ap[i - 1]; } norm = sqrt(norm); for (i = 1; i <= 3; ++i) { ap[i - 1] /= norm; } get_angles(thetaa, thetae, ap, sazi, sele, delta_a, azi0, delta_m, ele0); return 0; } /*************************************************************** get_angles Calculates the rotation of the siderostat axes. All angles are in degrees, except star hour angle. AZIMUTH angle of azimuth axis (output) ELEVATION angle of elevation axis (output)/ AP is the mirror normal in the local coordinate system. Siderostat model parameters. SAZI siderostat azimuth (0=north, 90 = west) SELE siderostat elevation (0=horizon, 90 = vertical) DELTA_A error in angle between axes (0=> axes are orthogonal) DELTA_M tilt of mirror in its cell (0=> mirror normal perp to axis) AZI0 zero point on azimuth axis ELE0 zero point on elevation axis With the siderostat axis roughly vertical, The motors are at azi=180, ele=90 when the normal to the mirror is pointing north and is normal to the azimuth axis. See Matt Kelly's thesis, appendix A for a description of this algorithm. ***************************************************************/ int get_angles(double *azimuth, double *elevation, double *wp, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) { double d__1, d__2; static double a, b, costa, coste, sinta, sinte, ap[3] ; /* Parameter adjustments */ --wp; /* Function Body */ d__1 = *sele - 90.; d__2 = *sazi - 90.; rotate ("ZT", &d__1, &d__2, &wp[1], ap ); coste = (ap[2] + sin(*delta_m * RAD) * sin(*delta_a * RAD)) / (cos(*delta_m * RAD) * cos(*delta_a * RAD)); sinte = sqrt(1. - coste * coste); a = cos(*delta_a * RAD) * sin(*delta_m * RAD) + sin(*delta_a * RAD) * cos(*delta_m * RAD) * coste; b = -cos(*delta_m * RAD) * sinte; sinta = a * ap[0] + b * ap[1]; costa = b * ap[0] - a * ap[1]; *elevation = atan2(sinte, coste)/RAD - *ele0; *azimuth = atan2(sinta, costa)/RAD - *azi0; if (*elevation < -180.) { *elevation += 360.; } if (*elevation > 180.) { *elevation += -360.; } if (*azimuth < -180.) { *azimuth += 360.; } if (*azimuth > 180.) { *azimuth += -360.; } return 0; } /* ----------------------------------------------------------------------- */ /* Calculate star position, HASTAR, DECSTAR, given the model */ /* and motor angles. */ int get_star(double *az, double *el, double *hastar, double *decstar, double *fazi, double *fele, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0) { double d__1; static double norm; static int i; static double s[3], fdotm; static double wf[3], wp[3]; /* Get the mirror normal */ get_norm(az, el, wp, sazi, sele, delta_a, azi0, delta_m, ele0); /* Convert mirror normal star direction. */ angle_to_vector(fazi, fele, wf); fdotm = wp[0] * wf[0] + wp[1] * wf[1] + wp[2] * wf[2]; norm = 0.; for (i = 1; i <= 3; ++i) { s[i - 1] = fdotm * 2. * wp[i - 1] - wf[i - 1]; norm += s[i - 1] * s[i - 1]; } if ((d__1 = norm - 1., abs(d__1)) > 1e-5) { printf ( "get_star: unit vector has length %f\n", norm ) ; } /* Remove atmospheric refraction. */ /* un_refrac(s, s); */ /* Convert alt and azi of star to hour angle and declination. */ d__1 = 90. - LATITUDE ; rotate ("YT", &c_b123, &d__1, s, s ); vector_to_angle(hastar, decstar, s); *hastar = -(*hastar) / 15.; return 0; } /************************************************************ get the normal to the mirror, given motor angles. ************************************************************/ int get_norm (double *az, double *el, double *wp, double *sazi, double *sele, double *delta_a, double *azi0, double *delta_m, double *ele0 ) { double d__1, d__2; static double ap[3], ep[3], mp[3]; static double cel, caz, sel, saz; /* Function Body */ caz = cos((*az + *azi0) * RAD); saz = sin((*az + *azi0) * RAD); cel = cos((*el + *ele0) * RAD); sel = sin((*el + *ele0) * RAD); /* Calculate the unit vector in the mirror frame */ mp[0] = 0.; mp[1] = 1.; mp[2] = 0.; d__1 = *el + *ele0; rotate ("Z ", delta_m, &d__1, mp, ep ); d__1 = *delta_a + 90.; d__2 = *az + *azi0; rotate ("Z ", &d__1, &d__2, ep, ap ); d__1 = *sele - 90.; d__2 = *sazi - 90.; rotate ("Z ", &d__1, &d__2, ap, &wp[0] ); return 0; } /************************************************************ Performs a coordinate transformation between two reference frames. initial frame unit vectors x0, y0, z0 final frame unit vectors x1, y1, z1 Description of transformation. if rotate by PHI rotate by THETA AXIS is about then about X X Y Y Y Z Z Z X If the axis designation is followed by an I or a T, then the transformstion is the inverse (=transpose). That is performed in the other order and with the opposite sign. The angles THETA and PHI are in degrees. Vectors IN and OUT can be the same. ****************************************/ int rotate (char *axis, double *theta, double *phi, double *in, double *out ) { static double temp[3]; static int i, j, i1, i2, i3; static double matrix[9] ; /* was [3][3] */ /* Parameter adjustments */ --out; --in; /* Function Body */ if (*axis == 'Z') { i1 = 1; i2 = 2; i3 = 3; } else if (*axis == 'X') { i1 = 2; i2 = 3; i3 = 1; } else if (*axis == 'Y') { i1 = 3; i2 = 1; i3 = 2; } else { printf ( "rotate: bad value of axis, \"%c\"\n", *axis ) ; exit(1) ; } /* ROW COL */ matrix[i1 + i1 * 3 - 4] = cos(*phi * RAD); matrix[i1 + i2 * 3 - 4] = -sin(*phi * RAD) * cos(*theta * RAD); matrix[i1 + i3 * 3 - 4] = sin(*phi * RAD) * sin(*theta * RAD); matrix[i2 + i1 * 3 - 4] = sin(*phi * RAD); matrix[i2 + i2 * 3 - 4] = cos(*phi * RAD) * cos(*theta * RAD); matrix[i2 + i3 * 3 - 4] = -cos(*phi * RAD) * sin(*theta * RAD); matrix[i3 + i1 * 3 - 4] = 0.; matrix[i3 + i2 * 3 - 4] = sin(*theta * RAD); matrix[i3 + i3 * 3 - 4] = cos(*theta * RAD); for (i = 1; i <= 3; ++i) { temp[i - 1] = in[i]; out[i] = 0.; } if (axis[1] == 'I' || axis[1] == 'T') { for (i = 1; i <= 3; ++i) { for (j = 1; j <= 3; ++j) { out[i] += matrix[j + i * 3 - 4] * temp[j - 1]; } } } else { for (i = 1; i <= 3; ++i) { for (j = 1; j <= 3; ++j) { out[i] += matrix[i + j * 3 - 4] * temp[j - 1]; } } } return 0; } /******************** angle_to_vector converts an elevation and azimuth into a unit vector pointing in that direction inputs: ELE = elevation. Angle 'above' the X-Y plane. AZI = azimuth. Angle in the X-Y plane, measured from the X-axis, counterclockwise toward the Y-axis. output VEC(1..3) = output vector ********************/ int angle_to_vector(double *azi, double *ele, double *vec) { /* Parameter adjustments */ --vec; vec[1] = cos(RAD * *ele) * cos(RAD * *azi); vec[2] = cos(RAD * *ele) * sin(RAD * *azi); vec[3] = sin(RAD * *ele); return 0; } /* This function is the inverse of ANGLE_TO_VECTOR. Given VEC, it returns ELE and AZI. */ int vector_to_angle (double *azi, double *ele, double *vec) { static double cel ; /* Parameter adjustments */ --vec; /* Function Body */ cel = sqrt(vec[1] * vec[1] + vec[2] * vec[2]); *ele = atan2(vec[3], cel) / RAD ; *azi = atan2(vec[2], vec[1]) / RAD ; return 0; } /*********************************************** refrac Correct a star vector for atmospheric refraction. The input unrefracted vector is S, the output refracted vector is SR. The vector must be in 'local' coordinates. (X, Y, Z = north, west, vertical). **********************************************/ int refrac (double *s, double *sr) { static double ele, azi; /* Parameter adjustments */ --sr; --s; /* Function Body */ vector_to_angle(&azi, &ele, &s[1]); ele += R0 / tan(ele * RAD); angle_to_vector(&azi, &ele, &sr[1]); return 0; } /*********************************************** unrefrac Remove atmospheric refraction from a star vector. The input refracted vector is Sr, the output unrefracted vector is S. The vector must be in 'local' coordinates; (X, Y, Z = north, west, vertical). **********************************************/ int un_refrac(double *s, double *sr) { static double ele, azi; /* Parameter adjustments */ --sr; --s; /* Function Body */ vector_to_angle(&azi, &ele, &sr[1]); ele -= R0 / tan(ele * RAD); angle_to_vector(&azi, &ele, &s[1]); return 0; }