1
$\begingroup$

if I am not using very odd programming language for Astrophysicists, here is my code for calculating path of sun. I am following this reference

I want to calculate the local time as a function of zenith (e.g. sunset, sunrise, golden hour etc)

import java.time.LocalTime;
import java.util.Calendar;
import static java.lang.Math.acos;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static java.lang.Math.tan;


/**
 * Created by rudra on 17/02/18.
 */

public class PathOfSun {

    static public double getTime() {
        LocalTime loacltime = LocalTime.now();//.getHour();
        double tzone = 5.5; //TODO: placeholder, will be calculated
        //    int localhour = time.getHour();
        int hour = loacltime.getHour();
        int min = loacltime.getMinute();
        int sec = loacltime.getSecond();
        double lat  =  28.6;//latlang.Lat;
        double lang =  77.2;//latlang.Lang;
        Calendar calendar = Calendar.getInstance();
        int dayofyear = calendar.get(Calendar.DAY_OF_YEAR);
        int daysthisyear = calendar.getActualMaximum(Calendar.DAY_OF_YEAR);
        double pi = Math.PI;

        double gamma = (2 * pi / daysthisyear) * (dayofyear - 1 + (hour - 12) / 24);
        double eqtime = 229.18 * (0.000075 + 0.001868 * cos(gamma) - 0.032077 * sin(gamma) -
                0.014615 * cos(2 * gamma) - 0.040849 * sin(2 *gamma));
        double decl = 0.006918 - 0.399912 * cos( (gamma)) + 0.070257 * sin( (gamma)) -
                0.006758 * cos(2 *  (gamma)) + 0.000907 * sin(2 *  (gamma)) -
                0.002697 * cos(3 *  (gamma)) + 0.00148 * sin(3 *  (gamma));
        double toffset = 4 * lang*60 + 60*(tzone) + eqtime;
        double tst = 60 * hour + min + sec / 60.0 + toffset;
        double ha = (tst/(4*60)) - 180;

        // Zenith
        double phi = acos(sin(Math.toRadians(lat)) * sin(decl) +
                      cos(Math.toRadians(lat)) * cos((decl)) * cos(Math.toRadians(ha*60)));
        // Azimuth
        double theta = pi - (acos((sin(Math.toRadians(lat)) * cos(phi) - sin(decl)) / (cos(Math.toRadians(lat)) * sin(phi))));

        // Calculating the hourAngle as in second page
        // At sunrise
        double ha2 = acos(cos(Math.toRadians(90.833))/cos( Math.toRadians(lat))*cos((decl))-
                tan(Math.toRadians(lat))*tan((decl)));

        double stime = (720 - 4*(lang + Math.toDegrees(ha)) - eqtime)/60.0;

        System.out.println("DayOfYear"+ dayofyear+ "  DaysThisYear"+daysthisyear);
        System.out.println("lat "+lat+"  lang "+lang+" ltime "+ loacltime);
        System.out.println("eqtime "+eqtime+"  ha "+ha+" , "+ Math.toRadians(ha)+ "  ha2 " + ha2+
                "\ndecl "+Math.toDegrees(decl)+" azimuth "+ Math.toDegrees(theta)+" Zenith "+Math.toDegrees(phi));
        System.out.println("Time of Sun "+ stime);

        return stime;
    }
}

The problem is its giving me more or less OK value for equationoftime and declination, but azimuth etc are grossly incorrect.

ANy idea where I am going wrong? I am following this page for reference.

$\endgroup$
  • 1
    $\begingroup$ This is more a debugging problem rather than an astronomical one. Split the function up into lots of single purpose functions. test each one. The javascript at the site you link to has examples. $\endgroup$ – James K Feb 18 '18 at 8:47
  • $\begingroup$ Actually the problem is I dont have any idea how hourangle ha should look like. I am getting respectable values for equation of tme, declination etc, but hourangle(ha2) is, as from the example, setting for sunrise is: 84.15 and the stime is 66.16. I dont know if those are correct value for sunrise at my coordinate. $\endgroup$ – BaRud Feb 18 '18 at 9:09
2
$\begingroup$

The NOAA formula sheet expresses some quantities in radians, some in degrees, and some in minutes of time. It does say which is which, but it's easy to miss. A few minor code changes can help expose bugs:

  • Store angles in radians. Trigonometric expressions free of Math.toRadians calls are easier to verify against a reference document. Put that document's URL in a comment so it's easy to find.
  • If you need to store something in other units, make it obvious with a naming convention such as long_deg for longitude in degrees.

The other quirks are the $\cos(\phi)$ and $\cos(180 - \theta)$ on the left hand sides of those equations. $\arccos$ has a restricted range of $[0, \pi]$, which is OK for zenith angle but leads to 180-degree errors in azimuth unless you code east/west cases or derive new formulas where you can use Math.atan2.

The solar hour angle should equal 0 at noon and increase by 1 degree every 4 minutes. Sunset hour angle should be about 90 degrees near an equinox, smaller in winter, larger in summer.

The sunrise/sunset formulas produce times in minutes past midnight UTC, so you'll need to apply a time zone offset to get local times.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.