Parent Directory
| 
Revision Log
Parent Directory
| Revision Log
New gps integration
| 1 | // http://answers.google.com/answers/threadview/id/782886.html |
| 2 | // gcc `pkg-config --cflags --libs glib-2.0` -lm -o astro astro.c |
| 3 | |
| 4 | #include <stdio.h> |
| 5 | #include <glib.h> |
| 6 | #include <math.h> |
| 7 | |
| 8 | #define DEG2RAD(a) ((a)/180.0*M_PI) |
| 9 | #define RAD2DEG(a) ((a)/M_PI*180.0) |
| 10 | |
| 11 | int main(void) { |
| 12 | puts("astro"); |
| 13 | |
| 14 | #if 0 |
| 15 | // For exemple, right now, I am at 48°49'N 02°17'23 E, and we are the |
| 16 | // 15th of november 2006, time is 1035 GMT. |
| 17 | guint day = 15, month = 11, year = 2006; |
| 18 | guint hour = 10, minute = 35; |
| 19 | gdouble lon = 02.0 + 17.0/60 + 23.0/3600; |
| 20 | gdouble lat = 48.0 + 49.0/60; |
| 21 | #else |
| 22 | // 26.5.2010 |
| 23 | // 11:50 GMT |
| 24 | guint day = 12, month = 1, year = 2001; |
| 25 | guint hour = 11, minute = 0; |
| 26 | gdouble lon = -48; |
| 27 | gdouble lat = -49; |
| 28 | #endif |
| 29 | |
| 30 | for(hour = 2;hour < 23;hour++) { |
| 31 | printf("%u: ", hour); |
| 32 | |
| 33 | GDate *now = g_date_new_dmy(day, month, year); |
| 34 | guint doy = g_date_get_day_of_year(now); |
| 35 | // printf("day of year = %u\n", doy); |
| 36 | |
| 37 | gdouble hourfrac = hour+(minute/60.0); |
| 38 | |
| 39 | double g = (360./365.25)*(doy + hourfrac/24); |
| 40 | // printf("g = %f\n", g); |
| 41 | |
| 42 | g = DEG2RAD(g); |
| 43 | double D = 0.396372-22.91327*cos(g)+4.02543*sin(g)-0.387205*cos(2*g)+ |
| 44 | +0.051967*sin(2*g)-0.154527*cos(3*g) + 0.084798*sin(3*g); |
| 45 | |
| 46 | // printf("D = %f\n", D); |
| 47 | |
| 48 | double TC = 0.004297+0.107029*cos(g)-1.837877*sin(g)-0.837378*cos(2*g) |
| 49 | -2.340475*sin(2*g); |
| 50 | |
| 51 | // printf("TC = %f\n", TC); |
| 52 | |
| 53 | double SHA = (hourfrac-12)*15 + lon + TC; |
| 54 | |
| 55 | // printf("SHA = %f\n", SHA); |
| 56 | |
| 57 | // printf("lat = %f\n", lat); |
| 58 | |
| 59 | double SZA = acos(sin(DEG2RAD(lat))*sin(DEG2RAD(D))+ |
| 60 | cos(DEG2RAD(lat))*cos(DEG2RAD(D))*cos(DEG2RAD(SHA))); |
| 61 | |
| 62 | // printf("SZA = %f\n", RAD2DEG(SZA)); |
| 63 | |
| 64 | printf("altitude = %f ", 90.0 - RAD2DEG(SZA)); |
| 65 | |
| 66 | |
| 67 | // cos(AZ) = (sin(D)-sin(Latitude)*cos(SZA))/(cos(Latitude)*sin(SZA)) = |
| 68 | |
| 69 | double AZ = acos((sin(DEG2RAD(D))-sin(DEG2RAD(lat))*cos(SZA))/(cos(DEG2RAD(lat))*sin(SZA))); |
| 70 | |
| 71 | if(hourfrac >= 12.000) |
| 72 | AZ = 2*M_PI - AZ; |
| 73 | |
| 74 | printf("AZ = %f\n", RAD2DEG(AZ)); |
| 75 | |
| 76 | } |
| 77 | |
| 78 | return 0; |
| 79 | } |