HELP!!!! Any one any good at Maths?

I need to write some code to do the following two things:

1) Given a date and location get back the time period during which it is astronomically dark on that day.

2) Given a date, location, Ra Dec Coord of an abject and a virtual horizon determine the time period that that object is above the virtual horizon.

I have implemented some code that I thought should do this but so far it keeps telling me that the sun neither rises nor sets in Maynooth in June and that Orion is up for most of the non-existent night!

The maths involved here has confused me beyod belief!

Bart.

The maths involved here has confused me beyod belief!

And the more I read on this the more confused I got also. There are algorithims to work out where the sun is, and there are algorithims to work out where that point rises etc. but then there are factors like refraction on the atmosphere, air pressure, temperature etc. etc.

Any chance you can post the code you do have and maybe someone can see an error?

I started writing some code to do this last night, but found it too hard without using libraries that were there already

Cheers,
~Al

I tihnk I found a library that I can use to do this for me ( archive.eso.org/JSky/ ) but I would of course rather find the bug in my code because god knows I spent long enough writing it!

My Code is attached below:

[code:1]
//astronomical.java
//This class contains a number of static methods that implement Astronomical formulae
//part of the vtie.util package
//Bart Busschots
//NUI Maynooth

package vtie.util;

//imports
import java.util.*; //for Dates etc

