satapi.c

00001 /****************************************************************************
00002  *    Copyright 2010 Andy Kirkham, Stellar Technologies Ltd
00003  *    
00004  *    This file is part of the Satellite Observers Workbench (SOWB).
00005  *
00006  *    SOWB is free software: you can redistribute it and/or modify
00007  *    it under the terms of the GNU General Public License as published by
00008  *    the Free Software Foundation, either version 3 of the License, or
00009  *    (at your option) any later version.
00010  *
00011  *    SOWB is distributed in the hope that it will be useful,
00012  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  *    GNU General Public License for more details.
00015  *
00016  *    You should have received a copy of the GNU General Public License
00017  *    along with SOWB.  If not, see <http://www.gnu.org/licenses/>.
00018  *
00019  *    $Id: main.cpp 5 2010-07-12 20:51:11Z ajk $
00020  *    
00021  ***************************************************************************/
00022  
00023 #include "sowb.h"
00024 #include "user.h"
00025 #include "satapi.h"
00026 #include "utils.h"
00027 #include "debug.h"
00028 #include "gpio.h"
00029 #include "osd.h"
00030 #include "nexstar.h"
00031 #include "utils.h"
00032 
00033 #ifndef M_PI
00034 #define M_PI 3.1415926535898
00035 #endif
00036 
00037 #define SCAN_INTERVAL 60 
00038 
00039 double satapi_aos(SAT_POS_DATA *q, bool goto_aos) {
00040     double tsince;
00041     char temp1[32], temp2[32];
00042         
00043     strcpy(q->elements[0], "Lacrosse 2");
00044     strcpy(q->elements[1], "1 21147U 91017A   10269.82093092 0.00000020  00000-0  28786-5 0    05");
00045     strcpy(q->elements[2], "2 21147  67.9820 220.2995 0002000 244.7811 115.2189 14.76261286    07");
00046         
00047     observer_now(q);
00048     
00049     P22_ASSERT;
00050        
00051     for (q->tsince = 0; q->tsince < (SCAN_INTERVAL * 90); q->tsince += SCAN_INTERVAL) {
00052         KICK_WATCHDOG; /* We are busy! */
00053         satallite_calculate(q);
00054         if (q->elevation >= 10.) {
00055             /* Above horizon viewing. Work back to AOS. */
00056             for (tsince = q->tsince, q->tsince--; q->elevation > 10. ; tsince = q->tsince--) {
00057                 satallite_calculate(q);
00058                 if (q->elevation < 10.) {
00059                     //sprintf(temp, "%03f Q AOS El:%.1f AZ:%.1f %dKm\r\n", q->tsince, q->elevation, q->azimuth, (int)q->range);
00060                     //debug_printf(temp);
00061                     q->tsince = tsince;
00062                     satallite_calculate(q);
00063                     sprintf(temp1, "%03f T AOS El:%.1f AZ:%.1f %dKm\r\n", q->tsince, q->elevation, q->azimuth, (int)q->range);
00064                     debug_printf(temp1);
00065                     P22_DEASSERT;
00066                     if (goto_aos) {
00067                         sprintf(temp1, "%s  T-%.2f", q->elements[0], tsince);                                     
00068                         osd_string_xy(1, 12, temp1);                        
00069                         sprintf(temp1, "AOS %.2f%c %s%c %dKm", q->elevation, 176, printDouble_3_2(temp2, q->azimuth), 176, (int)q->range);
00070                         osd_string_xy(1, 13, temp1);
00071                         _nexstar_goto((uint32_t)((q->elevation / 360.) * 65536), (uint32_t)((q->azimuth / 360.) * 65536));
00072                     }
00073                     return tsince;
00074                 }
00075             }
00076         }        
00077     }
00078            
00079     P22_DEASSERT;
00080     return 0.;     
00081 }
00082 
00083 int satallite_calculate(SAT_POS_DATA *q) {
00084     double tsince;
00085 
00086     /* Ensure the time and place are valid. */
00087     if (!q->time.is_valid)      return -1;
00088     if (!q->location.is_valid)  return -2;
00089         
00090     ClearFlag(ALL_FLAGS);
00091     
00092     Get_Next_Tle_Set(q->elements, &q->tle);
00093 
00094     select_ephemeris(&q->tle);
00095     
00096     q->jd_utc = gps_julian_date(&q->time);
00097     q->jd_epoch = Julian_Date_of_Epoch(q->tle.epoch);
00098     
00099     tsince = ((q->jd_utc + (q->tsince * (1 / 86400.))) - q->jd_epoch) * xmnpda;
00100     
00101     if (isFlagSet(DEEP_SPACE_EPHEM_FLAG)) {
00102         SDP4(tsince, &q->tle, &q->pos, &q->vel, &q->phase);
00103     }
00104     else {
00105         SGP4(tsince, &q->tle, &q->pos, &q->vel, &q->phase);
00106     }
00107 
00108     Convert_Sat_State(&q->pos, &q->vel);
00109     SgpMagnitude(&q->vel); // scalar magnitude, not brightness...
00110     q->velocity = q->vel.w;
00111 
00112     /* Populate the geodetic_t struct from data supplied. */
00113     q->observer.lat   = q->location.latitude  * de2ra;
00114     q->observer.lon   = q->location.longitude * de2ra;
00115     q->observer.alt   = q->location.height / 1000.;
00116     if (q->location.north_south == 'S') q->observer.lat *= -1.;
00117     if (q->location.east_west   == 'W') q->observer.lon *= -1.;
00118     
00119     Calculate_Obs(q->jd_utc, &q->pos, &q->vel, &q->observer, &q->obs_set);
00120     Calculate_LatLonAlt(q->jd_utc, &q->pos, &q->sat_geodetic);
00121 
00122     q->azimuth   = Degrees(q->obs_set.x);
00123     q->elevation = Degrees(q->obs_set.y);
00124     q->range     = q->obs_set.z;
00125     q->rangeRate = q->obs_set.w;
00126     q->height    = q->sat_geodetic.alt;
00127         
00128     return 0;
00129 }
00130 
00131 /** observer_now
00132  *
00133  * Fills the data structure with the observers time and position.
00134  * 
00135  * @param SAT_POS_DATA * A pointer to the data structure.
00136  */
00137 SAT_POS_DATA * observer_now(SAT_POS_DATA *q) {
00138     gps_get_time(&q->time);
00139     gps_get_location_average(&q->location);
00140     return q;    
00141 }
00142 
00143 AltAz * radec2altaz(double siderealDegrees, GPS_LOCATION_AVERAGE *location, RaDec *radec, AltAz *altaz) {
00144     double HA, DEC, LAT, mul, altitude, azimuth;
00145     
00146     mul = location->north_south == 'S' ? -1.0 : 1.0;
00147      
00148     /* Convert to radians. */
00149     HA = siderealDegrees * (M_PI / 180.0) - (radec->ra * (M_PI / 180));
00150     DEC = radec->dec * (M_PI / 180.0);
00151     LAT = (location->latitude * mul) * (M_PI / 180.0);
00152     
00153     altitude = atan2(- sin(HA) * cos(DEC), cos(LAT) * sin(DEC) - sin(LAT) * cos(DEC) * cos(HA));
00154     azimuth = asin(sin(LAT) * sin(DEC) + cos(LAT) * cos(DEC) * cos(HA));
00155 
00156     // Convert to degrees and swing azimuth around if needed.
00157     altaz->alt = azimuth * 180.0 / M_PI;
00158     altaz->azm = altitude * 180.0 / M_PI;
00159     if (altaz->azm < 0) altaz->azm += 360.0;
00160   
00161     return altaz;
00162 }
00163 
00164 RaDec * altaz2radec(double siderealDegrees, GPS_LOCATION_AVERAGE *location, AltAz *altaz, RaDec *radec) {
00165     double ALT, AZM, LAT, HA, DEC, mul;
00166     
00167     mul = location->north_south == 'S' ? -1.0 : 1.0;
00168     
00169     /* Convert to radians. */
00170     LAT = (location->latitude * mul) * (M_PI / 180.0);
00171     ALT = altaz->alt * (M_PI / 180.0);
00172     AZM = altaz->azm * (M_PI / 180.0);
00173     
00174     /* Calculate the declination. */
00175     DEC = asin( ( sin(ALT) * sin(LAT) ) + ( cos(ALT) * cos(LAT) * cos(AZM) ) );
00176     radec->dec = DEC * 180.0 / M_PI;
00177     while (radec->dec < 0.0)   radec->dec += 360.0;
00178     while (radec->dec > 360.0) radec->dec -= 360.0;
00179     
00180     /* Calculate the hour angle. */
00181     HA = ( acos((sin(ALT) - sin(LAT) * sin(DEC)) / (cos(LAT) * cos(DEC)))) * 180.0 / M_PI;
00182     if (sin(AZM) > 0.0) HA = 360.0 - HA;
00183     
00184     /* Correct the HA for our sidereal time. */    
00185     HA = (siderealDegrees / 360.0 * 24.0) - (HA / 15.0);
00186     if (HA < 0.0) HA += 24.0;
00187     
00188     /* Convert the HA into degrees for the return. */
00189     radec->ra = HA / 24.0 * 360.0;
00190     
00191     return radec;
00192 }
00193