home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / XBBS7200.ZIP / XBBS7200.TAR / today / sun.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-10-02  |  10.8 KB  |  551 lines
  1. /***** hpfcla:net.sources / nsc-pdc!rgb / 10:24 am  May 16, 1985
  2. *
  3. * Changed constants to Fort Collins, Colorado.  (ajs, 850520)
  4. * Made other minor output format improvements also.
  5. * 
  6. *        sun <options>
  7. *
  8. *        options:        -t hh:mm:ss    time (default is current system time)
  9. *             -d mm/dd/yy    date (default is current system date)
  10. *                        -a lat        decimal latitude (default = 45.5333)
  11. *                        -o lon        decimal longitude (default = 122.8333) 
  12. *             -z tz        timezone (default = 8, pst)
  13. *             -p        show position of sun (azimuth)
  14. *             -v        turn on debugging
  15. *        
  16. *        All output is to standard io.  
  17. *
  18. *     Compile with cc -O -o sun sun.c -lm
  19. *     Non 4.2 systems may have to change <sys/time.h> to <time.h> below.
  20. *    (yes, done)
  21. *
  22. *     Note that the latitude, longitude, time zone correction and
  23. *     time zone string are all defaulted in the global variable section.
  24. *
  25. *     Most of the code in this program is adapted from algorithms
  26. *     presented in "Practical Astronomy With Your Calculator" by
  27. *     Peter Duffet-Smith.
  28. *
  29. *     The GST and ALT-AZIMUTH algorithms are from Sky and Telescope,
  30. *     June, 1984 by Roger W. Sinnott
  31. *
  32. *     Author Robert Bond - Beaverton Oregon.
  33. *    
  34. */
  35.  
  36. #include <stdio.h>
  37. #include <math.h>
  38. #include <sys/types.h>
  39. #include <time.h>
  40.  
  41. #define PI       3.141592654
  42. #define EPOCH     1980
  43. #define JDE     2444238.5    /* Julian date of EPOCH */
  44.  
  45. static double dtor();
  46. static double adj360();
  47. double adj24();
  48. double julian_date();
  49. double hms_to_dh();
  50. double solar_lon();
  51. double acos_deg();
  52. double asin_deg();
  53. double atan_q_deg();
  54. double atan_deg();
  55. double sin_deg();
  56. double cos_deg();
  57. double tan_deg();
  58. double gmst();
  59.  
  60. long time();
  61. struct tm *localtime();
  62.  
  63. int th;
  64. int tm;
  65. int ts;
  66. int mo;
  67. int day;
  68. int yr;
  69. int tz=8;            /* Default time zone */
  70. char *tzs  = "(PST)";        /* Default time zone string */
  71. char *dtzs = "(PDT)";        /* Default daylight savings time string */
  72. int debug = 0;
  73. int popt = 0;
  74.  
  75. double lat =  35.0000;        /* Default latitude (Fort Collins, Colorado) */
  76. double lon = 118.0000;        /* Default Longitude (Degrees west) */ 
  77.  
  78. sun(sunrh, sunrm, sunsh, sunsm)
  79. int *sunrh, *sunrm, *sunsh, *sunsm;
  80. {
  81.     double ed, jd;
  82.     double alpha1, delta1, alpha2, delta2, st1r, st1s, st2r, st2s;
  83.     double a1r, a1s, a2r, a2s, dt, dh, x, y;
  84.     double trise, tset, ar, as, alpha, delta, tri, da;
  85.     double lambda1, lambda2;
  86.     double alt, az, gst, m1;
  87.     double hsm, ratio;
  88.     time_t sec_1970;
  89.     int h, m;
  90.     struct tm *pt;
  91.  
  92.     time(&sec_1970);
  93.     pt = localtime(&sec_1970);  
  94.  
  95.     th = pt->tm_hour;
  96.     tm = pt->tm_min;
  97.     ts = pt->tm_sec;
  98.     yr = pt->tm_year + 1900;
  99.     mo = pt->tm_mon + 1;
  100.     day = pt->tm_mday;
  101.     if (pt->tm_isdst) {        /* convert tz to daylight savings time */
  102.     tz--;
  103.     tzs = dtzs;    
  104.     }
  105.  
  106.  
  107.     if (debug)
  108.         printf("Date: %d/%d/%d,  Time: %d:%d:%d, Tz: %d, Lat: %lf, Lon: %lf \n",
  109.         mo,day,yr,th,tm,ts,tz,lat,lon);
  110.  
  111.     jd = julian_date(mo,day,yr);
  112.     ed = jd - JDE;
  113.  
  114.     lambda1 = solar_lon(ed);
  115.     lambda2 = solar_lon(ed + 1.0);
  116.  
  117.     lon_to_eq(lambda1, &alpha1, &delta1);
  118.     lon_to_eq(lambda2, &alpha2, &delta2);
  119.  
  120.     rise_set(alpha1, delta1, &st1r, &st1s, &a1r, &a1s);
  121.     rise_set(alpha2, delta2, &st2r, &st2s, &a2r, &a2s);
  122.  
  123.     m1 = adj24(gmst(jd - 0.5, 0.5 + tz / 24.0) - lon / 15); /* lst midnight */
  124.  
  125.     if (debug)
  126.     printf ("local sidereal time of midnight is %lf \n", m1);
  127.  
  128.     hsm = adj24(st1r - m1);
  129.  
  130.     if (debug)
  131.     printf ("about %lf hours from midnight to dawn \n", hsm);
  132.  
  133.     ratio = hsm / 24.07;
  134.  
  135.     if (debug)
  136.     printf("%lf is how far dawn is into the day \n", ratio);
  137.  
  138.     if (fabs(st2r - st1r) > 1.0) {
  139.     st2r += 24.0;
  140.     if (debug)
  141.         printf("st2r corrected from %lf to %lf \n", st2r-24.0, st2r);
  142.     }
  143.  
  144.     trise = adj24((1.0 - ratio) * st1r + ratio * st2r);
  145.  
  146.     hsm = adj24(st1s - m1);
  147.  
  148.     if (debug)
  149.     printf ("about %lf hours from midnight to sunset \n", hsm);
  150.  
  151.     ratio = hsm / 24.07;
  152.  
  153.     if (debug)
  154.     printf("%lf is how far sunset is into the day \n", ratio);
  155.  
  156.     if (fabs(st2s - st1s) > 1.0) {
  157.     st2s += 24.0;
  158.     if (debug)
  159.         printf("st2s corrected from %lf to %lf \n", st2s-24.0, st2s);
  160.     }
  161.  
  162.     tset = adj24((1.0 - ratio) * st1s + ratio * st2s);
  163.  
  164.     if (debug)
  165.     printf("Uncorrected rise = %lf, set = %lf \n", trise, tset);
  166.  
  167.     ar = a1r * 360.0 / (360.0 + a1r - a2r);
  168.     as = a1s * 360.0 / (360.0 + a1s - a2s);
  169.  
  170.     delta = (delta1 + delta2) / 2.0;
  171.     tri = acos_deg(sin_deg(lat)/cos_deg(delta));
  172.  
  173.     x = 0.835608;        /* correction for refraction, parallax, ? */
  174.     y = asin_deg(sin_deg(x)/sin_deg(tri));
  175.     da = asin_deg(tan_deg(x)/tan_deg(tri));
  176.     dt = 240.0 * y / cos_deg(delta) / 3600;
  177.  
  178.     if (debug)
  179.     printf("Corrections: dt = %lf, da = %lf \n", dt, da);
  180.  
  181.     lst_to_hm(trise - dt, jd, &h, &m);
  182.     *sunrh = h;
  183.     *sunrm = m;
  184.  
  185.     if (popt) {
  186.         dh_to_hm(ar - da, &h, &m);
  187.         printf("Azimuth: %3d %02d'\n", h, m);
  188.     }
  189.  
  190.     lst_to_hm(tset + dt, jd, &h, &m);
  191.     *sunsh = h;
  192.     *sunsm = m;
  193.  
  194.     if (popt) {
  195.         dh_to_hm(as + da, &h, &m);
  196.         printf("Azimuth: %3d %02d'\n", h, m);
  197.     } 
  198.      
  199.  
  200.     if (popt) {
  201.  
  202.     if (alpha1 < alpha2)
  203.         alpha = (alpha1 + alpha2) / 2.0;
  204.     else
  205.         alpha = (alpha1 + 24.0 + alpha2) / 2.0;
  206.     
  207.     if (alpha > 24.0)
  208.         alpha -= 24.0;
  209.  
  210.     dh = (hms_to_dh(th, tm, ts) + tz) / 24.0;
  211.     if (dh > 0.5) {
  212.         dh -= 0.5;
  213.         jd += 0.5;
  214.     } else {
  215.         dh += 0.5;
  216.         jd -= 0.5;
  217.     }
  218.  
  219.     gst = gmst(jd, dh);
  220.  
  221.     eq_to_altaz(alpha, delta, gst, &alt, &az);
  222.  
  223.     printf     ("The sun is at:   ");
  224.     dh_to_hm (az, &h, &m);
  225.     printf     ("Azimuth: %3d %02d'  ", h, m);
  226.     dh_to_hm (alt, &h, &m);
  227.     printf     ("Altitude: %3d %02d'\n", h, m);
  228.     }
  229. }
  230.  
  231. static double
  232. dtor(deg)
  233. double deg;
  234. {
  235.     return (deg * PI / 180.0);
  236. }
  237.  
  238. double
  239. rtod(deg)
  240. double deg;
  241. {
  242.     return (deg * 180.0 / PI);
  243. }
  244.  
  245.  
  246. static double 
  247. adj360(deg)
  248. double deg;
  249. {
  250.     while (deg < 0.0) 
  251.     deg += 360.0;
  252.     while (deg > 360.0)
  253.     deg -= 360.0;
  254.     return(deg);
  255. }
  256.  
  257. double 
  258. adj24(hrs)
  259. double hrs;
  260. {
  261.     while (hrs < 0.0) 
  262.     hrs += 24.0;
  263.     while (hrs > 24.0)
  264.     hrs -= 24.0;
  265.     return(hrs);
  266. }
  267.  
  268. double 
  269. julian_date(m, d, y) int m, d, y;
  270. {
  271.     long a, b;
  272.     double jd;
  273.  
  274.     if (m == 1 || m == 2) {
  275.     --y;
  276.     m += 12;
  277.     }
  278.     if (y < 1583) {
  279.     printf("Can't handle dates before 1583\n");
  280.     exit(1);
  281.     }
  282.     a = (long)y/100;
  283.     b = 2 - a + a/4;
  284.     b += (long)((double)y * 365.25);
  285.     b += (long)(30.6001 * ((double)m + 1.0));
  286.     jd = (double)d + (double)b + 1720994.5;
  287.  
  288.     if (debug) 
  289.     printf("Julian date for %d/%d/%d is %lf \n", m, d, y, jd);
  290.  
  291.     return(jd);
  292. }
  293.  
  294. double 
  295. hms_to_dh(h, m, s) int h, m, s;
  296. {
  297.     double rv;
  298.     rv = h + m / 60.0 + s / 3600.0;
  299.  
  300.     if (debug)
  301.     printf("For time %d:%d:%d frac hours are: %lf \n", h, m, s, rv);
  302.  
  303.     return rv;
  304. }
  305.  
  306. double 
  307. solar_lon(ed)
  308. double ed;
  309. {
  310.     double n, m, e, ect, errt, v;
  311.  
  312.     n = 360.0 * ed / 365.2422;
  313.     n = adj360(n);
  314.     m = n + 278.83354 - 282.596403;
  315.     m = adj360(m);
  316.     m = dtor(m);
  317.     e = m; ect = 0.016718;
  318.     while ((errt = e - ect * sin(e) - m) > 0.0000001) 
  319.         e = e - errt / (1 - ect * cos(e));
  320.     v = 2 * atan(1.0168601 * tan(e/2));
  321.     v = adj360(v * 180.0 / PI + 282.596403);
  322.  
  323.     if (debug)
  324.     printf("Solar Longitude for %lf days is %lf \n", ed, v); 
  325.  
  326.     return(v);
  327. }
  328.  
  329. double 
  330. acos_deg(x)
  331. double x;
  332. {
  333.     return rtod(acos(x));
  334. }
  335.  
  336. double 
  337. asin_deg(x)
  338. double x;
  339. {
  340.     return rtod(asin(x));
  341. }
  342.  
  343. double 
  344. atan_q_deg(y,x)
  345. double y,x;
  346. {
  347.     double rv;
  348.  
  349.     if (y == 0)
  350.         rv = 0;
  351.     else if (x == 0)
  352.         rv = y>0 ? 90.0 : -90.0;
  353.     else rv = atan_deg(y/x);
  354.  
  355.     if (x<0) return rv+180.0;
  356.     if (y<0) return rv+360.0;
  357.     return(rv);
  358. }
  359.  
  360. double
  361. atan_deg(x)
  362. double x;
  363. {
  364.     return rtod(atan(x));
  365. }
  366.  
  367. double 
  368. sin_deg(x)
  369. double x;
  370. {
  371.     return sin(dtor(x));
  372. }
  373.  
  374. double 
  375. cos_deg(x)
  376. double x;
  377. {
  378.     return cos(dtor(x));
  379. }
  380.  
  381. double 
  382. tan_deg(x)
  383. double x;
  384. {
  385.     return tan(dtor(x));
  386. }
  387.  
  388. lon_to_eq(lambda, alpha, delta)
  389. double lambda;
  390. double *alpha;
  391. double *delta;
  392. {
  393.     double tlam,epsilon;
  394.  
  395.     tlam = dtor(lambda);
  396.     epsilon = dtor((double)23.441884);
  397.     *alpha = atan_q_deg((sin(tlam))*cos(epsilon),cos(tlam)) / 15.0;
  398.     *delta = asin_deg(sin(epsilon)*sin(tlam));
  399.  
  400.     if (debug)
  401.     printf("Right ascension, declination for lon %lf is %lf, %lf \n",
  402.         lambda, *alpha, *delta);
  403. }
  404.  
  405. rise_set(alpha, delta, lstr, lsts, ar, as)
  406. double alpha, delta, *lstr, *lsts, *ar, *as;
  407. {
  408.     double tar;
  409.     double h;
  410.  
  411.     tar = sin_deg(delta)/cos_deg(lat);
  412.     if (tar < -1.0 || tar > 1.0) {
  413.     printf("The object is circumpolar\n");
  414.     exit (1);
  415.     }
  416.     *ar = acos_deg(tar);
  417.     *as = 360.0 - *ar;
  418.  
  419.     h = acos_deg(-tan_deg(lat) * tan_deg(delta)) / 15.0;
  420.     *lstr = 24.0 + alpha - h;
  421.     if (*lstr > 24.0)
  422.     *lstr -= 24.0;
  423.     *lsts = alpha + h;
  424.     if (*lsts > 24.0)
  425.     *lsts -= 24.0;
  426.  
  427.     if (debug) {
  428.     printf("For ra, decl. of %lf, %lf: \n", alpha, delta);
  429.     printf("lstr = %lf, lsts = %lf, \n", *lstr, *lsts);
  430.     printf("ar =   %lf, as =   %lf \n", *ar, *as);
  431.     }
  432. }
  433.  
  434. lst_to_hm(lst, jd, h, m)
  435. double lst, jd;
  436. int *h, *m;
  437. {
  438.     double ed, gst, jzjd, t, r, b, t0, gmt;
  439.  
  440.     gst = lst + lon / 15.0;
  441.     if (gst > 24.0)
  442.     gst -= 24.0;
  443.     jzjd = julian_date(1,0,yr);
  444.     ed = jd-jzjd;
  445.     t = (jzjd -2415020.0)/36525.0;
  446.     r = 6.6460656+2400.05126*t+2.58E-05*t*t;
  447.     b = 24.0-(r-24.0*(yr-1900));
  448.     t0 = ed * 0.0657098 - b;
  449.     if (t0 < 0.0)
  450.     t0 += 24;
  451.     gmt = gst-t0;
  452.     if (gmt<0) 
  453.     gmt += 24.0;
  454.     gmt = gmt * 0.99727 - tz;;
  455.     if (gmt < 0)
  456.     gmt +=24.0;
  457.     dh_to_hm(gmt, h, m);
  458. }
  459.  
  460. dh_to_hm(dh, h, m)
  461. double dh;
  462. int *h, *m;
  463. {
  464.     double tempsec;
  465.  
  466.     *h = dh;
  467.  /* *m = (dh - *h) * 60; 
  468.     tempsec = (dh - *h) * 60 - *m; */
  469.     *m = fmod(dh, 1.0) * 60.0; 
  470.     tempsec = fmod(dh, 1.0) * 60.0 - *m;
  471.     tempsec = tempsec * 60 + 0.5;
  472.     if (tempsec > 30.0)
  473.     (*m)++;
  474.     if (*m == 60) {
  475.     *m = 0;
  476.     (*h)++;
  477.     }
  478. }
  479.  
  480. eq_to_altaz(r, d, t, alt, az)
  481. double r, d, t;
  482. double *alt, *az;
  483. {
  484.     double p = 3.14159265;
  485.     double r1 = p / 180.0;
  486.     double b = lat * r1;
  487.     double l = (360 - lon) * r1;
  488.     double t5, s1, c1, c2, s2, a, h;
  489.  
  490.     if (debug)
  491.     printf("Given R. A. = %lf, DECL. = %lf, gmt = %lf \n", r, d, t);
  492.  
  493.     r = r * 15.0 * r1;
  494.     d = d * r1;
  495.     t = t * 15.0 * r1;
  496.     t5 = t - r + l;
  497.     s1 = sin(b) * sin(d) + cos(b) * cos(d) * cos(t5);
  498.     c1 = 1 - s1 * s1;
  499.     if (c1 > 0) {
  500.     c1 = sqrt(c1);
  501.     h = atan(s1 / c1);
  502.     } else {
  503.     h = (s1 / fabs(s1)) * (p / 2.0);
  504.     }
  505.     c2 = cos(b) * sin(d) - sin(b) * cos(d) * cos(t5);
  506.     s2 = -cos(d) * sin(t5);
  507.     if (c2 == 0) 
  508.     a = (s2/fabs(s2)) * (p/2);
  509.     else {
  510.     a = atan(s2/c2);
  511.     if (c2 < 0)
  512.         a=a+p;
  513.     }
  514.     if (a<0)
  515.         a=a+2*p;
  516.     *alt = h / r1;
  517.     *az = a / r1;
  518.  
  519.     if (debug)
  520.     printf("alt = %lf, az = %lf \n",*alt,*az);
  521. }
  522.  
  523. double
  524. gmst(j, f)
  525. double j,f;
  526. {
  527.     double d, j0, t, t1, t2, s;
  528.  
  529.     d = j - 2451545.0;
  530.     t = d / 36525.0;
  531.     t1 = floor(t);
  532.     j0 = t1 * 36525.0 + 2451545.0;
  533.     t2 = (j - j0 + 0.5)/36525.0;
  534.     s = 24110.54841 + 184.812866 * t1; 
  535.     s += 8640184.812866 * t2;
  536.     s += 0.093104 * t * t;
  537.     s -= 0.0000062 * t * t * t;
  538.     s /= 86400.0;
  539.     s -= floor(s);
  540.     s = 24 * (s + (f - 0.5) * 1.002737909);
  541.     if (s < 0)
  542.     s += 24.0;
  543.     if (s > 24.0)
  544.     s -= 24.0;
  545.  
  546.     if (debug)
  547.     printf("For jd = %lf, f = %lf, gst = %lf \n", j, f, s);
  548.  
  549.     return(s);
  550. }
  551.