/**
* This class defines a number of static methods that implement
* astronomical functions. <br />
* <br />
* <strong>WARNING:</strong> This class is currently returning invalid answers.
* Do not use it untill tihs warning is removed from the documentation. Use
* <code>vtie.util.JSkyInterface</code> instead, it is a wrapper for the JSky
* (<a href="http://archive.eso.org/JSky/" target="_blank">http://archive.eso.org/JSky/</a>)
* package.
*
* @author Bart Busschots
*/
public abstract class astronomical{
//static variables
/** A constant for the calculation twilight times representing Sun Rise.*/
public static final double SUN_RISE = -0.833;

/** A constant for the calculation twilight times representing Civil Twilight.*/
public static final double CIVIL_TWILIGHT = -6.0;

/** A constant for the calculation twilight times representing Naultical Twilight.*/
public static final double NAUTICAL_TWILIGHT = -12.0;

/** A constant for the calculation twilight times representing Astronomical Twilight.*/
public static final double ASTRONOMICAL_TWILIGHT = -18.0;

/**
* Returns a Time Window Representing the times it is Astronomically Dark on a given Day in a given location.
*/
public static VTTimeWindow getAstronomicalDarkness(Date d, LatLongCoord loc){
//exctracting the time in UT from the date
Calendar theD = new GregorianCalendar();
theD.setTime(d);
theD.setTimeZone(VTLocale.timeZone);

//get the rise set times for this Day (There is no guarantee BOTH the returned rise and set times
//will be in this UT day)
Date tempRS[] = null;
try{
tempRS = riseSetTimes(d, loc, ASTRONOMICAL_TWILIGHT);
}catch(noEventException nee){
//if it is perpetual day return a time window representing a full day otherwise return an empty one
if(nee.reason == noEventException.PERPETUAL_NIGHT){
return new VTTimeWindow(VTTimePeriod.getStartDay(d), VTTimePeriod.getEndDay(d));
}else{
return new VTTimeWindow();
}
}
Date rise = tempRS[0];
Date set = tempRS[1];

//there are now three posibilities to deal with.

VTTimeWindow ans = new VTTimeWindow();
//Both could be in this Day (in terms of UT)
if(VTTimePeriod.inSameDay(rise, set)){
//there are two time periods of darkness inthis case, one from midnight till Dawn and one from
//sun set till the End time of the day.
ans.addTimePeriod(new VTTimePeriod(VTTimePeriod.getStartDay(d), rise));
ans.addTimePeriod(new VTTimePeriod(set, VTTimePeriod.getEndDay(d)));
return ans;
}else{
//they are in differnet days
//there are now two possibliities, the rise could be in this day or the set could be
if(VTTimePeriod.inSameDay(rise, d)){
//if the rise is in this day get the set time from the day before
//there is only one Time Period of darkness in this case
//this period is from yesterdays Set to todays Rise
Calendar yesterday = new GregorianCalendar();
yesterday.setTime(d);
yesterday.setTimeZone(VTLocale.timeZone);
yesterday.add(Calendar.DATE, -1);
Date yrs[] = null;
try{
yrs = riseSetTimes(yesterday.getTime(), loc, ASTRONOMICAL_TWILIGHT);
}catch(Exception e){
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getAstronomicalDarkness(Date, LatLongCoord):");
System.err.println(" A blank Time Window has been returned.");
return new VTTimeWindow();
}
ans.addTimePeriod(new VTTimePeriod(yrs[1], rise));
return ans;
}
if(VTTimePeriod.inSameDay(set, d)){
//if the set is in this day get the rise from the day after
//there is only one Time Period of Darkness inthis case
//this period is from todays set to tomorrows rise
Calendar tomorrow = new GregorianCalendar();
tomorrow.setTime(d);
tomorrow.setTimeZone(VTLocale.timeZone);
tomorrow.add(Calendar.DATE, 1);
Date trs[] = null;
try{
trs = riseSetTimes(tomorrow.getTime(), loc, ASTRONOMICAL_TWILIGHT);
}catch(Exception e){
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getAstronomicalDarkness(Date, LatLongCoord):");
System.err.println(" A blank Time Window has been returned.");
return new VTTimeWindow();
}
ans.addTimePeriod(new VTTimePeriod(set, trs[0]));
return ans;
}
}

//unless something has gone very funny at this point we should have returned
//print out an error warning and return an empty Time Window (ans)
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getAstronomicalDarkness(Date, LatLongCoord):");
System.err.println(" A blank Time Window has been returned.");

return ans;
}

/**
* Returns the rise and set times of the Sun on a given Date in a given place.
* The third argument is the distance above or below the horizon we are actually interested in.
* If it is set to zero we get Sun Rise (center of disk on Mathematical Horizon), however we can
* also use this method to caluculate the Various twilights as follows:
* <ul>
* <li>Sun Rise: -0.833, astronomical.SUN_RISE (Takes atmospheric refraction into account)</li>
* <li>Civil Twilight: -6.0, astronomical.CIVIL_TWILIGHT</li>
* <li>Nautical Twilight: -12.0, astronomical.NAUTICAL_TWILIGHT</li>
* <li>Astronomical Twilight: -18.0, astronomical.ASTRONOMICAL_TWILIGHT</li>
* </ul>
* * The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/riset.html" target="_blank">http://hotel04.ausys.se/pausch/comp/riset.html</a>.
*/
public static Date[] riseSetTimes(Date x, LatLongCoord loc, double h)throws noEventException{
//First calculate the local noon time by getting the noon time in Greenwich and adding or
//subtracting the appropriate number of hours depending on the Longitude
Calendar tempC = new GregorianCalendar();
tempC.setTime(x);
tempC.setTimeZone(VTLocale.UT);
tempC.set(Calendar.AM_PM, Calendar.PM);
tempC.set(Calendar.HOUR, 0);
tempC.set(Calendar.MINUTE, 0);
tempC.set(Calendar.SECOND, 0);
tempC.set(Calendar.MILLISECOND, 0);
double OffsetHours = loc.getDecimalLongitude() * 15.0;
int MSOffset = (int)(OffsetHours * 3600000.0);
tempC.add(Calendar.MILLISECOND, MSOffset);
Date localNoon = tempC.getTime();

//Calculate the RA and Dec of the Sun at local Noon
double tempD[] = getSunRaDec(localNoon);
double RA = tempD[0];
double Dec = tempD[1];

//calculate the UT at which the sun is in the south
double UT_Sun_in_south = (RA - GMST0(localNoon) - loc.getDecimalLongitude()) / 15.0;

//calculate the Cos of LHA (the Sun's Local Hour Angle)
double Cos_LAH = (Math.sin(Math.toRadians(h)) - Math.sin(Math.toRadians(loc.getDecimalLatitude())) * Math.sin(Math.toRadians(Dec)))/(Math.cos(Math.toRadians(loc.getDecimalLatitude())) * Math.cos(Math.toRadians(Dec)));

//if Cos_LAH > 1 then sun always above horizon so throw an exception
if(Cos_LAH > 1.0) throw new noEventException("The Sun is always above the requested horizon on the given day from the given location.", noEventException.PERPETUAL_DAY);

//if Cos_LAH < -1 then the sun is always below the horizon so throw an exception
if(Cos_LAH < -1.0) throw new noEventException("The Sun is always below the requested horizon on the given day from the given location.", noEventException.PERPETUAL_NIGHT);

//otherwise go ahead and calculate LAH
double LAH = Math.toDegrees(Math.acos(Cos_LAH));

//convert LAH from degrees to hours
LAH /= 15.0;

//convert the UTSun_in_south time to a Date
Calendar UTSis = new GregorianCalendar();
UTSis.setTime(x);
UTSis.setTimeZone(VTLocale.UT);
UTSis.set(Calendar.AM_PM, Calendar.AM);
UTSis.set(Calendar.HOUR, 0);
UTSis.set(Calendar.MINUTE, 0);
UTSis.set(Calendar.SECOND, 0);
UTSis.set(Calendar.MILLISECOND, 0);
double UTSisMS = UT_Sun_in_south * 3600000.0;
UTSis.add(Calendar.MILLISECOND, (int)UTSisMS);
Date UTSisDate = UTSis.getTime();

//get the rise and set times (in UT)
double LAHMS = LAH * 3600000.0;
Calendar rise = new GregorianCalendar();
rise.setTime(UTSisDate);
rise.setTimeZone(VTLocale.UT);
rise.add(Calendar.MILLISECOND, ((int)(-LAHMS)));

Calendar set = new GregorianCalendar();
set.setTime(UTSisDate);
set.setTimeZone(VTLocale.UT);
set.add(Calendar.MILLISECOND, (int)LAHMS);

//return the answers
Date ans[] = new Date[2];
ans[0] = rise.getTime();
ans[1] = set.getTime();
return ans;
}

/**
* Returns the Time Period that a given RaDec Coordinates is above a given Virtual Horizon on a given day.
*
* @param loc The Location of the Observer.
* @param target The Ra Dec Coordinates in question.
* @param virtHor The Height in Degrees of the Virtual Horizonm or the height an object has to rise above
* the true horizon in order to be considered visible.
* @param d The Date in question (in UT).
*/
public static VTTimeWindow getVisibleWindow(LatLongCoord loc, RaDecCoord target, int virtHor, Date d){
//exctracting the time in UT from the date
Calendar theD = new GregorianCalendar();
theD.setTime(d);
theD.setTimeZone(VTLocale.timeZone);

//get the rise set times for this Day (There is no guarantee BOTH the returned rise and set times
//will be in this UT day)
Date tempRS[] = null;
try{
tempRS = riseSetTimes(d, loc, target, virtHor);
}catch(noEventException nee){
//if it is perpetual day return a time window representing a full day otherwise return an empty one
if(nee.reason == noEventException.PERPETUAL_NIGHT){
return new VTTimeWindow(VTTimePeriod.getStartDay(d), VTTimePeriod.getEndDay(d));
}else{
return new VTTimeWindow();
}
}
Date rise = tempRS[0];
Date set = tempRS[1];

//there are now three posibilities to deal with.

VTTimeWindow ans = new VTTimeWindow();
//Both could be in this Day (in terms of UT)
if(VTTimePeriod.inSameDay(rise, set)){
//there are two time periods of darkness inthis case, one from midnight till Dawn and one from
//sun set till the End time of the day.
ans.addTimePeriod(new VTTimePeriod(VTTimePeriod.getStartDay(d), rise));
ans.addTimePeriod(new VTTimePeriod(set, VTTimePeriod.getEndDay(d)));
return ans;
}else{
//they are in differnet days
//there are now two possibliities, the rise could be in this day or the set could be
if(VTTimePeriod.inSameDay(rise, d)){
//if the rise is in this day get the set time from the day before
//there is only one Time Period of darkness in this case
//this period is from yesterdays Set to todays Rise
Calendar yesterday = new GregorianCalendar();
yesterday.setTime(d);
yesterday.setTimeZone(VTLocale.timeZone);
yesterday.add(Calendar.DATE, -1);
Date yrs[] = null;
try{
yrs = riseSetTimes(yesterday.getTime(), loc, target, virtHor);
}catch(Exception e){
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getVisibleWindow(RaDecCoord, int, Date):");
System.err.println(" A blank Time Window has been returned.");
return new VTTimeWindow();
}
ans.addTimePeriod(new VTTimePeriod(yrs[1], rise));
return ans;
}
if(VTTimePeriod.inSameDay(set, d)){
//if the set is in this day get the rise from the day after
//there is only one Time Period of Darkness inthis case
//this period is from todays set to tomorrows rise
Calendar tomorrow = new GregorianCalendar();
tomorrow.setTime(d);
tomorrow.setTimeZone(VTLocale.timeZone);
tomorrow.add(Calendar.DATE, 1);
Date trs[] = null;
try{
trs = riseSetTimes(tomorrow.getTime(), loc, target, virtHor);
}catch(Exception e){
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getVisibleWindow(RaDecCoord, int, Date):");
System.err.println(" A blank Time Window has been returned.");
return new VTTimeWindow();
}
ans.addTimePeriod(new VTTimePeriod(set, trs[0]));
return ans;
}
}

//unless something has gone very funny at this point we should have returned
//print out an error warning and return an empty Time Window (ans)
System.err.print("WARNING: an error occoured in method vtie.util.astronomical.getVisibleWindow(RaDecCoord, int, Date):");
System.err.println(" A blank Time Window has been returned.");

return ans;
}

/**
* Returns the rise and set times of an object (given it's Ra Dec Coordinates) on a given Date in a given place.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/riset.html" target="_blank">http://hotel04.ausys.se/pausch/comp/riset.html</a>.
*
* @param x The Date in question.
* @param loc The Location in question.
* @param taret The Ra Dec Coordinate of the object.
* @param h The hight above the horizon the object must rise to be considered visible (in degrees).
*/
public static Date[] riseSetTimes(Date x, LatLongCoord loc, RaDecCoord target, double h)throws noEventException{
//First calculate the local noon time
Date localNoon = getLocalNoon(x, loc);

//Calculate the RA and Dec of the target at local Noon
double RA = target.getDecimalRADegrees();
double Dec = target.getDecimalDec();

//calculate the UT at which the target is in the south
double UT_Sun_in_south = (RA - GMST0(localNoon) - loc.getDecimalLongitude()) / 15.0;

//calculate the Cos of LHA (the Target's Local Hour Angle)
double Cos_LAH = (Math.sin(Math.toRadians(h)) - Math.sin(Math.toRadians(loc.getDecimalLatitude())) * Math.sin(Math.toRadians(Dec)))/(Math.cos(Math.toRadians(loc.getDecimalLatitude())) * Math.cos(Math.toRadians(Dec)));

//if Cos_LAH > 1 then target always above horizon so throw an exception
if(Cos_LAH > 1.0) throw new noEventException("The Target is always above the requested horizon on the given day from the given location.", noEventException.PERPETUAL_DAY);

//if Cos_LAH < -1 then the target is always below the horizon so throw an exception
if(Cos_LAH < -1.0) throw new noEventException("The Target is always below the requested horizon on the given day from the given location.", noEventException.PERPETUAL_NIGHT);

//otherwise go ahead and calculate LAH
double LAH = Math.toDegrees(Math.acos(Cos_LAH));

//convert LAH from degrees to hours
LAH /= 15.0;

//convert the UTSun_in_south time to a Date
Calendar UTSis = new GregorianCalendar();
UTSis.setTime(x);
UTSis.setTimeZone(VTLocale.UT);
UTSis.set(Calendar.AM_PM, Calendar.AM);
UTSis.set(Calendar.HOUR, 0);
UTSis.set(Calendar.MINUTE, 0);
UTSis.set(Calendar.SECOND, 0);
UTSis.set(Calendar.MILLISECOND, 0);
double UTSisMS = UT_Sun_in_south * 3600000.0;
UTSis.add(Calendar.MILLISECOND, (int)UTSisMS);
Date UTSisDate = UTSis.getTime();

//get the rise and set times (in UT)
double LAHMS = LAH * 3600000.0;
Calendar rise = new GregorianCalendar();
rise.setTime(UTSisDate);
rise.setTimeZone(VTLocale.UT);
rise.add(Calendar.MILLISECOND, ((int)(-LAHMS)));

Calendar set = new GregorianCalendar();
set.setTime(UTSisDate);
set.setTimeZone(VTLocale.UT);
set.add(Calendar.MILLISECOND, (int)LAHMS);

//return the answers
Date ans[] = new Date[2];
ans[0] = rise.getTime();
ans[1] = set.getTime();
return ans;
}

/**
* Returns the Day Number of a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double d(Date x){
//convert x to a Calendar to allow easy extraction of infomration
Calendar tempX = new GregorianCalendar();
tempX.setTime(x);
tempX.setTimeZone(VTLocale.UT);

//temporary variables to make the equation easier to read
int y = tempX.get(Calendar.YEAR);
int m = tempX.get(Calendar.MONTH) + 1;
int D = tempX.get(Calendar.DATE);

//calculating the day number for midnight on this date.
double ans = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + D - 730530;

//adding the time of day to this and returning the answer
return ans += UT(x);
}

/**
* Returns the M orbital element of the Sun on a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double MSun(Date x){
return 168.6562 + 4.0923344368 * d(x);
}

/**
* Returns the w orbital element of the Sun on a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double wSun(Date x){
return 282.9404 + 4.70935E-5 * d(x);
}

/**
* Returns the e orbital element of the Sun on a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double eSun(Date x){
return 0.016709 - 1.151E-9 * d(x);
}

/**
* Returns the Sideral time at Greenwich at 00:00 UT on a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/riset.html" target="_blank">http://hotel04.ausys.se/pausch/comp/riset.html</a>.
*/
public static double GMST0(Date x){
return L(x) + 180.0;
}

/**
* Returns the Sun's mean longitude on a given Date.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/riset.html" target="_blank">http://hotel04.ausys.se/pausch/comp/riset.html</a>.
*/
public static double L(Date x){
return MSun(x) + wSun(x);
}

/**
* Returns the Local Sideral Time given a Date and a Location.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/riset.html" target="_blank">http://hotel04.ausys.se/pausch/comp/riset.html</a>.
*/
public static double LST(Date x, LatLongCoord l){
return GMST0(x) + UT(x)*15.0 + l.getDecimalLongitude();
}

/**
* Returns the Universal Time, expressed in hours+decimals given a Date.
*/
public static double UT(Date x){
//convert x to a Calendar to allow easy extraction of infomration
Calendar tempX = new GregorianCalendar();
tempX.setTime(x);
tempX.setTimeZone(VTLocale.UT);

//doing the calculation
double UTime = (double)tempX.get(Calendar.HOUR);
UTime += ((double)tempX.get(Calendar.MINUTE))/60.0;
UTime += ((double)tempX.get(Calendar.SECOND))/3600.0;
UTime += ((double)tempX.get(Calendar.MILLISECOND)/3600000.0);

return UTime;
}

/**
* Returns the Obliquity of the Ecliptic in Degrees.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double ecl(Date x){
return 23.4393 - 3.563E-7 * d(x);
}

/**
* Returns the Sun's RA & Dec on a given Date.
* The first element of the Array returned is the RA and the second element is the Dec.
* Both are expressed in Degrees.
* The formula used can be found at
* <a href="http://hotel04.ausys.se/pausch/comp/ppcomp.html" target="_blank">http://hotel04.ausys.se/pausch/comp/ppcomp.html</a>.
*/
public static double[] getSunRaDec(Date x){
//First calculate the Eccentric Anomoly E in degrees
double E = MSun(x) + eSun(x)*(180/Math.PI) * Math.sin(Math.toRadians(MSun(x))) * (1.0 + eSun(x) * Math.cos(Math.toRadians(MSun(x))));

//Next we need to get the Sun's distance r and it;s true anomoly v. To do this we must first get xv and yv
//and then convert using the atan2 function.
double xv = Math.cos(Math.toRadians(E)) - eSun(x);
double yv = Math.sqrt(1.0 - Math.pow(eSun(x), 2)) * Math.sin(Math.toRadians(E));

double v = atan2(yv, xv);
double r = Math.sqrt(xv*xv + yv*yv);

//now we must compute the Sun's true logitude
double lonsun = v + wSun(x);

//now convert these to rectangular geocentric coordinates xs ans ys
//(zs is zero as we are always in teh place of the ecliptic)
double xs = r * Math.cos(Math.toRadians(lonsun));
double ys = r * Math.sin(Math.toRadians(lonsun));

//convert these to the ecliptic plane xe, ye and ze.
double xe = xs;
double ye = ys * Math.cos(Math.toRadians(ecl(x)));
double ze = ys * Math.sin(Math.toRadians(ecl(x)));

//finally fromthis we can work out the Sun's RA & Dec
double RA = atan2(ye, xe);
double Dec = atan2(ze, Math.sqrt(xe*xe + ye*ye));

//return the answer
double ans[] = {RA, Dec};
return ans;
}

/**
* Converts an x,y coordinate pair (in degrees) to the correct angle in all four quadrants.
* The result is returned in Degrees.
*/
public static double atan2(double x, double y){
if(x > 0.0){
return Math.toDegrees(Math.atan(y/x));
}
if(x < 0.0){
if(y >= 0.0){
return Math.toDegrees(Math.atan(y/x)) + 180.0;
}else{
return Math.toDegrees(Math.atan(y/x)) - 180.0;
}
}

//otherwise (x == 0)
if(y > 0){
return 90.0;
}
if(y < 0){
return -90.0;
}

//if we have still not returned there is a problem so return NaN
return Double.NaN;
}

/**
* Returns the UT Time of Local Noon for a given Location on a given Date.
*/
public static Date getLocalNoon(Date x, LatLongCoord loc){
//Calculate the local noon time by getting the noon time in Greenwich and adding or
//subtracting the appropriate number of hours depending on the Longitude
Calendar tempC = new GregorianCalendar();
tempC.setTime(x);
tempC.setTimeZone(VTLocale.UT);
tempC.set(Calendar.AM_PM, Calendar.PM);
tempC.set(Calendar.HOUR, 0);
tempC.set(Calendar.MINUTE, 0);
tempC.set(Calendar.SECOND, 0);
tempC.set(Calendar.MILLISECOND, 0);
double OffsetHours = loc.getDecimalLongitude() * 15.0;
int MSOffset = (int)(OffsetHours * 3600000.0);
tempC.add(Calendar.MILLISECOND, MSOffset);
return tempC.getTime();
}

}
[/code:1][/code]

:shock: . Is there another site around that you could see how it is done