Calculate the Time When the Sun Is X Degrees Below/Above the Horizon

Calculate whether it is close to dawn/dusk based on sunrise/sunset?

there are some simplified equations for that but the real thing is this:

• solving Kepler's equation
• btw your question is very similar to this: time when the sun is X degrees below/above Horizon

if you want easy solution then obtain sunrise calendar for your location

it should look something like this:

Just pick the closest latitude or generate the right diagram for your location. The picture was taken from here so there should be also the generation applet somewhere. Can't find any image of the astronomy sun calendar on the WEB so here is scan from really old book of mine:

Booth are the same thing but the scan is just for single long/lat geo-location.

• x axis is the date

• y axis is the time sun reach (sun rise/set or dusk/dawn,...)

• Noc means night and Den means day

• all the marks are unimportant for you

• dash dash line is sun rise/set (leading visual edge at mathematical horizont)

• dash dot line is the civil twilight (visual center of sun is 6 deg below horizont)

• filled line is astronomical twilight (visual center of sun is 18 deg below horizont)

  Diagram is for latitudes 48,49,50 deg, nautical twilight (visual center of sun is 12 deg below horizon) is not on the image. When you approximate this curve by polynomial or piecewise polyline then you can easily compute what you need per any date unless you want to use this for thousands of years then you can take this as periodic function so you need just one year table and all the other years are the same.

[Notes]

visual

• means geometric position + atmospheric refraction + aberation + precession/nutation + speed of light
• near horizon is the difference little more then 0.5 degree !!! (one Sun disc)
• and the time difference is around 8 minutes

Calculations of the path of the Sun

You pass the longitude to localSolarTime() as degrees, and then you divide that by 60, with a comment claiming this is in order to convert to minutes of arc. This is wrong; your later calculations require degrees, and even if you needed minutes of arc, you'd multiply by 60, not divide.

This mistaken division results in a longitude of -1.3°, and when you find the angle between your local time meridian and your position, you get a large angle (about 75°). It should be a small angle, generally ±7.5°. The large angle results in a large time correction, and throws everything off.

---

Update: In the updated version of the azimuth() method, the quadrant selection should be based on the hour angle of the sun, or, equivalently, on local solar time, rather than standard wall clock time. And, the hour angle used in all calculations should not be rounded. Rather than testing four different quadrants, the method could look like this:

public static double azimuth(double lat, double declination, double zenith, double hourAngle)
{
  double elevation = Math.toRadians(90 - zenith);
  lat = Math.toRadians(lat);
  declination = Math.toRadians(declination);
  hourAngle = Math.toRadians(hourAngle);
  double azimuthRadian = acos(((sin(declination) * cos(lat)) - (cos(hourAngle) * cos(declination) * sin(lat))) / cos(elevation));
  double azimuthDegree = Math.toDegrees(azimuthRadian);
  if (hourAngle > 0)
    azimuthDegree = 360 - azimuthDegree;
  System.out.println("Azimuth: " + df.format(azimuthDegree));
  return azimuthDegree;
}

Finally, you are passing dcLong in as the lat parameter of the azimuth() method; this should be dcLat.

I'd recommend using radians internally throughout, and only converting from and to degrees on input and output. This will help prevent mistakes, and cut down on rounding errors and unnecessary clutter.

libnova odd behavior for 89.5 degrees north latitude

I've contacted the authors of libnova, who have confirmed this is a bug and are working to correct it.

libnova incorrectly assumes bodies are circumpolar if they are above the horizon when due north. This is untrue: https://astronomy.stackexchange.com/q/963

I've written a fix which doesn't assume this, but may still be inaccurate if a body's declination is non-unimodal: https://astronomy.stackexchange.com/questions/962/is-lunar-elevation-at-a-given-location-for-a-given-day-unimodal

I need an equation for equal movement along an ellipse

If you want the real thing

then the planets closer to the primary focus point (center of mass of stellar system ... very close to star) are moving faster so use Kepler's equation here: C++ implementation of mine. Do not forget to check out all the sub-links in that answer you can find there everything you need.

If you want constant speed instead

