// Convert radian angle to degrees
function radToDeg(angleRad)
{
return (180.0 * angleRad / Math.PI);
}
//*********************************************************************/
// Convert degree angle to radians
function degToRad(angleDeg)
{
return (Math.PI * angleDeg / 180.0);
}
//*********************************************************************/
//***********************************************************************/
//* Name: calcDayOfYear */
//* Type: Function */
//* Purpose: Finds numerical day-of-year from mn, day and lp year info */
//* Arguments: */
//* month: January = 1 */
//* day : 1 - 31 */
//* lpyr : 1 if leap year, 0 if not */
//* Return value: */
//* The numerical day of year */
//***********************************************************************/
function calcDayOfYear(mn, dy, lpyr)
{
var k = (lpyr ? 1 : 2);
var doy = Math.floor((275 * mn)/9) - k * Math.floor((mn + 9)/12) + dy -30;
return doy;
}
//***********************************************************************/
//* Name: calcDayOfWeek */
//* Type: Function */
//* Purpose: Derives weekday from Julian Day */
//* Arguments: */
//* juld : Julian Day */
//* Return value: */
//* String containing name of weekday */
//***********************************************************************/
function calcDayOfWeek(juld)
{
var A = (juld + 1.5) % 7;
var DOW = (A==0)?"Sunday":(A==1)?"Monday":(A==2)?"Tuesday":(A==3)?"Wednesday":(A==4)?"Thursday":(A==5)?"Friday":"Saturday";
return DOW;
}
//***********************************************************************/
//* Name: calcJD */
//* Type: Function */
//* Purpose: Julian day from calendar day */
//* Arguments: */
//* year : 4 digit year */
//* month: January = 1 */
//* day : 1 - 31 */
//* Return value: */
//* The Julian day corresponding to the date */
//* Note: */
//* Number is returned for start of day. Fractional days should be */
//* added later. */
//***********************************************************************/
function calcJD(year, month, day)
{
if (month <= 2) {
year -= 1;
month += 12;
}
var A = Math.floor(year/100);
var B = 2 - A + Math.floor(A/4);
var JD = Math.floor(365.25*(year + 4716)) + Math.floor(30.6001*(month+1)) + day + B - 1524.5;
return JD;
}
//***********************************************************************/
//* Name: calcDateFromJD */
//* Type: Function */
//* Purpose: Calendar date from Julian Day */
//* Arguments: */
//* jd : Julian Day */
//* Return value: */
//* String date in the form DD-MONTHNAME-YYYY */
//* Note: */
//***********************************************************************/
function calcDateFromJD(jd)
{
var z = Math.floor(jd + 0.5);
var f = (jd + 0.5) - z;
if (z < 2299161) {
var A = z;
} else {
alpha = Math.floor((z - 1867216.25)/36524.25);
var A = z + 1 + alpha - Math.floor(alpha/4);
}
var B = A + 1524;
var C = Math.floor((B - 122.1)/365.25);
var D = Math.floor(365.25 * C);
var E = Math.floor((B - D)/30.6001);
var day = B - D - Math.floor(30.6001 * E) + f;
var month = (E < 14) ? E - 1 : E - 13;
var year = (month > 2) ? C - 4716 : C - 4715;
// alert ("date: " + day + "-" + monthList[month-1].name + "-" + year);
return (day + "-" + monthList[month-1].name + "-" + year);
}
//***********************************************************************/
//* Name: calcDayFromJD */
//* Type: Function */
//* Purpose: Calendar day (minus year) from Julian Day */
//* Arguments: */
//* jd : Julian Day */
//* Return value: */
//* String date in the form DD-MONTH */
//***********************************************************************/
function calcDayFromJD(jd)
{
var z = Math.floor(jd + 0.5);
var f = (jd + 0.5) - z;
if (z < 2299161) {
var A = z;
} else {
alpha = Math.floor((z - 1867216.25)/36524.25);
var A = z + 1 + alpha - Math.floor(alpha/4);
}
var B = A + 1524;
var C = Math.floor((B - 122.1)/365.25);
var D = Math.floor(365.25 * C);
var E = Math.floor((B - D)/30.6001);
var day = B - D - Math.floor(30.6001 * E) + f;
var month = (E < 14) ? E - 1 : E - 13;
var year = (month > 2) ? C - 4716 : C - 4715;
return ((day<10 ? "0" : "") + day + monthList[month-1].abbr);
}
//***********************************************************************/
//* Name: calcTimeJulianCent */
//* Type: Function */
//* Purpose: convert Julian Day to centuries since J2000.0. */
//* Arguments: */
//* jd : the Julian Day to convert */
//* Return value: */
//* the T value corresponding to the Julian Day */
//***********************************************************************/
function calcTimeJulianCent(jd)
{
var T = (jd - 2451545.0)/36525.0;
return T;
}
//***********************************************************************/
//* Name: calcJDFromJulianCent */
//* Type: Function */
//* Purpose: convert centuries since J2000.0 to Julian Day. */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* the Julian Day corresponding to the t value */
//***********************************************************************/
function calcJDFromJulianCent(t)
{
var JD = t * 36525.0 + 2451545.0;
return JD;
}
//***********************************************************************/
//* Name: calGeomMeanLongSun */
//* Type: Function */
//* Purpose: calculate the Geometric Mean Longitude of the Sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* the Geometric Mean Longitude of the Sun in degrees */
//***********************************************************************/
function calcGeomMeanLongSun(t)
{
var L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
while(L0 > 360.0)
{
L0 -= 360.0;
}
while(L0 < 0.0)
{
L0 += 360.0;
}
return L0; // in degrees
}
//***********************************************************************/
//* Name: calGeomAnomalySun */
//* Type: Function */
//* Purpose: calculate the Geometric Mean Anomaly of the Sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* the Geometric Mean Anomaly of the Sun in degrees */
//***********************************************************************/
function calcGeomMeanAnomalySun(t)
{
var M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
return M; // in degrees
}
//***********************************************************************/
//* Name: calcEccentricityEarthOrbit */
//* Type: Function */
//* Purpose: calculate the eccentricity of earth's orbit */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* the unitless eccentricity */
//***********************************************************************/
function calcEccentricityEarthOrbit(t)
{
var e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
return e; // unitless
}
//***********************************************************************/
//* Name: calcSunEqOfCenter */
//* Type: Function */
//* Purpose: calculate the equation of center for the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* in degrees */
//***********************************************************************/
function calcSunEqOfCenter(t)
{
var m = calcGeomMeanAnomalySun(t);
var mrad = degToRad(m);
var sinm = Math.sin(mrad);
var sin2m = Math.sin(mrad+mrad);
var sin3m = Math.sin(mrad+mrad+mrad);
var C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
return C; // in degrees
}
//***********************************************************************/
//* Name: calcSunTrueLong */
//* Type: Function */
//* Purpose: calculate the true longitude of the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun's true longitude in degrees */
//***********************************************************************/
function calcSunTrueLong(t)
{
var l0 = calcGeomMeanLongSun(t);
var c = calcSunEqOfCenter(t);
var O = l0 + c;
return O; // in degrees
}
//***********************************************************************/
//* Name: calcSunTrueAnomaly */
//* Type: Function */
//* Purpose: calculate the true anamoly of the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun's true anamoly in degrees */
//***********************************************************************/
function calcSunTrueAnomaly(t)
{
var m = calcGeomMeanAnomalySun(t);
var c = calcSunEqOfCenter(t);
var v = m + c;
return v; // in degrees
}
//***********************************************************************/
//* Name: calcSunRadVector */
//* Type: Function */
//* Purpose: calculate the distance to the sun in AU */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun radius vector in AUs */
//***********************************************************************/
function calcSunRadVector(t)
{
var v = calcSunTrueAnomaly(t);
var e = calcEccentricityEarthOrbit(t);
var R = (1.000001018 * (1 - e * e)) / (1 + e * Math.cos(degToRad(v)));
return R; // in AUs
}
//***********************************************************************/
//* Name: calcSunApparentLong */
//* Type: Function */
//* Purpose: calculate the apparent longitude of the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun's apparent longitude in degrees */
//***********************************************************************/
function calcSunApparentLong(t)
{
var o = calcSunTrueLong(t);
var omega = 125.04 - 1934.136 * t;
var lambda = o - 0.00569 - 0.00478 * Math.sin(degToRad(omega));
return lambda; // in degrees
}
//***********************************************************************/
//* Name: calcMeanObliquityOfEcliptic */
//* Type: Function */
//* Purpose: calculate the mean obliquity of the ecliptic */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* mean obliquity in degrees */
//***********************************************************************/
function calcMeanObliquityOfEcliptic(t)
{
var seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
var e0 = 23.0 + (26.0 + (seconds/60.0))/60.0;
return e0; // in degrees
}
//***********************************************************************/
//* Name: calcObliquityCorrection */
//* Type: Function */
//* Purpose: calculate the corrected obliquity of the ecliptic */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* corrected obliquity in degrees */
//***********************************************************************/
function calcObliquityCorrection(t)
{
var e0 = calcMeanObliquityOfEcliptic(t);
var omega = 125.04 - 1934.136 * t;
var e = e0 + 0.00256 * Math.cos(degToRad(omega));
return e; // in degrees
}
//***********************************************************************/
//* Name: calcSunRtAscension */
//* Type: Function */
//* Purpose: calculate the right ascension of the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun's right ascension in degrees */
//***********************************************************************/
function calcSunRtAscension(t)
{
var e = calcObliquityCorrection(t);
var lambda = calcSunApparentLong(t);
var tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(lambda)));
var tanadenom = (Math.cos(degToRad(lambda)));
var alpha = radToDeg(Math.atan2(tananum, tanadenom));
return alpha; // in degrees
}
//***********************************************************************/
//* Name: calcSunDeclination */
//* Type: Function */
//* Purpose: calculate the declination of the sun */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* sun's declination in degrees */
//***********************************************************************/
function calcSunDeclination(t)
{
var e = calcObliquityCorrection(t);
var lambda = calcSunApparentLong(t);
var sint = Math.sin(degToRad(e)) * Math.sin(degToRad(lambda));
var theta = radToDeg(Math.asin(sint));
return theta; // in degrees
}
//***********************************************************************/
//* Name: calcEquationOfTime */
//* Type: Function */
//* Purpose: calculate the difference between true solar time and mean */
//* solar time */
//* Arguments: */
//* t : number of Julian centuries since J2000.0 */
//* Return value: */
//* equation of time in minutes of time */
//***********************************************************************/
function calcEquationOfTime(t)
{
var epsilon = calcObliquityCorrection(t);
var l0 = calcGeomMeanLongSun(t);
var e = calcEccentricityEarthOrbit(t);
var m = calcGeomMeanAnomalySun(t);
var y = Math.tan(degToRad(epsilon)/2.0);
y *= y;
var sin2l0 = Math.sin(2.0 * degToRad(l0));
var sinm = Math.sin(degToRad(m));
var cos2l0 = Math.cos(2.0 * degToRad(l0));
var sin4l0 = Math.sin(4.0 * degToRad(l0));
var sin2m = Math.sin(2.0 * degToRad(m));
var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
return radToDeg(Etime)*4.0; // in minutes of time
}
// Return the hour angle for the given location, decl, and time of day
function calcHourAngle(time, longitude, eqtime)
{
return (15.0*(time - (longitude/15.0) - (eqtime/60.0)));
// in degrees
}
// Returns decimal latitude from deg, min, sec entered in the form
function getLatitude(latLongForm)
{
var neg = 0;
var degs = parseFloat(latLongForm["latDeg"].value);
if (latLongForm["latDeg"].value.charAt(0) == '-') {
neg = 1;
}
var mins = parseFloat(latLongForm["latMin"].value);
var secs = parseFloat(latLongForm["latSec"].value);
if(neg != 1){
var decLat = degs + (mins / 60) + (secs / 3600);
} else if(neg == 1){
var decLat = degs - (mins / 60) - (secs / 3600);
} else {
return -9999;
}
return decLat;
}
// Returns decimal longitude from the deg, min, sec entered in the form
function getLongitude(latLongForm)
{
var neg = 0;
var degs = parseFloat(latLongForm["lonDeg"].value);
if (latLongForm["lonDeg"].value.charAt(0) == '-') {
neg = 1;
}
var mins = parseFloat(latLongForm["lonMin"].value);
var secs = parseFloat(latLongForm["lonSec"].value);
var decLon = degs + (mins / 60) + (secs / 3600);
if(neg != 1){
var decLon = degs + (mins / 60) + (secs / 3600);
} else if(neg == 1){
var decLon = degs - (mins / 60) - (secs / 3600);
} else {
return -9999;
}
return decLon;
}
// *****************************************************
// *****************************************************
// *********** Main calculation routine **************
// *****************************************************
// *****************************************************
//***********************************************************************/
//* Name: calcSun */
//* Type: Main Function called by form controls */
//* Purpose: calculate solar position for the entered date, time and */
//* location. Results are reported in azimuth and elevation */
//* (in degrees) and cosine of solar zenith angle. */
//* Arguments: */
//* riseSetForm : for displaying results */
//* latLongForm : for reading latitude and longitude data */
//* index : daylight saving yes/no select */
//* index2 : city select index */
//* Return value: */
//* none */
//* (fills riseSetForm text fields with results of calculations) */
//***********************************************************************/
function calcSun(riseSetForm, latLongForm, index, index2)
{
if(index2 != 0)
{
setLatLong(latLongForm, index2);
}
var latitude = getLatitude(latLongForm);
var longitude = getLongitude(latLongForm);
var indexRS = riseSetForm.mos.selectedIndex
if (isValidInput(riseSetForm, indexRS, latLongForm))
{
if((latitude >= -90) && (latitude < -89.8))
{
alert("All latitudes between 89.8 and 90 S\n will be set to -89.8.");
latLongForm["latDeg"].value = -89.8;
latitude = -89.8;
}
if ((latitude <= 90) && (latitude > 89.8))
{
alert("All latitudes between 89.8 and 90 N\n will be set to 89.8.");
latLongForm["latDeg"].value = 89.8;
latitude = 89.8;
}
//***** Get calc date/time
// var julDay = calcJulianDay(indexRS,
// parseFloat(riseSetForm["day"].value),
// isLeapYear(riseSetForm["year"].value));
var zone = parseFloat(latLongForm["hrsToGMT"].value);
var daySavings = YesNo[index].value; // = 0 (no) or 60 (yes)
if(zone > 12 || zone < -12.5)
{
alert("The offset must be between -12.5 and 12. \n Setting \"Off-Set\"=0");
zone = "0";
latLongForm["hrsToGMT"].value = zone;
}
var ss = parseFloat(riseSetForm["secs"].value);
var mm = parseFloat(riseSetForm["mins"].value);
var hh = parseFloat(riseSetForm["hour"].value) - (daySavings/60);
if(riseSetForm.ampm[1].checked)
{
hh += 12;
}
while (hh > 23)
{
hh -= 24;
}
riseSetForm["hour"].value = hh + (daySavings/60);
if (mm > 9)
{
riseSetForm["mins"].value = mm;
}
else
{
riseSetForm["mins"].value = "0" + mm;
}
if (ss > 9)
{
riseSetForm["secs"].value = ss;
}
else
{
riseSetForm["secs"].value = "0" + ss;
}
riseSetForm.ampm[2].checked = 1;
// timenow is GMT time for calculation
timenow = hh + mm/60 + ss/3600 + zone; // in hours since 0Z
//alert("timenow = " + timenow);
var JD = (calcJD(parseFloat(riseSetForm["year"].value), indexRS + 1, parseFloat(riseSetForm["day"].value)));
var dow = calcDayOfWeek(JD);
var doy = calcDayOfYear(indexRS + 1, parseFloat(riseSetForm["day"].value), isLeapYear(riseSetForm["year"].value));
var T = calcTimeJulianCent(JD + timenow/24.0);
//var L0 = calcGeomMeanLongSun(T);
//var M = calcGeomMeanAnomalySun(T);
//var e = calcEccentricityEarthOrbit(T);
//var C = calcSunEqOfCenter(T);
//var O = calcSunTrueLong(T);
//var v = calcSunTrueAnomaly(T);
var R = calcSunRadVector(T);
//var lambda = calcSunApparentLong(T);
//var epsilon0 = calcMeanObliquityOfEcliptic(T);
//var epsilon = calcObliquityCorrection(T);
var alpha = calcSunRtAscension(T);
var theta = calcSunDeclination(T);
var Etime = calcEquationOfTime(T);
var eqTime = Etime;
var solarDec = theta; // in degrees
var earthRadVec = R;
riseSetForm["eqTime"].value = (Math.floor(100*eqTime))/100;
riseSetForm["solarDec"].value = (Math.floor(100*(solarDec)))/100;
var solarTimeFix = eqTime - 4.0 * longitude + 60.0 * zone;
var trueSolarTime = hh * 60.0 + mm + ss/60.0 + solarTimeFix;
// in minutes
while (trueSolarTime > 1440)
{
trueSolarTime -= 1440;
}
//var hourAngle = calcHourAngle(timenow, longitude, eqTime);
var hourAngle = trueSolarTime / 4.0 - 180.0;
// Thanks to Louis Schwarzmayr for finding our error,
// and providing the following 4 lines to fix it:
if (hourAngle < -180)
{
hourAngle += 360.0;
}
// alert ("Hour Angle = " + hourAngle);
var haRad = degToRad(hourAngle);
var csz = Math.sin(degToRad(latitude)) *
Math.sin(degToRad(solarDec)) +
Math.cos(degToRad(latitude)) *
Math.cos(degToRad(solarDec)) * Math.cos(haRad);
if (csz > 1.0)
{
csz = 1.0;
} else if (csz < -1.0)
{
csz = -1.0;
}
var zenith = radToDeg(Math.acos(csz));
var azDenom = ( Math.cos(degToRad(latitude)) *
Math.sin(degToRad(zenith)) );
if (Math.abs(azDenom) > 0.001) {
azRad = (( Math.sin(degToRad(latitude)) *
Math.cos(degToRad(zenith)) ) -
Math.sin(degToRad(solarDec))) / azDenom;
if (Math.abs(azRad) > 1.0) {
if (azRad < 0) {
azRad = -1.0;
} else {
azRad = 1.0;
}
}
var azimuth = 180.0 - radToDeg(Math.acos(azRad));
if (hourAngle > 0.0) {
azimuth = -azimuth;
}
} else {
if (latitude > 0.0) {
azimuth = 180.0;
} else {
azimuth = 0.0;
}
}
if (azimuth < 0.0) {
azimuth += 360.0;
}
exoatmElevation = 90.0 - zenith;
if (exoatmElevation > 85.0) {
refractionCorrection = 0.0;
} else {
te = Math.tan (degToRad(exoatmElevation));
if (exoatmElevation > 5.0) {
refractionCorrection = 58.1 / te - 0.07 / (te*te*te) +
0.000086 / (te*te*te*te*te);
} else if (exoatmElevation > -0.575) {
refractionCorrection = 1735.0 + exoatmElevation *
(-518.2 + exoatmElevation * (103.4 +
exoatmElevation * (-12.79 +
exoatmElevation * 0.711) ) );
} else {
refractionCorrection = -20.774 / te;
}
refractionCorrection = refractionCorrection / 3600.0;
}
solarZen = zenith - refractionCorrection;
if(solarZen < 108.0) { // astronomical twilight
riseSetForm["azimuth"].value = (Math.floor(100*
azimuth))/100;
riseSetForm["elevation"].value = (Math.floor(100*
(90.0 - solarZen)))/100;
if (solarZen < 90.0) {
riseSetForm["coszen"].value =
(Math.floor(10000.0*(Math.cos(degToRad(solarZen)))))/10000.0;
} else {
riseSetForm["coszen"].value = 0.0;
}
} else { // do not report az & el after astro twilight
riseSetForm["azimuth"].value = "dark";
riseSetForm["elevation"].value = "dark";
riseSetForm["coszen"].value = 0.0;
}
//***********Conv lat and long
convLatLong(latLongForm);
} else { // end of IF ISVALIDINPUT
riseSetForm["azimuth"].value = "error";
riseSetForm["elevation"].value = "error";
riseSetForm["eqTime"].value = "error";
riseSetForm["solarDec"].value = "error";
riseSetForm["coszen"].value = "error";
// alert("Invalid Input");
}
}