/******************************************************************************!
* \file sundial.c
* \author Sebastien Beaugrand
* \sa http://beaugrand.chez.com/
* \copyright CeCILL 2.1 Free Software license
* \note Source 1: Astronomical algorithms - Jean Meeus - 1991
* Source 2: Astronomical algorithms - Jean Meeus - 1998
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "algos.h"
double gX;
double gY;
double gZ;
/******************************************************************************!
* \fn rotx
******************************************************************************/
void rotx(double a)
{
double lY;
a = RAD(a);
lY = gY * cos(a) + gZ * sin(a);
gZ = gZ * cos(a) - gY * sin(a);
gY = lY;
}
/******************************************************************************!
* \fn roty
******************************************************************************/
void roty(double a)
{
double lZ;
a = RAD(a);
lZ = gZ * cos(a) + gX * sin(a);
gX = gX * cos(a) - gZ * sin(a);
gZ = lZ;
}
/******************************************************************************!
* \fn rotz
******************************************************************************/
void rotz(double a)
{
double lX;
a = RAD(a);
lX = gX * cos(a) + gY * sin(a);
gY = gY * cos(a) - gX * sin(a);
gX = lX;
}
/******************************************************************************!
* \fn main
******************************************************************************/
int main(int argc, char* argv[])
{
int year;
int month;
int day;
double lat;
double lon;
double hour;
double straightStylusLength; // cm
double gnomonicDeclination; // west > 0, est < 0
double stylusAngle; // 0: horizontal sundial, 90: vertical
if (argc != 8) {
printf("Usage: %s <YYYY-MM-DD> <latitude> <longitude>"
" <hour> <straight-stylus-length> <declination> <stylus-angle>\n",
argv[0]);
return EXIT_FAILURE;
}
sscanf(argv[1], "%d-%d-%d", &year, &month, &day);
sscanf(argv[2], "%lf", &lat);
sscanf(argv[3], "%lf", &lon);
sscanf(argv[4], "%lf", &hour);
sscanf(argv[5], "%lf", &straightStylusLength);
sscanf(argv[6], "%lf", &gnomonicDeclination);
sscanf(argv[7], "%lf", &stylusAngle);
// Julian Ephemeris Day
double jd = julianDay(year, month, day + hour / 24.0);
// Julian time
double t = julianTime(jd);
// Geometric mean longitude of the Sun
double L = sunGeometricMeanLongitude(t);
// Mean anomaly of the Sun
double M = sunMeanAnomaly(t);
// Sun's equation of center C
double C = sunCenterEquation(t, M);
// True geometric longitude referred to the mean equinox of the date
double theta = sunTrueLongitude(L, C);
// Apparent longitude referred to the true equinox of the date
double lambda = sunApparentLongitude(t, theta);
// Obliquity of the ecliptic
double eps0 = eclipticObliquity(t);
#if \
defined(WITH_EQUATION_OF_TIME) || \
defined(WITH_AZIMUT) || \
defined(WITH_DECLINATION)
// Nutation in longitude and in obliquity
double nutationInLongitude;
double nutationInObliquity;
nutation(t, &nutationInLongitude, &nutationInObliquity);
// Corrected obliquity of the ecliptic
double eps = eps0 + deg(0, 0, nutationInObliquity);
// Apparent declination
double dec = sunApparentDeclination(eps, lambda);
// Apparent right ascension
double ra = sunApparentRightAscension(eps, lambda);
# ifndef NDEBUG
fprintf(stderr, "debug: ra = %.7f, dec = %.7f\n", ra, dec);
# endif
#endif
#if defined(WITH_EQUATION_OF_TIME)
// 1: (56) p372
// 2: (58) p402
double H = (hour - 12) * 15 + lon +
equationOfTime(L, ra, nutationInLongitude, eps);
double P =
SIN(lat) * COS(stylusAngle) -
COS(lat) * SIN(stylusAngle) * COS(gnomonicDeclination);
double Q = SIN(gnomonicDeclination) * SIN(stylusAngle) * SIN(H) +
(COS(lat) * COS(stylusAngle) +
SIN(lat) * SIN(stylusAngle) * COS(gnomonicDeclination)) * COS(H) +
P * TAN(dec);
if (Q > 0.1) {
double nx =
COS(gnomonicDeclination) * SIN(H) -
SIN(gnomonicDeclination) *
(SIN(lat) * COS(H) - COS(lat) * TAN(dec));
double ny = COS(stylusAngle) * SIN(gnomonicDeclination) * SIN(H) -
(COS(lat) * SIN(stylusAngle) - SIN(lat) *
COS(stylusAngle) * COS(gnomonicDeclination)) * COS(H) -
(SIN(lat) * SIN(stylusAngle) + COS(lat) *
COS(stylusAngle) * COS(gnomonicDeclination)) * TAN(dec);
gX = -nx;
gY = -ny;
gZ = Q;
}
#elif defined(WITH_AZIMUT)
// Apparent Sideral Time
double theta0 = apparentSideralTime(reduceAngle(sideralTime(jd, t)),
nutationInLongitude, eps);
// Altitude
double h = altitude(ra, dec, theta0, lat, -lon);
// Azimut
double a = azimut(ra, dec, theta0, lat, -lon);
# ifndef NDEBUG
fprintf(stderr,
"debug: theta0 = %.7f, azi = %.7f, ele = %.7f\n",
theta0,
a,
h);
# endif
gX = 0;
gY = 0;
gZ = 1;
rotx(h);
roty(a - gnomonicDeclination);
rotx(stylusAngle - 90);
#else
# if defined(WITH_DECLINATION)
gX = COS(ra) * COS(dec);
gY = SIN(ra) * COS(dec);
gZ = SIN(dec);
# else
double eps = eclipticObliquityCorrected(t, eps0);
gX = COS(lambda);
gY = COS(eps) * SIN(lambda);
gZ = SIN(eps) * SIN(lambda);
# endif
rotz(L);
rotz((hour - 24) * 15 + lon - 180);
roty(90 - lat);
rotz(-gnomonicDeclination);
roty(stylusAngle);
rotz(90);
#endif
if (gZ > 0) {
double x = -straightStylusLength * gX / gZ;
double y = -straightStylusLength * gY / gZ;
double sousStylaire = straightStylusLength * TAN(lat);
if (((month >= 3 && month < 6) || (month == 6 && day < 21)) &&
(hour - truncf(hour) < 0.1)) {
if (x > -29.5 && x < 29.5 &&
y > -28.5 + sousStylaire && y < sousStylaire) {
printf("%.7f %.7f\n", x, y);
}
} else {
if (x > -29.9 && x < 29.9 &&
y > -29.5 + sousStylaire && y < sousStylaire) {
printf("%.7f %.7f\n", x, y);
}
}
}
return EXIT_SUCCESS;
}