Close

The Mathematics

A project log for Daylight Geolocator Remix

This is a remix of https://hackaday.io/project/28550-light-level-geolocator

agpcooperagp.cooper • 02/04/2018 at 09:04•0 Comments

Julian Day (JD) number

Most astronomy based formulas use the Julian day number as an input.

I use the a modified 2000 Julian Day number, the number of days since mid-night of the 1st of January 2000:

double J2000(long Year,long Month,long Day,double UT) {
  double JD;
  
  // Julian (2000) Day Number
  JD=(367*Year-(7*(Year+(Month+9)/12))/4+275*Month/9+Day);
  JD=JD-730531.5;
  JD=JD+UT/24.0;
  
  return JD;
}

The formular is correct even for dated before 1/1/2000.

One of the inputs to the function is UT or Universal Time. UT is the time in Greenwich UK. I will write more about this later.

Note that Year, Month and Day are type long. Therefore integer division rules apply.

Reverse of the Julian Day number

I found a source for calculating the reverse of the 2000 Julian Day Number:

Source: http://mathforum.org/library/drmath/view/51907.htm

  unsigned long JD,p,q,r,s,t,u,v;

  p=JD+2451545+68569;
  q=4*p/146097;
  r=p-(146097*q+3)/4;
  s=4000*(r+1)/1461001;
  t=r-1461*s/4+31;
  u=80*t/2447;
  v=u/11;
  Year=100*(q-49)+s+v;
  Month=u+2-12*v;
  Day=t-2447*u/80;

It may be useful when checking/calculating valid dates. I checked the code and it good from 1/1/2000 to 31/12/2099 but then the Julian Day number calculation is only good to 31/12/2099.

Greenwich Mean Sideral Time (GMST)

This is the time of the stars! You would be aware that there are 365.25 days in a solar year. But the stars rotate around the earth 366.25 times a year. The extra rotation is because the earth also rotates around the sun as well as rotating around it axis. Here is the GMST function:

double GMST(double JD,double UT) {
  double gmst;
  
  // Greenwich Mean Sidereal Time (~36625/36525)
  gmst=18.697374558+24.06570982441908*(JD+UT/24.0);
  gmst=fmod(fmod(gmst,24.0)+24.0,24.0);
  
  return gmst;
}

 We need GMST even for the sun as the sun is cosidered just another start. But the Right Ascension and Declination moves a lot more than ordinary stars.

Sun Right Ascension (RA) and Declination (Dec)

Now we can calculate the position of the sun amongst the stars. I use a series of formulas from my copy of the "The Astronomical Amanac (1986)" but I have seen the same formulas on the Internet:

void SunData(double JD,double *RA,double *Dec,double *EoT) {
  double L,g,a,b,c;
  
  /* Source: Page C24 from "The Astronomical Almanac" (1986) */
  // Mean longitude of Sun, corrected for aberration
  L=280.460+0.9856474*JD;
  L=fmod(fmod(L,360.0)+360.0,360.0);
  // Mean anomaly
  g=357.528+0.9856003*JD;
  g=fmod(fmod(g,360.0)+360.0,360.0);
  // Ecliptic longitude
  a=L+1.915*sin(g*M_PI/180)+0.020*sin(2*g*M_PI/180);
  a=fmod(fmod(a,360.0)+360.0,360.0);
  // Ecliptic latitude
  b=0;
  // Obliquity of ecliptic
  c=23.439-0.0000004*JD;
  c=fmod(fmod(c,360.0)+360.0,360.0);
  // Right ascension
  (*RA)=atan2(cos(c*M_PI/180)*sin(a*M_PI/180),cos(a*M_PI/180))*12/M_PI;
  if ((*RA)<0) (*RA)=(*RA)+24.0;
  // Equation of Time
  (*EoT)=(L/15.0-(*RA));
  // Declination
  (*Dec)=asin(sin(c*M_PI/180)*sin(a*M_PI/180))*180/M_PI;
}

The procedure accepts modified Julian Day number (JD) and returns Right Ascension (RA), Declination (Dec) and the Equarion of Time (EoT). The units for RA and EoT are hours but the units for Declination is degrees.

Checking the formulas

It is really easy to make a minor mistake with these formulas so it is important to check the results against a worked example. The Internet has many calculators that you can check your results against.

Sunrise and sunset

Sunrise and sunset times are published on the Internet for any place that you live.

On the 20th of January 2018 in Perth Western Australia the Sunrise was 5:30 and Sunset was 19:24. Local noon was therefore 12:27. The local time is 8 hours ahead of UT.

Calculating  longitude

