wmweather+/sunzenith.c

   1 #include "config.h"
   2
   3 /*  Copyright (C) 2002  Brad Jorsch <anomie@users.sourceforge.net>
   4
   5     This program is free software; you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation; either version 2 of the License, or
   8     (at your option) any later version.
   9
  10     This program is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License
  16     along with this program; if not, write to the Free Software
  17     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  18 */
  19
  20 /* Algorithms from http://www.srrb.noaa.gov/highlights/sunrise/azel.html */
  21
  22 #include <stdio.h>
  23 #include <math.h>
  24 #include <time.h>
  25 #include <stdlib.h>
  26 #include <limits.h>
  27
  28 #include "convert.h"
  29
  30 /*  Purpose: calculate the Geometric Mean Longitude of the Sun (degrees) */
  31 double calcGeomMeanLongSun(double t) {
  32     double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
  33     while(L0 > 360.0) {
  34         L0 -= 360.0;
  35     }
  36     while(L0 < 0.0) {
  37         L0 += 360.0;
  38     }
  39     return L0;
  40 }
  41
  42
  43 /*  Purpose: calculate the Geometric Mean Anomaly of the Sun (degrees) */
  44 double calcGeomMeanAnomalySun(double t) {
  45     return 357.52911 + t * (35999.05029 - 0.0001537 * t);
  46 }
  47
  48
  49 /*  Purpose: calculate the eccentricity of earth's orbit */
  50 double calcEccentricityEarthOrbit(double t) {
  51     return 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  52 }
  53
  54
  55 /*  Purpose: calculate the equation of center for the sun (degrees) */
  56 double calcSunEqOfCenter(double t) {
  57     double m = deg2rad(calcGeomMeanAnomalySun(t));
  58
  59     return sin(m) * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin(m+m) * (0.019993 - 0.000101 * t) + sin(m+m+m) * 0.000289;
  60 }
  61
  62
  63 /*  Purpose: calculate the true longitude of the sun (degrees) */
  64 double calcSunTrueLong(double t) {
  65     return calcGeomMeanLongSun(t) + calcSunEqOfCenter(t);
  66 }
  67
  68
  69 /*  Purpose: calculate the apparent longitude of the sun (degrees) */
  70 double calcSunApparentLong(double t) {
  71     return calcSunTrueLong(t) - 0.00569 - 0.00478 * sin(deg2rad(125.04-1934.136*t));
  72 }
  73
  74
  75 /*  Purpose: calculate the mean obliquity of the ecliptic (degrees) */
  76 double calcMeanObliquityOfEcliptic(double t) {
  77     return 23.0 + (26.0 + ((21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813))))/60.0))/60.0;
  78 }
  79
  80
  81 /*  Purpose: calculate the corrected obliquity of the ecliptic (degrees) */
  82 double calcObliquityCorrection(double t) {
  83     return calcMeanObliquityOfEcliptic(t) + 0.00256*cos(deg2rad(125.04-1934.136*t));
  84 }
  85
  86
  87 /*  Purpose: calculate the declination of the sun (degrees) */
  88 double calcSunDeclination(double t) {
  89     return rad2deg(asin(sin(deg2rad(calcObliquityCorrection(t))) *
  90                          sin(deg2rad(calcSunApparentLong(t)))));
  91 }
  92
  93
  94 /*  Purpose: calculate the difference between true solar time and mean
  95  *              solar time (minutes)
  96  */
  97 double calcEquationOfTime(double t) {
  98     double l0 = deg2rad(calcGeomMeanLongSun(t));
  99     double e = calcEccentricityEarthOrbit(t);
 100     double m = deg2rad(calcGeomMeanAnomalySun(t));
 101     double y = tan(deg2rad(calcObliquityCorrection(t))/2.0);
 102     double sinm = sin(m);
 103
 104     y *= y;
 105
 106     return rad2deg(y*sin(2.0*l0) - 2.0*e*sinm + 4.0*e*y*sinm*cos(2.0*l0)
 107                     - 0.5*y*y*sin(4.0*l0) - 1.25*e*e*sin(2.0*m))*4.0;
 108 }
 109
 110
 111 double calcSolarZenith(double latitude, double longitude, int year, int month, int day, int timeUTC){
 112     double T, trueSolarTime, hourAngle, solarDec, csz, zenith, exoatmElevation, te, refractionCorrection;
 113
 114     T=jd2jcent(mdy2jd(year, month, day) + timeUTC/1440.0);
 115     trueSolarTime = timeUTC + calcEquationOfTime(T) - 4.0 * longitude;
 116     hourAngle = trueSolarTime / 4.0 - 180.0;
 117     solarDec = calcSunDeclination(T);
 118     csz = sin(deg2rad(latitude)) * sin(deg2rad(solarDec)) +
 119         cos(deg2rad(latitude)) * cos(deg2rad(solarDec)) *
 120         cos(deg2rad(hourAngle));
 121     zenith=rad2deg(acos(csz));
 122     exoatmElevation = 90.0 - zenith;
 123     if (exoatmElevation > 85.0) {
 124         refractionCorrection = 0.0;
 125     } else {
 126         te = tan(deg2rad(exoatmElevation));
 127         if (exoatmElevation > 5.0) {
 128             refractionCorrection = 58.1/te - 0.07/(te*te*te) +
 129                 0.000086/(te*te*te*te*te);
 130         } else if (exoatmElevation > -0.575) {
 131             refractionCorrection = 1735.0 + exoatmElevation*(-518.2 + exoatmElevation*(103.4 + exoatmElevation*(-12.79 + exoatmElevation*0.711)));
 132         } else {
 133             refractionCorrection = -20.774 / te;
 134         }
 135         refractionCorrection = refractionCorrection / 3600.0;
 136     }
 137     return zenith - refractionCorrection;
 138 }