Then use parametric ellipse equation

x(a)=x0+rx*cos(a)
y(a)=y0+ry*sin(a)

where a is angle <0,2.0*PI> (x0,y0) is the ellipse center and (rx,ry) are the ellipse semi axises (radii).

if a is incremented with constant speed then the area increase is constant so the a is the mean circular angle not the visual on ellipse !!! For more info look here:

• Issue with ellipse angle calculation and point calculation

[edit1] as MartinR pointed out the speed is not constant

so here is approximation with his formula for speed. Ellipse is axis aligned defined by x0,y0,rx,ry (rx>=ry) the perimeter aproximation l:

h=(rx-ry)/(rx+ry); h*=3.0*h; l=M_PI*(rx+ry)*(1.0+(h/(10.0+sqrt(4.0-h))));

if you want to have n chunks of equal sized steps along the perimeter then

l/=n;

initial computations:

double x0,y0,rx,ry,n,l,h;
x0=Form1->ClientWidth>>1;  // center is centered on form
y0=Form1->ClientHeight>>1;
rx=200; // semiaxises rx>=ry !!!
ry=75;
n=40.0; // number of chunks per ellipse (1/speed)
//l=2.0*M_PI*sqrt(0.5*((rx*rx)+(ry*ry)));  // not accurate enough
h=(rx-ry)/(rx+ry); h*=3.0*h; l=M_PI*(rx+ry)*(1.0+(h/(10.0+sqrt(4.0-h)))); // this is more precise
l/=n; // single step size in units,pixels,or whatever

first the slow bruteforce attack (black):

int i;
double a,da,x,y,xx,yy,ll;
a=0.0;
x=x0+rx*cos(a);
y=y0+ry*sin(a);
for (i=n;i>0;i--)
    {
    xx=x; yy=y;
    for (da=a;;)
        {
        a+=0.001;
        x=x0+rx*cos(a);
        y=y0+ry*sin(a);
        ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
        if (ll>=l) break;
        } da=a-da;
    scr->MoveTo(5.0+50.0*a,5.0);
    scr->LineTo(5.0+50.0*a,5.0+300.0*da);

    scr->MoveTo(x0,y0);
    scr->LineTo(xx,yy);
    scr->LineTo(x ,y );
    ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
    scr->TextOutA(0.5*(x+xx)+20.0*cos(a),0.5*(y+yy)+20.0*sin(a),floor(ll));
    }

Now the approximation (Blue):

a=0.0; da=0;
x=x0+rx*cos(a);
y=y0+ry*sin(a);
for (i=n;i>0;i--)
    {
    scr->MoveTo(5.0+50.0*a,5.0+300.0*da);
    xx=rx*sin(a);
    yy=ry*cos(a);
    da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
    scr->LineTo(5.0+50.0*a,5.0+300.0*da);

    xx=x; yy=y;
    x=x0+rx*cos(a);
    y=y0+ry*sin(a);

    scr->MoveTo(x0,y0);
    scr->LineTo(xx,yy);
    scr->LineTo(x ,y );
    ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
    scr->TextOutA(0.5*(x+xx)+40.0*cos(a),0.5*(y+yy)+40.0*sin(a),floor(ll));
    }

This is clean ellipse step (no debug draws)

a=???; // some initial angle
// point on ellipse 
x=x0+rx*cos(a);
y=y0+ry*sin(a);
// next angle by almost constant speed
xx=rx*sin(a);
yy=ry*cos(a);
da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
// next point on ellipse ...
x=x0+rx*cos(a);
y=y0+ry*sin(a);

Here the output of comparison bruteforce and approximation:

[edit2] little precision boost

 a,da=???; // some initial angle and step (last)
 x=x0+rx*cos(a);
 y=y0+ry*sin(a);
 // next angle by almost constant speed
 xx=rx*sin(a+0.5*da); // use half step angle for aproximation ....
 yy=ry*cos(a+0.5*da);
 da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
 // next point on ellipse ...
 x=x0+rx*cos(a);
 y=y0+ry*sin(a);

the half step angle in approximation lead to much closer result to bruteforce attack

---

Related Topics