Andy K / Mbed 2 deprecated SOWB

Satellite Observers Workbench. NOT yet complete, just published for forum posters to \"cherry pick\" pieces of code as requiered as an example.

Dependencies:   mbed

satapi/satapi.c@0:0a841b89d614, 2010-10-11 (annotated)

Committer:
AjK
Date:
Mon Oct 11 10:34:55 2010 +0000
Revision:
0:0a841b89d614
Totally Alpha quality as this project isn\t completed. Just publishing it as it answers many questions asked in the forums

Who changed what in which revision?

UserRevisionLine numberNew contents of line
AjK 0:0a841b89d614 1 /****************************************************************************
AjK 0:0a841b89d614 2 * Copyright 2010 Andy Kirkham, Stellar Technologies Ltd
AjK 0:0a841b89d614 3 *
AjK 0:0a841b89d614 4 * This file is part of the Satellite Observers Workbench (SOWB).
AjK 0:0a841b89d614 5 *
AjK 0:0a841b89d614 6 * SOWB is free software: you can redistribute it and/or modify
AjK 0:0a841b89d614 7 * it under the terms of the GNU General Public License as published by
AjK 0:0a841b89d614 8 * the Free Software Foundation, either version 3 of the License, or
AjK 0:0a841b89d614 9 * (at your option) any later version.
AjK 0:0a841b89d614 10 *
AjK 0:0a841b89d614 11 * SOWB is distributed in the hope that it will be useful,
AjK 0:0a841b89d614 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
AjK 0:0a841b89d614 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
AjK 0:0a841b89d614 14 * GNU General Public License for more details.
AjK 0:0a841b89d614 15 *
AjK 0:0a841b89d614 16 * You should have received a copy of the GNU General Public License
AjK 0:0a841b89d614 17 * along with SOWB. If not, see <http://www.gnu.org/licenses/>.
AjK 0:0a841b89d614 18 *
AjK 0:0a841b89d614 19 * $Id: main.cpp 5 2010-07-12 20:51:11Z ajk $
AjK 0:0a841b89d614 20 *
AjK 0:0a841b89d614 21 ***************************************************************************/
AjK 0:0a841b89d614 22
AjK 0:0a841b89d614 23 #include "sowb.h"
AjK 0:0a841b89d614 24 #include "user.h"
AjK 0:0a841b89d614 25 #include "satapi.h"
AjK 0:0a841b89d614 26 #include "utils.h"
AjK 0:0a841b89d614 27 #include "debug.h"
AjK 0:0a841b89d614 28 #include "gpio.h"
AjK 0:0a841b89d614 29 #include "osd.h"
AjK 0:0a841b89d614 30 #include "nexstar.h"
AjK 0:0a841b89d614 31 #include "utils.h"
AjK 0:0a841b89d614 32
AjK 0:0a841b89d614 33 #ifndef M_PI
AjK 0:0a841b89d614 34 #define M_PI 3.1415926535898
AjK 0:0a841b89d614 35 #endif
AjK 0:0a841b89d614 36
AjK 0:0a841b89d614 37 #define SCAN_INTERVAL 60
AjK 0:0a841b89d614 38
AjK 0:0a841b89d614 39 double satapi_aos(SAT_POS_DATA *q, bool goto_aos) {
AjK 0:0a841b89d614 40 double tsince;
AjK 0:0a841b89d614 41 char temp1[32], temp2[32];
AjK 0:0a841b89d614 42
AjK 0:0a841b89d614 43 strcpy(q->elements[0], "Lacrosse 2");
AjK 0:0a841b89d614 44 strcpy(q->elements[1], "1 21147U 91017A 10269.82093092 0.00000020 00000-0 28786-5 0 05");
AjK 0:0a841b89d614 45 strcpy(q->elements[2], "2 21147 67.9820 220.2995 0002000 244.7811 115.2189 14.76261286 07");
AjK 0:0a841b89d614 46
AjK 0:0a841b89d614 47 observer_now(q);
AjK 0:0a841b89d614 48
AjK 0:0a841b89d614 49 P22_ASSERT;
AjK 0:0a841b89d614 50
AjK 0:0a841b89d614 51 for (q->tsince = 0; q->tsince < (SCAN_INTERVAL * 90); q->tsince += SCAN_INTERVAL) {
AjK 0:0a841b89d614 52 KICK_WATCHDOG; /* We are busy! */
AjK 0:0a841b89d614 53 satallite_calculate(q);
AjK 0:0a841b89d614 54 if (q->elevation >= 10.) {
AjK 0:0a841b89d614 55 /* Above horizon viewing. Work back to AOS. */
AjK 0:0a841b89d614 56 for (tsince = q->tsince, q->tsince--; q->elevation > 10. ; tsince = q->tsince--) {
AjK 0:0a841b89d614 57 satallite_calculate(q);
AjK 0:0a841b89d614 58 if (q->elevation < 10.) {
AjK 0:0a841b89d614 59 //sprintf(temp, "%03f Q AOS El:%.1f AZ:%.1f %dKm\r\n", q->tsince, q->elevation, q->azimuth, (int)q->range);
AjK 0:0a841b89d614 60 //debug_printf(temp);
AjK 0:0a841b89d614 61 q->tsince = tsince;
AjK 0:0a841b89d614 62 satallite_calculate(q);
AjK 0:0a841b89d614 63 sprintf(temp1, "%03f T AOS El:%.1f AZ:%.1f %dKm\r\n", q->tsince, q->elevation, q->azimuth, (int)q->range);
AjK 0:0a841b89d614 64 debug_printf(temp1);
AjK 0:0a841b89d614 65 P22_DEASSERT;
AjK 0:0a841b89d614 66 if (goto_aos) {
AjK 0:0a841b89d614 67 sprintf(temp1, "%s T-%.2f", q->elements[0], tsince);
AjK 0:0a841b89d614 68 osd_string_xy(1, 12, temp1);
AjK 0:0a841b89d614 69 sprintf(temp1, "AOS %.2f%c %s%c %dKm", q->elevation, 176, printDouble_3_2(temp2, q->azimuth), 176, (int)q->range);
AjK 0:0a841b89d614 70 osd_string_xy(1, 13, temp1);
AjK 0:0a841b89d614 71 _nexstar_goto((uint32_t)((q->elevation / 360.) * 65536), (uint32_t)((q->azimuth / 360.) * 65536));
AjK 0:0a841b89d614 72 }
AjK 0:0a841b89d614 73 return tsince;
AjK 0:0a841b89d614 74 }
AjK 0:0a841b89d614 75 }
AjK 0:0a841b89d614 76 }
AjK 0:0a841b89d614 77 }
AjK 0:0a841b89d614 78
AjK 0:0a841b89d614 79 P22_DEASSERT;
AjK 0:0a841b89d614 80 return 0.;
AjK 0:0a841b89d614 81 }
AjK 0:0a841b89d614 82
AjK 0:0a841b89d614 83 int satallite_calculate(SAT_POS_DATA *q) {
AjK 0:0a841b89d614 84 double tsince;
AjK 0:0a841b89d614 85
AjK 0:0a841b89d614 86 /* Ensure the time and place are valid. */
AjK 0:0a841b89d614 87 if (!q->time.is_valid) return -1;
AjK 0:0a841b89d614 88 if (!q->location.is_valid) return -2;
AjK 0:0a841b89d614 89
AjK 0:0a841b89d614 90 ClearFlag(ALL_FLAGS);
AjK 0:0a841b89d614 91
AjK 0:0a841b89d614 92 Get_Next_Tle_Set(q->elements, &q->tle);
AjK 0:0a841b89d614 93
AjK 0:0a841b89d614 94 select_ephemeris(&q->tle);
AjK 0:0a841b89d614 95
AjK 0:0a841b89d614 96 q->jd_utc = gps_julian_date(&q->time);
AjK 0:0a841b89d614 97 q->jd_epoch = Julian_Date_of_Epoch(q->tle.epoch);
AjK 0:0a841b89d614 98
AjK 0:0a841b89d614 99 tsince = ((q->jd_utc + (q->tsince * (1 / 86400.))) - q->jd_epoch) * xmnpda;
AjK 0:0a841b89d614 100
AjK 0:0a841b89d614 101 if (isFlagSet(DEEP_SPACE_EPHEM_FLAG)) {
AjK 0:0a841b89d614 102 SDP4(tsince, &q->tle, &q->pos, &q->vel, &q->phase);
AjK 0:0a841b89d614 103 }
AjK 0:0a841b89d614 104 else {
AjK 0:0a841b89d614 105 SGP4(tsince, &q->tle, &q->pos, &q->vel, &q->phase);
AjK 0:0a841b89d614 106 }
AjK 0:0a841b89d614 107
AjK 0:0a841b89d614 108 Convert_Sat_State(&q->pos, &q->vel);
AjK 0:0a841b89d614 109 SgpMagnitude(&q->vel); // scalar magnitude, not brightness...
AjK 0:0a841b89d614 110 q->velocity = q->vel.w;
AjK 0:0a841b89d614 111
AjK 0:0a841b89d614 112 /* Populate the geodetic_t struct from data supplied. */
AjK 0:0a841b89d614 113 q->observer.lat = q->location.latitude * de2ra;
AjK 0:0a841b89d614 114 q->observer.lon = q->location.longitude * de2ra;
AjK 0:0a841b89d614 115 q->observer.alt = q->location.height / 1000.;
AjK 0:0a841b89d614 116 if (q->location.north_south == 'S') q->observer.lat *= -1.;
AjK 0:0a841b89d614 117 if (q->location.east_west == 'W') q->observer.lon *= -1.;
AjK 0:0a841b89d614 118
AjK 0:0a841b89d614 119 Calculate_Obs(q->jd_utc, &q->pos, &q->vel, &q->observer, &q->obs_set);
AjK 0:0a841b89d614 120 Calculate_LatLonAlt(q->jd_utc, &q->pos, &q->sat_geodetic);
AjK 0:0a841b89d614 121
AjK 0:0a841b89d614 122 q->azimuth = Degrees(q->obs_set.x);
AjK 0:0a841b89d614 123 q->elevation = Degrees(q->obs_set.y);
AjK 0:0a841b89d614 124 q->range = q->obs_set.z;
AjK 0:0a841b89d614 125 q->rangeRate = q->obs_set.w;
AjK 0:0a841b89d614 126 q->height = q->sat_geodetic.alt;
AjK 0:0a841b89d614 127
AjK 0:0a841b89d614 128 return 0;
AjK 0:0a841b89d614 129 }
AjK 0:0a841b89d614 130
AjK 0:0a841b89d614 131 /** observer_now
AjK 0:0a841b89d614 132 *
AjK 0:0a841b89d614 133 * Fills the data structure with the observers time and position.
AjK 0:0a841b89d614 134 *
AjK 0:0a841b89d614 135 * @param SAT_POS_DATA * A pointer to the data structure.
AjK 0:0a841b89d614 136 */
AjK 0:0a841b89d614 137 SAT_POS_DATA * observer_now(SAT_POS_DATA *q) {
AjK 0:0a841b89d614 138 gps_get_time(&q->time);
AjK 0:0a841b89d614 139 gps_get_location_average(&q->location);
AjK 0:0a841b89d614 140 return q;
AjK 0:0a841b89d614 141 }
AjK 0:0a841b89d614 142
AjK 0:0a841b89d614 143 AltAz * radec2altaz(double siderealDegrees, GPS_LOCATION_AVERAGE *location, RaDec *radec, AltAz *altaz) {
AjK 0:0a841b89d614 144 double HA, DEC, LAT, mul, altitude, azimuth;
AjK 0:0a841b89d614 145
AjK 0:0a841b89d614 146 mul = location->north_south == 'S' ? -1.0 : 1.0;
AjK 0:0a841b89d614 147
AjK 0:0a841b89d614 148 /* Convert to radians. */
AjK 0:0a841b89d614 149 HA = siderealDegrees * (M_PI / 180.0) - (radec->ra * (M_PI / 180));
AjK 0:0a841b89d614 150 DEC = radec->dec * (M_PI / 180.0);
AjK 0:0a841b89d614 151 LAT = (location->latitude * mul) * (M_PI / 180.0);
AjK 0:0a841b89d614 152
AjK 0:0a841b89d614 153 altitude = atan2(- sin(HA) * cos(DEC), cos(LAT) * sin(DEC) - sin(LAT) * cos(DEC) * cos(HA));
AjK 0:0a841b89d614 154 azimuth = asin(sin(LAT) * sin(DEC) + cos(LAT) * cos(DEC) * cos(HA));
AjK 0:0a841b89d614 155
AjK 0:0a841b89d614 156 // Convert to degrees and swing azimuth around if needed.
AjK 0:0a841b89d614 157 altaz->alt = azimuth * 180.0 / M_PI;
AjK 0:0a841b89d614 158 altaz->azm = altitude * 180.0 / M_PI;
AjK 0:0a841b89d614 159 if (altaz->azm < 0) altaz->azm += 360.0;
AjK 0:0a841b89d614 160
AjK 0:0a841b89d614 161 return altaz;
AjK 0:0a841b89d614 162 }
AjK 0:0a841b89d614 163
AjK 0:0a841b89d614 164 RaDec * altaz2radec(double siderealDegrees, GPS_LOCATION_AVERAGE *location, AltAz *altaz, RaDec *radec) {
AjK 0:0a841b89d614 165 double ALT, AZM, LAT, HA, DEC, mul;
AjK 0:0a841b89d614 166
AjK 0:0a841b89d614 167 mul = location->north_south == 'S' ? -1.0 : 1.0;
AjK 0:0a841b89d614 168
AjK 0:0a841b89d614 169 /* Convert to radians. */
AjK 0:0a841b89d614 170 LAT = (location->latitude * mul) * (M_PI / 180.0);
AjK 0:0a841b89d614 171 ALT = altaz->alt * (M_PI / 180.0);
AjK 0:0a841b89d614 172 AZM = altaz->azm * (M_PI / 180.0);
AjK 0:0a841b89d614 173
AjK 0:0a841b89d614 174 /* Calculate the declination. */
AjK 0:0a841b89d614 175 DEC = asin( ( sin(ALT) * sin(LAT) ) + ( cos(ALT) * cos(LAT) * cos(AZM) ) );
AjK 0:0a841b89d614 176 radec->dec = DEC * 180.0 / M_PI;
AjK 0:0a841b89d614 177 while (radec->dec < 0.0) radec->dec += 360.0;
AjK 0:0a841b89d614 178 while (radec->dec > 360.0) radec->dec -= 360.0;
AjK 0:0a841b89d614 179
AjK 0:0a841b89d614 180 /* Calculate the hour angle. */
AjK 0:0a841b89d614 181 HA = ( acos((sin(ALT) - sin(LAT) * sin(DEC)) / (cos(LAT) * cos(DEC)))) * 180.0 / M_PI;
AjK 0:0a841b89d614 182 if (sin(AZM) > 0.0) HA = 360.0 - HA;
AjK 0:0a841b89d614 183
AjK 0:0a841b89d614 184 /* Correct the HA for our sidereal time. */
AjK 0:0a841b89d614 185 HA = (siderealDegrees / 360.0 * 24.0) - (HA / 15.0);
AjK 0:0a841b89d614 186 if (HA < 0.0) HA += 24.0;
AjK 0:0a841b89d614 187
AjK 0:0a841b89d614 188 /* Convert the HA into degrees for the return. */
AjK 0:0a841b89d614 189 radec->ra = HA / 24.0 * 360.0;
AjK 0:0a841b89d614 190
AjK 0:0a841b89d614 191 return radec;
AjK 0:0a841b89d614 192 }
AjK 0:0a841b89d614 193