The steps are:

  1.   SunRise=5.50-8;
  2.   SunSet=19.40-8;
  3.   JD=J2000(Year,Month,Day,(SunRise+SunSet)/2.0);
  4.   SunData(JD,&RA,&Dec,&EoT);
  5.   Lon=15.0*(RA-GMST(JD,0.0));

The values were:

Which are correct.

Calculating latitude

Originally I used an iterative method but Jaromir Sukuba work out an exact answer:

double latitude(double SunRise,double SunSet,double Alt,double Dec) {
  // The mathematics:  
  //   Cos(15*DL/2)=Sin(Alt)/Cos(Dec)/Cos(Lat)-Tan(Dec)*Tan(Lat)
  //   Use this model:
  //     C=B*Cos(Lat)+A*Sin(Lat)
  //     D=Sqrt(A*A+B*B);
  //     Sin(U)=B/D;
  //     Sin(V)=C/D;
  //     U=ASIN(B/D);
  //     V=ASIN(C/D);
  //     Lat=V-U;
  //  Rearrange:
  //     Sin(Alt)/Cos(Dec)=Cos(15*DL/2)*Cos(Lat)+Tan(Dec)*Sin(Lat))
  //  Test data:
  //     SunRise       5.50 (need not be local time)
  //     SunSet       19.40 (need not be local time)
  //     DayLight     13.90
  //     Altitude     -0.833 (definition for Sun rise/set)
  //     Declination -20.066095
  //     Latitude    -31.9615789037
  double DL,A,B,C,D,U,V,Lat;

  if (fabs(Dec-Lat)>90.0) return -90.0; // Check if sun rises 
  if (fabs(Dec+Lat)>90.0) return +90.0; // Check if sun sets
  DL=SunSet-SunRise;                    //  13.90
  A=tan(Dec*PI/180);                    //  -0.3652771781
  B=cos(DL*PI/24);                      //  -0.246153293  
  C=sin(Alt*PI/180)/cos(Dec*PI/180);    //  -0.015477611  
  D=sqrt(A*A+B*B);                      //   0.4404757207
  if (fabs(C)<=D) {                     // Check if out of range
    U=asin(B/D)*180/PI;                 // -33.9752753119  
    V=asin(C/D)*180/PI;                 //  -2.0136964081
    Lat=V-U;                            // -31.9615789037
    if (Dec<0) Lat=-Lat;                // -31.9615789037  
  } else {
    Lat=-90;
  }
  return Lat;
}

If we apply our numbers to this function then the answer is:

  6. Lat=latitude(SunRise,SunSet,Alt,Dec);

Which is correct.

The only variable not discussed here is Alititude (Alt). Altitude is the angle from the horizon. For sunrise and sunset the offical altitude (nainlt due to reflaction) is -0.833 degrees (i.e. below the horizon), not 0 degrees.

Some code to try

Here is some arduino code:

double J2000(long Year,long Month,long Day,double UT) {
  double JD;
  
  // Julian (2000) Day Number
  JD=(367*Year-(7*(Year+(Month+9)/12))/4+275*Month/9+Day);
  JD=JD-730531.5;
  JD=JD+UT/24.0;
  
  return JD;
}

double GMST(double JD,double UT) {
  double gmst;
  
  // Greenwich Mean Sidereal Time (~36625/36525)
  gmst=18.697374558+24.06570982441908*(JD+UT/24.0);
  gmst=fmod(fmod(gmst,24.0)+24.0,24.0);
  
  return gmst;
}

void SunData(double JD,double *RA,double *Dec,double *EoT) {
  double L,g,a,b,c;
  
  /* Source: Page C24 from "The Astronomical Almanac" (1986) */
  // Mean longitude of Sun, corrected for aberration
  L=280.460+0.9856474*JD;
  L=fmod(fmod(L,360.0)+360.0,360.0);
  // Mean anomaly
  g=357.528+0.9856003*JD;
  g=fmod(fmod(g,360.0)+360.0,360.0);
  // Ecliptic longitude
  a=L+1.915*sin(g*M_PI/180)+0.020*sin(2*g*M_PI/180);
  a=fmod(fmod(a,360.0)+360.0,360.0);
  // Ecliptic latitude
  b=0;
  // Obliquity of ecliptic
  c=23.439-0.0000004*JD;
  c=fmod(fmod(c,360.0)+360.0,360.0);
  // Right ascension
  (*RA)=atan2(cos(c*M_PI/180)*sin(a*M_PI/180),cos(a*M_PI/180))*12/M_PI;
  if ((*RA)<0) (*RA)=(*RA)+24.0;
  // Equation of Time
  (*EoT)=(L/15.0-(*RA));
  // Declination
  (*Dec)=asin(sin(c*M_PI/180)*sin(a*M_PI/180))*180/M_PI;
}

