/********************************
* Projekt: Funkuhr DCF77 *
* Mondalter-Berechnung zur *
* korrekten Mondphasen-Anzeige *
* sowie Berechnung von *
* Auf- und Untergangszeiten *
********************************/
#include "Funkuhr.h"
#include <math.h>
#define CS(x) cos((x)*PI180)
#define SN(x) sin((x)*PI180)
#define FRAK(x) fmod(x,1)
PointF geoloc; // Geoposition in Grad, für AufUnterBerechnung()
int geotz; // wie Windows in Minuten, östlich = negativ (bekloppt, aber Redmont ist nicht Moskau)
TIME_ZONE_INFORMATION geotzi; // mit entsprechend umgerechneten Umschalt-Tagen
static double cphi;
static double sphi;
/* JavaScript-Routinen geklaut von www.computus.de und leider nicht ganz verstanden */
double Mondalter(double JDE){
double t, d, m, m1, i;
t = JDE / 36525.0;
d = 297.8502042 + 445267.11151686 * t - .00163 * t * t + t * t * t / 545868.0 - t * t * t * t / 113065000.0;
m = 357.5291092 + 35999.0502909 * t - .0001536 * t * t + t * t * t / 24490000.0;
m1 = 134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699.0 - t * t * t * t / 14712000.0;
i = d + 6.289 * SN(m1) - 2.1 * SN(m) + 1.274 * SN(2 * d - m1) + .658 * SN(2 * d) + .241 * SN(2 * m1) + .110 * SN(d);
return fmod(i,360); // Beleuchtungsgrad dazu wäre: (1-CS(i))/2;
}
typedef struct{
double dec,ra;
}DECRA;
// Position (Sternenposition?) des Mondes zur Zeit <t>
static void Mondberechnung(double t, DECRA*decra){
double p2 = 6.283185307;
double arc = 206264.8062;
double coseps = .91748;
double sineps = .39778;
double lo = FRAK(.606433 + 1336.855225 * t);
double l = p2 * FRAK(.374897 + 1325.55241 * t);
double ls = p2 * FRAK(.993133 + 99.997361 * t);
double d = p2 * FRAK(.827361 + 1236.853086 * t);
double f = p2 * FRAK(.259086 + 1342.227825 * t);
double dl = 22640 * sin(l) - 4586 * sin(l - 2 * d) + 2370 * sin(2 * d) + 769 * sin(2 * l) - 668 * sin(ls) - 412 * sin(2 * f) - 212 * sin(2 * l - 2 * d) - 206 * sin(l + ls - 2 * d) + 192 * sin(l + 2 * d) - 165 * sin(ls - 2 * d) - 125 * sin(d) - 110 * sin(l + ls) + 148 * sin(l - ls) - 55 * sin(2 * f - 2 * d);
double s = f + (dl + 412 * sin(2 * f) + 541 * sin(ls)) / arc;
double h = f - 2 * d;
double n = -526 * sin(h) + 44 * sin(l + h) - 31 * sin(-l + h) - 23 * sin(ls + h) + 11 * sin(-ls + h) - 25 * sin(-2 * l + f) + 21 * sin(-l + f);
double lmoon = p2 * FRAK(lo + dl / 1296000);
double bmoon = (18520 * sin(s) + n) / arc;
double cb = cos(bmoon);
double x = cb * cos(lmoon);
double v = cb * sin(lmoon);
double w = sin(bmoon);
double y = coseps * v - sineps * w;
double z = sineps * v + coseps * w;
double rho = sqrt(1 - z * z);
decra->dec = (360 / p2) * atan(z / rho);
decra->ra = (48 / p2) * atan(y / (x + rho));
if (decra->ra < 0) decra->ra += 24;
}
static void Sonnenberechnung(double t, DECRA*decra){
double p2 = 6.283185307; // s.o.
double coseps = .91748; // s.o.
double sineps = .39778; // s.o.
double m = p2 * FRAK(.993133 + 99.997361 * t); // s.o.
double dl = 6893 * sin(m) + 72 * sin(2 * m);
double l = p2 * FRAK(.7859453 + m / p2 + (6191.2 * t + dl) / 1296000);
double sl = sin(l);
double x = cos(l);
double y = coseps * sl;
double z = sineps * sl;
double rho = sqrt(1 - z * z);
decra->dec = (360 / p2) * atan(z / rho);
decra->ra = (48 / p2) * atan(y / (x + rho));
if (decra->ra < 0) decra->ra += 24;
}
static double LMST(double MJD) {
double MJDO = lrint(MJD-0.5);
double ut = (MJD - MJDO) * 24;
double t = (MJDO - 51544.5) / 36525; // 1 ~ 100 Jahre
double gmst = 6.697374558 + 1.0027379093 * ut + (8640184.812866 + (.093104 - .0000062 * t) * t) * t / 3600;
return 24 * FRAK((gmst + geoloc.X / 15) / 24);
}
static double SINALT(bool bSonne, double MJD, double Hour) {
double MJDO = MJD + Hour / 24; // Gebrochenes Modifiziertes Julianisches Datum
double t = (MJDO - 51544.5) / 36525; // 1 ~ 100 Jahre
DECRA decra;
double tau;
(bSonne?Sonnenberechnung:Mondberechnung)(t,&decra);
tau = 15 * (LMST(MJDO) - decra.ra);
return sphi * SN(decra.dec) + cphi * CS(decra.dec) * CS(tau);
}
typedef struct{
double ye,zero1,zero2;
int nz;
}QUADRET;
static void QUAD(double yminus, double yo, double yplus, QUADRET*quadret) {
double a = .5 * (yminus + yplus) - yo;
double b = .5 * (yplus - yminus);
double c = yo;
double xe = -b / (2 * a);
double dis = b * b - 4 * a * c;
quadret->ye = (a * xe + b) * xe + c;
quadret->nz = 0;
quadret->zero1 = 0;
quadret->zero2 = 0;
if (dis >= 0){
double dx = .5 * sqrt(dis) / fabs(a);
quadret->zero1 = xe - dx;
quadret->zero2 = xe + dx;
if (fabs(quadret->zero1) <= 1) quadret->nz++;
if (fabs(quadret->zero2) <= 1) quadret->nz++;
if (quadret->zero1 < -1) quadret->zero1 = quadret->zero2;
}
}
typedef struct{
double rise,set;
}RISESET;
// liefert <true> wenn um 00:00 Uhr sichtbar
// Routine schwächelt wenn (Mond)Aufgangszeit exakt 00:00 Uhr ist
static bool AufUnterBerechnung(bool bSonne, double MJD, double sinho, RISESET*RiseSet) {
double Hour = 1;
double yo, yplus;
QUADRET QuatRet;
double yminus = SINALT(bSonne, MJD, 0) - sinho;
bool above = yminus > 0;
RiseSet->rise=0;
RiseSet->set=0;
do {
yo = SINALT(bSonne, MJD, Hour) - sinho;
yplus = SINALT(bSonne, MJD, Hour + 1) - sinho;
QUAD(yminus, yo, yplus, &QuatRet);
switch (QuatRet.nz) { // Anzahl detektierter Nulldurchgänge
case 1: {
if (yminus < 0) RiseSet->rise = Hour + QuatRet.zero1;
else RiseSet->set = Hour + QuatRet.zero1;
}break;
case 2: {
if (QuatRet.ye < 0) {
RiseSet->rise = Hour + QuatRet.zero2;
RiseSet->set = Hour + QuatRet.zero1;
}else{
RiseSet->rise = Hour + QuatRet.zero1;
RiseSet->set = Hour + QuatRet.zero2;
}
}break;
}
yminus = yplus;
Hour += 2;
}while (Hour < 25 && (!RiseSet->rise || !RiseSet->set));
return above;
}
static int HH_MM(double hour, TCHAR*s, int len) {
if (hour) {
int m=lrint(hour*60);
return wnsprintf(s,len,T("%02d%s%02d"),m/60,sTime,m%60);
}else{
return wnsprintf(s,len,T("--%s--"),sTime);
}
}
bool AufUntergang(bool bSonne, double MJD, TCHAR*Auf, TCHAR*Unter) {
double sinho; // Sinus des Horizonts
RISESET RiseSet;
bool ret;
static const float horizont[4]={
8.0f/60, //Mond: 8 Bogenminuten über geometrischem Horizont
-50.0f/60, //Sonne: 50 Bogenminuten unter geom. Horizont
-6.0f, //Dämmerung (bürgerlich): 6° unterm Horizont (hier nicht verwendet)
-12.0f}; //Dämmerung (nautisch): 12°
// Modifiziertes Julianisches Datum == Julianisches Datum - 2400000.5
// Mittwoch, 17. Nov 00:00 UT1 1858
MJD += geotz / 1440.0f;
sphi = SN(geoloc.Y); // nicht thread-sicher
cphi = CS(geoloc.Y);
sinho = SN(horizont[bSonne&3]);
ret=AufUnterBerechnung(bSonne, MJD, sinho, &RiseSet);
HH_MM(RiseSet.rise,Auf,6);
HH_MM(RiseSet.set,Unter,6);
return ret;
}
// EINGABE:
// date.wMonth = Umschalt-Monat
// date.wDayOfWeek = Wochentag des Umschaltens
// date.wDay = Nummer des Wochentags im Monat:
// 1..4 = 1. bis 4. Wochentag
// 5..8 = letzter, vorletzter usw. Wochentag
// date.wYear = Jahr, das abgefragt werden soll (einsetzen!)
// AUSGABE
// date.wDay aktualisiert, also 1..31
BOOL GetTimezoneChangeDate(MYSYSTEMTIME *st) {
FILETIME ft;
BYTE nWeek=st->bDay-1; // Nummer der Woche im Monat -1
BYTE day;
if (!st->bMonth || !st->st.wYear || (unsigned)nWeek>=8) return FALSE;
day=st->bDayOfWeek; // retten (0..6)
if (nWeek>=4) { // letzter, vorletzter usw. Wochentag im Monat
BYTE last;
nWeek=3-nWeek; // negativ machen (typisch: -1)
// Letzten Tag des Monats berechnen (geht nur so)
for (st->bDay=31; !SystemTimeToFileTime(&st->st,&ft); st->bDay--) {
if (st->bDay==28) return FALSE; // Problem! Notbremse
}
FileTimeToSystemTime(&ft,&st->st);
last=st->bDayOfWeek;
st->bDayOfWeek=day;
day-=last;
if (day<=0) nWeek++;
}else{ // Wochentag in 1.-4. Woche im Monat
BYTE first;
// Wochentag des ersten Tages im Monat berechnen
st->bDay=1;
if (!SystemTimeToFileTime(&st->st,&ft)) return FALSE;
FileTimeToSystemTime(&ft,&st->st);
first=st->bDayOfWeek;
st->bDayOfWeek=day;
day-=first;
if (day<0) nWeek++;
}
st->bDay+=day+nWeek*7;
return TRUE;
}
// Initialisiert <geotz> und <geotzi> u.a. mit echten Umschalt-Daten für das aktuelle Jahr
void InitTZ() {
geotz = 0;
switch (GetTimeZoneInformation(&geotzi)) {
case 1: geotz = geotzi.StandardBias; break;
case 2: geotz = geotzi.DaylightBias; break;
}
geotz += geotzi.Bias;
if (geotzi.StandardDate.wMonth) {
SYSTEMTIME st;
GetLocalTime(&st);
geotzi.StandardDate.wYear=st.wYear;
GetTimezoneChangeDate((MYSYSTEMTIME*)&geotzi.StandardDate);
geotzi.DaylightDate.wYear=st.wYear;
GetTimezoneChangeDate((MYSYSTEMTIME*)&geotzi.DaylightDate);
}
}
/*
int mainCRTStartup() {
FILETIME ft;
GetSystemTimeAsFileTime(&ft);
// zeitzone=int((*(__int64*)&ft-*(__int64*)&ftl)>>8)/2343750;
double MJD=*(__int64*)(&ft)/(24*60*60*10000000.0)-94187;
char sx[8],sy[8],sz[8];
d2s(sx,geolocation.x);
d2s(sy,geolocation.y);
d2s(sz,zeitzone/60.0);
myprintf("Koordinaten: %s° / %s°, Zeitzone: %s\n",sx,sy,sz);
myprintf("Mondalter jetzt:%4d°\n", (int)RoundInt(Mondalter(MJD-51544.5)));
MJD=myfloor(MJD);
for (int i=0; i<45; i++) {
double NachtMJD = MJD+i;
char sa[6],su[6],ma[6],mu[6],db[6],de[6],date[14];
SYSTEMTIME st;
*(__int64*)&ft=RoundInt((NachtMJD+94187)*(24*60*60*10000000.0));
FileTimeToSystemTime(&ft,&st);
if (tzi.StandardDate.wMonth) {
if (st.wMonth==tzi.StandardDate.wMonth
&& st.wDay==tzi.StandardDate.wDay) {
zeitzone=tzi.Bias+tzi.StandardBias;
d2s(sz,zeitzone/60.0);
myprintf("Zeitzone: %s\n",sz);
}
if (st.wMonth==tzi.DaylightDate.wMonth
&& st.wDay==tzi.DaylightDate.wDay) {
zeitzone=tzi.Bias+tzi.DaylightBias;
d2s(sz,zeitzone/60.0);
myprintf("Zeitzone: %s\n",sz);
}
}
AufUntergang(2,NachtMJD,sa,su);
AufUntergang(1,NachtMJD,ma,mu);
AufUntergang(3,NachtMJD,db,de);
GetDateFormat(LOCALE_USER_DEFAULT,0,&st,"d. MMMM",date,14);
myprintf("%13s:%4d°, SA=%s, SU=%s, MA=%s, MU=%s, D=%s-%s\n",
date, (int)RoundInt(Mondalter(NachtMJD-51544.5)),sa,su,ma,mu,db,de);
}
return 0;
}
*/
| Vorgefundene Kodierung: ANSI (CP1252) | 4
|
|
|