double latitude(double SunRise,double SunSet,double Alt,double Dec) {
  // The mathematics:  
  //   Cos(15*DL/2)=Sin(Alt)/Cos(Dec)/Cos(Lat)-Tan(Dec)*Tan(Lat)
  //   Use this model:
  //     C=B*Cos(Lat)+A*Sin(Lat)
  //     D=Sqrt(A*A+B*B);
  //     Sin(U)=B/D;
  //     Sin(V)=C/D;
  //     U=ASIN(B/D);
  //     V=ASIN(C/D);
  //     Lat=V-U;
  //  Rearrange:
  //     Sin(Alt)/Cos(Dec)=Cos(15*DL/2)*Cos(Lat)+Tan(Dec)*Sin(Lat))
  //  Test data:
  //     SunRise       5.50 (need not be local time)
  //     SunSet       19.40 (need not be local time)
  //     DayLight     13.90
  //     Altitude     -0.833 (definition for Sun rise/set)
  //     Declination -20.066095
  //     Latitude    -31.9615789037
  double DL,A,B,C,D,U,V,Lat;

  DL=SunSet-SunRise;                    //  13.90
  A=tan(Dec*PI/180);                    //  -0.3652771781
  B=cos(DL*PI/24);                      //  -0.246153293  
  C=sin(Alt*PI/180)/cos(Dec*PI/180);    //  -0.015477611  
  D=sqrt(A*A+B*B);                      //   0.4404757207
  if (fabs(C)<=D) {                     // Check if out of range
    U=asin(B/D)*180/PI;                 // -33.9752753119  
    V=asin(C/D)*180/PI;                 //  -2.0136964081
    Lat=V-U;                            // -31.9615789037
    if (Dec<0) Lat=-Lat;                // -31.9615789037  
  } else {
    Lat=0.0;
  }
  return Lat;
}

void setup()  {
  Serial.begin(9600);
  while (!Serial);

  // Day Light to Lat/Long
  int Year,Month,Day;
  double LTC,UT,JD,SunRise,SunSet,RA,Dec,EoT,Lat,Lon,Alt;

  LTC=8.0;  // Local Time Correction
  Year=2018;
  Month=1;
  Day=20; 
  Alt=-0.833;                // Sunrise and Sunset definition
  SunRise=5.0+30.0/60.0-LTC; // Convert Local Time to Universal Time
  SunSet=19.0+24.0/60.0-LTC; // Convert Local Time to Universal Time

  // Problem parameters:
  Serial.println("DayLight to Latitude:");
  Serial.print("Year        ");Serial.println(Year);
  Serial.print("Month       ");Serial.println(Month);
  Serial.print("Day         ");Serial.println(Day);
  Serial.print("SunRise(UT) ");Serial.println(SunRise,4);
  Serial.print("SunSet (UT) ");Serial.println(SunSet,4); 
   
  // Solve the problem:
  JD=J2000(Year,Month,Day,(SunRise+SunSet)/2.0);
  SunData(JD,&RA,&Dec,&EoT);
  Lat=latitude(SunRise,SunSet,Alt,Dec);
  Lon=15.0*(RA-GMST(JD,0.0));

  // Debug:
  Serial.print("JD          ");Serial.println(JD,8);
  Serial.print("RA          ");Serial.println(RA,8);
  Serial.print("Dec         ");Serial.println(Dec,8);
  Serial.print("EoT         ");Serial.println(EoT,8);

  // Solution:
  Serial.print("Latitude    ");Serial.println(Lat,4);
  Serial.print("Longitude   ");Serial.println(Lon,4);
  Serial.println(); 
}

void loop(){    

}

And here is the run:

DayLight to Latitude:
Year        2018
Month       1
Day         20
SunRise(UT) -2.5000
SunSet (UT) 11.4000
JD          6593.68554687
RA          20.14925956
Dec         -20.13604354
EoT         -0.18200683
Latitude    -31.8655
Longitude   115.9108

Conclusion

For many places in the world if you can measure sunrse and sunset accuarately you can calculate your latitude and longitude accurately. There are places/dates where the sun does not set or rise, in these places/dates the mathematics does not work.

Magic
Next log

The Box

02/04/2018 at 11:37 • 0 comments

Discussions

Yes, delete it Cancel

Going up?

About Us Contact Hackaday.io Give Feedback Terms of Use Privacy Policy Hackaday API

© 2024 Hackaday