//=========================================================
//
// WEATHER_FUNCTIONS.JS is a collection of scripts used by the weather pages
//
//=========================================================

//=========================================================
//
// Define latitude/longitude and other constants
//
//=========================================================
	latitude = 41.1586
	longitude = -96.1922
	degrees_to_radians = Math.PI / 180
	radians_to_degrees = 180 / Math.PI
	PI_2 = 2 * Math.PI
	spacer = "&nbsp;&nbsp;&nbsp;&nbsp;"
	small_font = "<font size=\"1\" face=\"Arial\" color=\"#000000\">&nbsp;&nbsp;"

//=========================================================
//
// Function to display popup windows for individual images
//
//=========================================================
	function popUp(Display_URL) {
		URL = "http://gretnawx.net/" + Display_URL
		window.open(URL,"Display_URL","toolbar=0,scrollbars=0,location=0,statusbar=0,menubar=0,resizable=0,width=500,height=420,top=50,left=50")}

//=========================================================
//
// Names for months
//
//=========================================================
	text_month = new Array(13)
	text_month[1]  = "January"
	text_month[2]  = "February"
	text_month[3]  = "March"
	text_month[4]  = "April"
	text_month[5]  = "May"
	text_month[6]  = "June"
	text_month[7]  = "July"
	text_month[8]  = "August"
	text_month[9]  = "September"
	text_month[10] = "October"
	text_month[11] = "November"
	text_month[12] = "December"

//=========================================================
//
// Number of days in each month
//
//=========================================================
	days_in_month = new Array(13)
	days_in_month[1]  = 31
	days_in_month[2]  = 28
	days_in_month[3]  = 31
	days_in_month[4]  = 30
	days_in_month[5]  = 31
	days_in_month[6]  = 30
	days_in_month[7]  = 31
	days_in_month[8]  = 31
	days_in_month[9]  = 30
	days_in_month[10] = 31
	days_in_month[11] = 30
	days_in_month[12] = 31
	
//=========================================================
//
// Names for days of week
//
//=========================================================
	text_day = new Array(8)
	text_day[1] = "Sunday"
	text_day[2] = "Monday"
	text_day[3] = "Tuesday"
	text_day[4] = "Wednesday"
	text_day[5] = "Thursday"
	text_day[6] = "Friday"
	text_day[7] = "Saturday"
	
//=========================================================
//
// Suffix for numeric days of month
//
//=========================================================
	function dateSuffix(numeral) {
	numeral = "" + numeral
	suffix = "th"
	if (numeral < 10 && numeral < 20) {
		if (numeral.substring(0,1) == "1") {suffix = "st"}
		if (numeral.substring(0,1) == "2") {suffix = "nd"}	
		if (numeral.substring(0,1) == "3") {suffix = "rd"}}
	if (numeral > 20 && numeral < 100) {
		if (numeral.substring(1,2) == "1") {suffix = "st"}
		if (numeral.substring(1,2) == "2") {suffix = "nd"}	
		if (numeral.substring(1,2) == "3") {suffix = "rd"}}	
	if (numeral > 100) {
		if (numeral.substring(2,3) == "1") {suffix = "st"}
		if (numeral.substring(2,3) == "2") {suffix = "nd"}	
		if (numeral.substring(2,3) == "3") {suffix = "rd"}}
	if (numeral > 100 && numeral.substring(1,2) == "1") {suffix = "th"}
	return suffix}
	
//=========================================================
//
// Calculates first day of the year (beginning in 2000)
//
//=========================================================
function calcFirstDay(this_year) {
	loop_year = 2000
	first_day = 7
	while (loop_year < this_year) {
		first_day = first_day + 1
		if (loop_year % 4 == 0) {first_day = first_day + 1}
		if ((loop_year % 100 == 0) && (loop_year % 400 != 0)) {first_day = first_day - 1}
		if (first_day > 7) {first_day = first_day - 7}
		loop_year++}
	return first_day}

//=========================================================
//
// Calculates Julian day of the year
//
//=========================================================
function calcJday(current_month, current_day, current_year) {
	i = 1
	days_in_month[2] = 28
	days_in_year = 365
	leap = "False"
	if (current_year % 4 == 0) {leap = "True"}
	if ((current_year % 100 == 0) && (current_year % 400 != 0)) {leap = "False"}
	if (leap == "True") {days_in_year = 366, days_in_month[2] = 29}
	prior_days = 0
	while (i < current_month) {prior_days = prior_days + days_in_month[i], i++}
	Jday = prior_days + current_day
	return Jday}

//=========================================================
//
// Calculates day of week from Julian date
//
//=========================================================
function calcDayOfWeek(jday, yy) {
	week = Math.floor(jday / 7)
	day_in_week = jday - (7 * week) + calcFirstDay(yy) - 1
	if (day_in_week == 0) {day_in_week = 7}
	if (day_in_week > 7) {day_in_week = day_in_week - 7}
	return text_day[day_in_week]}

//=========================================================
//
// Calculates start and end dates for DST 
//
//=========================================================
function calcDSTDates(cur_year) {
	DST_start_month = 3; start_sunday = 2
	DST_end_month = 11; end_sunday = 1
	start_date = 1
	sundays = 0
	loop_end = start_sunday * 7
	counter = 0
	while (counter <= loop_end) {
		date = start_date + counter
		jul_date = calcJday(DST_start_month, date, cur_year)
		day_name = calcDayOfWeek(jul_date, cur_year)
		if (day_name == "Sunday") {sundays = sundays + 1}
		if (sundays == start_sunday) {DST_start_date = date, counter = loop_end}
	counter++}
	sundays = 0
	loop_end = end_sunday * 7
	counter = 0
	while (counter <= loop_end) {
		date = start_date + counter
		jul_date = calcJday(DST_end_month, date, cur_year)
		day_name = calcDayOfWeek(jul_date, cur_year)
		if (day_name == "Sunday") {sundays = sundays + 1}
		if (sundays == end_sunday) {DST_end_date = date, counter = loop_end}
	counter++}
	return}

//=========================================================
//
// Truncates numeric value to specified number of places
//
//=========================================================
function setDecimal(value, places) {
	if (value < 0) {sign = "–"} else {sign = ""}
	value = Math.abs(value)
	factor = Math.pow(10, places)
	valInt = Math.round(value * factor)
	if (places == 0 && valInt == 0) {sign = ""}
	if (places == 0) {return sign + valInt}
	valStr = valInt.toString(10)
	len = valStr.length
	radix = len - places - 1
	radShift = 0
	if (radix < 0) {radShift = -radix, radix = 0}
	for (i = 0; i < radShift; i++) {valStr = "0" + valStr}
	intStr = valStr.substring(0, radix + 1)
	decStr = valStr.substring(radix + 1, len + radShift)
	valStr = sign + intStr + "." + decStr
	return valStr}

//=========================================================
//
// Adds commas to numeric strings over 1000
//
//=========================================================
function addCommas(nStr) {
	nStr += ''
	x = nStr.split('.')
	x1 = x[0]
	x2 = x.length > 1 ? '.' + x[1] : ''
	rgx = /(\d+)(\d{3})/;
	while (rgx.test(x1)) {x1 = x1.replace(rgx, '$1' + ',' + '$2')}
	return x1 + x2}

//=========================================================
//
// Converts Farenheit to Celsius
//
//=========================================================
function Far2Cel(tFar) {
	tCel = (tFar-32) * (5/9)
	return tCel}

//=========================================================
//
// Converts Celsius to Farenheit
//
//=========================================================
function Cel2Far (tCel) {
	tFar = tCel * (9/5) + 32
	return tFar}

//=========================================================
//
// Converts inches to millimeters
//
//=========================================================
function Inches2mm (inches) {
	mm = inches * 25.4
	return mm}

//=========================================================
//
// Converts millibars to inches of mercury
//
//=========================================================
function mb2inches(pressMB) {
	pressIN = pressMB/33.864
	return pressIN}

//=========================================================
//
// Converts miles per hour to knots
//
//=========================================================
function mph2kts(speedMPH) {
	speedKTS = speedMPH*0.869
	return speedKTS}

//=========================================================
//
// Converts miles per hour to kilometers per hour
//
//=========================================================
function mph2kph(speedMPH) {
	speedKPH = speedMPH*1.609344
	return speedKPH}

//=========================================================
//
// Converts kilometers to miles
//
//=========================================================
function km2miles(input_km) {
	miles_distance = input_km / 1.609344
	return miles_distance}

//=========================================================
//
// Computes wet-bulb temperature
//
//=========================================================
function calcWetbulb(pressure_mb, temp_far, dewpoint_far) {
	temp_celsius = Far2Cel(temp_far)
	dewpoint_celsius = Far2Cel(dewpoint_far)
	tmin = Math.min(dewpoint_celsius, temp_celsius)
	tmax = Math.max(dewpoint_celsius, temp_celsius)
	vapor_pressure = 6.112 * Math.pow(10, (7.5 * dewpoint_celsius) / (237.7 + dewpoint_celsius))
	while (true)
	   	{tcur = (tmax + tmin) / 2
	   	vpcur = 6.112 * Math.pow(10, (7.5 * tcur) / (237.7 + tcur))
	   	peq = 0.00066 * (1+0.00155 * tcur) * pressure_mb * (temp_celsius - tcur)
	   	diff = peq - vpcur + vapor_pressure
	   	if (Math.abs(diff) < 0.01) break
	   	if (diff < 0) {tmax = tcur} else {tmin = tcur}}
	wetbulb_far = Cel2Far(tcur)    
	return wetbulb_far}

//=========================================================
//
// Calculates monthly precip departure
//
//=========================================================
function monthlyDepart(monthtotal, dayofmonth, monthofyear)  {
	normaltodate = 0
	i = 1
	while (i <= dayofmonth) {normaltodate = normaltodate + parseFloat(normals_records [monthofyear] [i].substring(6,9)), i++}
	departmonth = monthtotal - normaltodate
	return departmonth}

//=========================================================
//
// Calculates yearly precip departure
//
//=========================================================
function yearlyDepart(monthofyear, yearprecip, norm_month_to_date) {
	normal_prior_rain = 0
	i = 1
	while (i < monthofyear) {normal_prior_rain = normal_prior_rain + monthly_avg_rain[i], i++}
	normaltotal = normal_prior_rain + norm_month_to_date
	departyear = yearprecip - normaltotal
	return departyear} 

//=========================================================
//
// Calculates monthly snowfall departure
//
//=========================================================
function monthlySnowDepart(monthsnow, dayofmonth, monthofyear)  {
	normalsnowtodate = 0
	i = 1
	while (i <= dayofmonth) {normalsnowtodate = normalsnowtodate + parseFloat(normals_records [monthofyear] [i].substring(10,12)), i++}
	departsnow = monthsnow - normalsnowtodate
	return departsnow}

//=========================================================
//
// Calculates yearly snowfall departure
//
//=========================================================
function seasonSnowDepart(monthofyear, seasonsnow, norm_month_to_date) {
	i = 7
	normal_through_dec = 0
	while (i <= 12) {normal_through_dec = normal_through_dec + monthly_avg_snow[i], i++}
	if (monthofyear < 7) {i = 1, normal_prior_snow = normal_through_dec} else {i = 7, normal_prior_snow = 0}
	while (i < monthofyear) {normal_prior_snow = normal_prior_snow + monthly_avg_snow[i], i++}
	normalsnow = normal_prior_snow + norm_month_to_date
	departseason = seasonsnow - normalsnow
	return departseason} 
	
//=========================================================
//
// Calculates solar altitude angle based on lat, date, time
//
//=========================================================
function solarAltitude(sun_declination, suntime, latitude) {
	sun_hour = (15 * suntime) * degrees_to_radians
	decl = sun_declination * degrees_to_radians
	lat = latitude * degrees_to_radians
	sin_alt = (Math.cos(lat) * Math.cos(decl) * Math.cos(sun_hour)) + (Math.sin(lat) * Math.sin(decl))
	altitude_angle = Math.asin(sin_alt) * radians_to_degrees
	return altitude_angle}

//=========================================================
//
// Calculates solar azimuth angle based on lat, date, time
//
//=========================================================
function solarAzimuth(sun_declination, suntime, latitude) {
	sun_hour = (15 * suntime) * degrees_to_radians
	lat = latitude * degrees_to_radians
	decl = sun_declination * degrees_to_radians
	y_azm = (-1 * (Math.cos(sun_hour)) * Math.cos(decl) * Math.sin(lat)) + (Math.cos(lat) * Math.sin(decl))
	x_azm = Math.sin(sun_hour) * Math.cos(decl)
	azimuth_angle = Math.atan(x_azm/y_azm) * radians_to_degrees
	if (x_azm > 0 && y_azm > 0) {azimuth_angle = azimuth_angle}
	if (x_azm >= 0 && y_azm <= 0 || x_azm < 0 && y_azm < 0) {azimuth_angle = 180 + azimuth_angle}
	if (x_azm <= 0 && y_azm >= 0) {azimuth_angle = 360 + azimuth_angle}
	return azimuth_angle}

//=========================================================
//
// Calculates max possible radiation for given sun angle
//
//=========================================================
function maxSolar(solar_elevation) {
	solar_constant = 1366
	solar_attenuation_factor = .22
	max_attenuation = .30
	maxradpossible = solar_constant * Math.sin(solar_elevation * degrees_to_radians) // Top of Atmosphere
	zenith_angle = 90 - solar_elevation
	attenuation_factor = solar_attenuation_factor * (1 / Math.cos(zenith_angle * degrees_to_radians))
	if (zenith_angle > 45.3) {attenuation_factor = max_attenuation}
	maxradpossible = maxradpossible * (1 - attenuation_factor)
	if (maxradpossible < 0) {maxradpossible = 0}
	return maxradpossible}

//=========================================================
//
// Computes moonrise and moonset for a given date
//
//=========================================================
function moontimes(Y, M, D, Time_Offset) {
	A5 = 0
	D5 = 0
	V0 = 0
	C = 0
	S = 0
	Rstring = ""
	Sstring = ""
	V1 = 0
	m = new Array(14)
	m[13] = 1.00021
	m[14] = 60.40974
	n = new Array
	x = new Array
	moon = 0
	moonrise = "00:-1"
	moonset = "00:-1"
	rise = ""
	set = ""
	Day = D
	Month = M
	m[M] = days_in_month[M]
	Year = Y
	mr = 13
	B5 = parseFloat(latitude)
	L5 = parseFloat(longitude)
	H = parseFloat(Time_Offset)
	Y = parseFloat(Y)
	M = parseFloat(M)
	D = parseFloat(D)
	P2 = PI_2
	DR = degrees_to_radians
	R1 = DR
	K1 = 15 * DR * 1.0027379
	P1 = Math.PI
	L5 = L5 / 360
	Z0 = H / 2
	subroutine6(Y, D, M)
	T = F + (J - 2451545)
	TT = T / 36525 + 1
	subroutine2(T, Z0, L5, DR)
	T = T + Z0
	subroutine5(T, TT, P2)
	AX = A5
	DX = D5
	T = T + 1
	subroutine5(T, TT, P2)
	AY = A5
	DY = D5
	if (AY < AX) {AY = AY + P2}
	Z1 = DR * 90.833
	S = Math.sin(B5 * DR)
	C = Math.cos(B5 * DR)
	Z = Math.cos(Z1)
	M8 = 0
	W8 = 0
	A0 = AX
	D0 = DX
	DA = AY - AX
	DD = DY - DX
	for (C0 = 0; C0 <= 23;) 
	   	{P = (C0 + 1) / 24
	   	A2 = AX + P * DA
	   	D2 = DX + P * DD
	   	subroutine3(T0, A0, A2, D0, D2, V0, K1, C0, C, S, Rstring, Z, V1, mr, P1)
	   	A0 = A2
	   	D0 = D2
	   	V0 = V2
	   	hs = A0
	   	if (rise == "") {rise = "00:-1"}
	   	if (set == "") {set = "00:-1"}
	   	if (alt < -.73) {rise = "00:-1", set = "00:-1"}
	   	alt = Math.floor(alt * 100 + .5) / 100
	   	if (alt < 0 && alt > -.73) {alt = 0}
	   	C0++}
	moon = 1
	mr = 14
	Z0 = H / 24
	subroutine6(Y, D, M)
	T = F + (J - 2451545)
	subroutine2(T, Z0, L5, DR)
	T = T + Z0
	subroutine4(T, P2)
	n[1] = A5
	n[2] = D5
	n[3] = R5
	T = T + 0.5
	subroutine4(T, P2)
	n[4] = A5
	n[5] = D5
	n[6] = R5
	T = T + 0.5
	subroutine4(T, P2)
	n[7] = A5
	n[8] = D5
	n[9] = R5
	if (n[4] <= n[1]) {n[4] = n[4] + P2}
	if (n[7] <= n[4]) {n[7] = n[7] + P2}
	Z1 = R1 * (90.567 - 41.685 / n[6])
	S = Math.sin(B5 * R1)
	C = Math.cos(B5 * R1)
	Z = Math.cos(Z1)
	M8 = 0
	W8 = 0
	A0 = n[1]
	D0 = n[2]
	for (var C0 = 0; C0 <= 23;)
	   	{P = (C0 + 1) / 24
	   	F0 = n[1]
	   	F1 = n[4]
	   	F2 = n[7]
	   	subroutine1(F1, F0, F2, P, B)
	   	A2 = F
	   	F0 = n[2]
	   	F1 = n[5]
	   	F2 = n[8]
	   	subroutine1(F1, F0, F2, P, B)
	   	D2 = F
	   	subroutine3(T0, A0, A2, D0, D2, V0, K1, C0, C, S, Rstring, Z, V1, mr, P1)
	   	A0 = A2
	   	D0 = D2
	   	V0 = V2
	   	if (ri == 1) {moonrise = rise}
	   	if (st == 1) {moonset = set}
	   	C0++}
	mrise = moonrise
	mset = moonset}

//=========================================================	 
//  ** moontimes subfunctions **
//=========================================================
function subroutine1(F1, F0, F2, P, B) {
	A = F1 - F0
	B = F2 - F1 - A
	F = F0 + P * (2 * A + B * (2 * P - 1))
	return}

//=========================================================
function subroutine2(T, Z0, L5, DR) {
	T0 = T / 36525
	S = 24110.5 + 8640184.813 * T0
	S = S + 86636.6 * Z0 + 86400 * L5
	S = S / 86400
	S = S - Math.floor(S)
	T0 = S * 360 * DR
	return}

//=========================================================	 
function subroutine3(T0, A0, A2, D0, D2, V0, K1, C0, C, S, Rstring, Z, V1, mr, P1) {
	ri = 0
	st = 0
	L0 = T0 + C0 * K1
	L2 = L0 + K1
	if (mr == 14) {if (A2 < A0) {A2 = A2 + PI_2}}
	H0 = L0 - A0
	H2 = L2 - A2
	H1 = (H2 + H0) / 2
	D1 = (D2 + D0) / 2
	if (moon == 0) {
		dxc = D1 * 57.29578 + .005
	   	dxc = Math.floor(dxc * 100) / 100; 
	   	alt = 90 - latitude + dxc
	   	if (alt > 90) {alt = 180 - alt}}
	if (C0 == 0) {V0 = S * Math.sin(D0) + C * Math.cos(D0) * Math.cos(H0) - Z} 
	V2 = S * Math.sin(D2) + C * Math.cos(D2) * Math.cos(H2) - Z
	VS = V0 * V2
	if (VS > 0) {return}
	V1 = S * Math.sin(D1) + C * Math.cos(D1) * Math.cos(H1) - Z
	A = 2 * V2 - 4 * V1 + 2 * V0
	B = 4 * V1 - 3 * V0 - V2
	D = B * B - 4 * A * V0
	if (D < 0) {return}
	D = Math.sqrt(D)
	if (V2 > 0 && V0 < 0) {ri = 1; M8 = 1}
	if (V2 < 0 && V0 > 0) {st = 1; W8 = 1}
	E = (-B + D) / (2 * A)
	if (E > 1 || E < 0) {E = (-B - D) / (2 * A)} 
	T3 = C0 + E + 1 / 120
	H7 = H0 + E * (H2 - H0)
	N7 = -1 * Math.cos(D1) * Math.sin(H7)
	D7 = C * Math.sin(D1) - S * Math.cos(D1) * Math.cos(H7)
	AZ = Math.atan(N7 / D7) / DR
	H3 = Math.floor(T3)
	M3 = Math.floor((T3 - H3) * 60)
	if (ri == 1) {rise = convertTime(H3, M3)}
	if (st == 1) {set = convertTime(H3, M3)}
	return}

//=========================================================	 
function subroutine4(T, P2) {
	L = 0.606434 + 0.03660110129 * T
	M = 0.374897 + 0.03629164709 * T
	F = 0.259091 + 0.03674819520 * T
	D = 0.827362 + 0.03386319198 * T
	N = 0.347343 - 0.00014709391 * T
	G = 0.993126 + 0.00273777850 * T
	L = L - Math.floor(L)
	M = M - Math.floor(M)
	F = F - Math.floor(F)
	D = D - Math.floor(D)
	N = N - Math.floor(N)
	G = G - Math.floor(G)
	L = L * P2
	M = M * P2
	F = F * P2
	D = D * P2
	N = N * P2
	G = G * P2
	V = 0.39558 * Math.sin(F + N)
	V = V + 0.08200 * Math.sin(F)
	V = V + 0.03257 * Math.sin(M - F - N)
	V = V + 0.01092 * Math.sin(M + F + N)
	V = V + 0.00666 * Math.sin(M - F)
	V = V - 0.00644 * Math.sin(M + F - 2 * D + N)
	V = V - 0.00331 * Math.sin(F - 2 * D + N)
	V = V - 0.00304 * Math.sin(F - 2 * D)
	V = V - 0.00240 * Math.sin(M - F - 2 * D - N)
	V = V + 0.00226 * Math.sin(M + F)
	V = V - 0.00108 * Math.sin(M + F - 2 * D)
	V = V - 0.00079 * Math.sin(F - N)
	V = V + 0.00078 * Math.sin(F + 2 * D + N)
	U = 1 - 0.10828 * Math.cos(M)
	U = U - 0.01880 * Math.cos(M - 2 * D)
	U = U - 0.01479 * Math.cos(2 * D)
	U = U + 0.00181 * Math.cos(2 * M - 2 * D)
	U = U - 0.00147 * Math.cos(2 * M)
	U = U - 0.00105 * Math.cos(2 * D - G)
	U = U - 0.00075 * Math.cos(M - 2 * D + G)
	W = 0.10478 * Math.sin(M)
	W = W - 0.04105 * Math.sin( 2 * F + 2 * N)
	W = W - 0.02130 * Math.sin(M - 2 * D)
	W = W - 0.01779 * Math.sin(2 * F + N)
	W = W + 0.01774 * Math.sin(N)
	W = W + 0.00987 * Math.sin(2 * D)
	W = W - 0.00338 * Math.sin(M - 2 * F - 2 * N)
	W = W - 0.00309 *Math.sin(G)
	W = W - 0.00190 * Math.sin(2 * F)
	W = W - 0.00144 * Math.sin(M + N)
	W = W - 0.00144 * Math.sin(M - 2 * F - N)
	W = W - 0.00113 * Math.sin(M + 2 * F + 2 * N)
	W = W - 0.00094 * Math.sin(M - 2 * D + G)
	W = W - 0.00092 * Math.sin(2 * M - 2 * D)
	S = W / Math.sqrt(U - V * V)
	A5 = L + Math.atan(S / Math.sqrt(1 - S * S))
	S = V / Math.sqrt(U)
	D5 = Math.atan(S / Math.sqrt(1 - S * S))
	R5 = 60.40974 * Math.sqrt(U)
	return}

//=========================================================	
function subroutine5(T, TT, P2) {
	L = .779072 + .00273790931 * T
	G = .993126 + .0027377785 * T
	L = L - Math.floor(L)
	G = G - Math.floor(G)
	L = L * P2
	G = G * P2
	V = .39785 * Math.sin(L)
	V = V - .01000 * Math.sin(L - G)
	V = V + .00333 * Math.sin(L + G)
	V = V - .00021 * TT * Math.sin(L)
	U = 1 - .03349 * Math.cos(G)
	U = U - .00014 * Math.cos(2 * L)
	U = U + .00008 * Math.cos(L)
	W = -.00010 - .04129 * Math.sin(2 * L)
	W = W + .03211 * Math.sin(G)
	W = W + .00104 * Math.sin(2 * L - G)
	W = W - .00035 * Math.sin(2 * L + G)
	W = W - .00008 * TT * Math.sin(G)
	return}

//=========================================================	
function subroutine6(Y, D, M) {
	G = 1
	if (Y < 1583) {G = 0}
	D1 = Math.floor(D)
	F = D - D1 - .5
	J = -1 * Math.floor(7 * (Math.floor((M + 9) / 12) + Y) / 4)
	if (G != 0) {S = (M - 9)
	   if (S < 0) {S = -1}
	   else if (S >= 0) {S = 1}
	   A = Math.abs(M - 9)
	   J3 = Math.floor(Y + S * Math.floor(A / 7))
	   J3 = -1 *Math.floor((Math.floor(J3 / 100) + 1) * 3 / 4)}
	J = J + Math.floor(275 * M / 9) + D1 + G * J3
	J = J + 1721027 + 2 * G + 367 * Y
	if (F < 0) {F = F + 1;  J = J - 1}
	return}

//=========================================================	
function convertTime(H3, M3) {
	if (H3 < 10) {H3 = "0" + H3}
	if (M3 < 10) {M3 = "0" + M3}
	timestring = H3 + ":" + M3
	return timestring}

//=========================================================	 
//  ** End moontimes subfunctions **
//=========================================================
	
//=========================================================
//
// Compute current phase and % illumination of the moon
//
//=========================================================
function calcMoonPhase(moon_year, moon_month, moon_day, moon_hour) {
	moondate = new Date(moon_year, moon_month - 1, moon_day, moon_hour, 00)
	moonmsec = moondate.getTime()
	GMT_time = moonmsec + (moondate.getTimezoneOffset() * 60 * 1000)
	startDate = new Date(1989, 11, 31, 00, 00)
	startMsec = startDate.getTime()
	dmsec = GMT_time - startMsec
	Dmoon  = ((((dmsec /1000) /60) /60) /24)
	Nmoon = Dmoon * (360 / 365.249)
	if (Nmoon > 0) {Nmoon = Nmoon - Math.floor(Math.abs(Nmoon / 360)) * 360} else {Nmoon = Nmoon + (360 + Math.floor(Math.abs(Nmoon / 360)) * 360)}
	Mo = Nmoon + 279.403303 - 282.768422
	if(Mo < 0) {Mo = Mo + 360}
	Ec = 360 * .016713 * Math.sin(Mo * degrees_to_radians) / Math.PI
	lamda = Nmoon + Ec + 279.403303
	if (lamda > 360) {lamda = lamda - 360}
	Lmoon = 13.1763966 * Dmoon + 318.351648
	if (Lmoon > 0) {Lmoon = Lmoon - Math.floor(Math.abs(Lmoon / 360)) * 360} else {Lmoon = Lmoon + (360 + Math.floor(Math.abs(Lmoon / 360)) * 360)}
	Mm = Lmoon - .1114041 * Dmoon - 36.34041
	if (Mm > 0) {Mm = Mm - Math.floor(Math.abs(Mm / 360)) * 360} else {Mm = Mm + (360 + Math.floor(Math.abs(Mm / 360)) * 360)}
	N65 = 318.510107 - .0529539 * Dmoon
	if (N65 > 0) {N65 = N65 - Math.floor(Math.abs(N65 / 360)) * 360} else {N65 = N65 + (360 + Math.floor(Math.abs(N65 / 360)) * 360)}
	Ev = 1.2739 * Math.sin((2 * (Lmoon - lamda) - Mm) * degrees_to_radians)
	Ae = .1858 * Math.sin(Mo * degrees_to_radians)
	A3 = .37 * Math.sin(Mo * degrees_to_radians)
	Mmp = Mm + Ev - Ae - A3
	Ec = 6.2886 * Math.sin(Mmp * degrees_to_radians)
	A4 = .214 * Math.sin((2 * Mmp) * degrees_to_radians)
	lp = Lmoon + Ev + Ec - Ae + A4
	Vmoon = .6583 * Math.sin((2 * (lp - lamda)) * degrees_to_radians)
	lpp = lp + Vmoon
	D67 = lpp - lamda
	percent_visible = .5 * (1 - Math.cos(D67 * degrees_to_radians)) * 100
	trend = (Math.sin(D67 * degrees_to_radians))
	if((percent_visible > 40) && (percent_visible < 60) && (trend > 0)) {MoonPhase = "First Quarter"}
	if((percent_visible > 40) && (percent_visible < 60) && (trend < 0)) {MoonPhase = "Last Quarter"}             
	if((percent_visible > .02) && (percent_visible <= 40) && (trend > 0)) {MoonPhase = "Waxing Crescent"}
	if((percent_visible > .02) && (percent_visible <= 40) && (trend < 0)) {MoonPhase = "Waning Crescent"}
	if((percent_visible >= 60) && (percent_visible < 99.8) && (trend > 0)) {MoonPhase = "Waxing Gibbous"}
	if((percent_visible >= 60) && (percent_visible < 99.8) && (trend < 0)) {MoonPhase = "Waning Gibbous"}
	if(percent_visible <= .02) {MoonPhase = "New Moon"}
	if(percent_visible >= 99.8) {MoonPhase = "Full Moon"}
	return MoonPhase}

//=========================================================
//
// Computes sunrise and sunset for a given date
//
//=========================================================
function sun_rise_set(Sun_Day, Sun_Month, Sun_Year, Event) {
	K = degrees_to_radians
	sh = Math.sin(-K * 0.8333)
	RISE = "False"
	SET = "False"
	if (Event == "Rise") {Start_Hour = 11} // earliest UTC sunrise of year
	if (Event == "Set") {Start_Hour = 22} // earliest UTC sunset of year
	Sun_Hour = Start_Hour
	while (Sun_Hour <= Start_Hour + 4) // Max diff between solstices for rise or set
		{jd=JulDay (Sun_Day, Sun_Month, Sun_Year, Sun_Hour)
		dec = sunDec(jd)
		gha = computeGHA (Sun_Day, Sun_Month, Sun_Year, Sun_Hour)	
		Y0 = Math.sin(K*computeHeight(dec, latitude, longitude, gha)) - sh
		gha = computeGHA (Sun_Day, Sun_Month, Sun_Year, Sun_Hour + 1.0)
		yPlus = Math.sin(K*computeHeight(dec, latitude, longitude, gha)) - sh
		gha = computeGHA (Sun_Day, Sun_Month, Sun_Year, Sun_Hour - 1.0)
		yMinus = Math.sin(K*computeHeight(dec, latitude, longitude, gha)) - sh
		QUAD(yMinus,yPlus)
		switch (NZ){
			case 1:	if (yMinus < 0.0) {UTRISE = Sun_Hour + zero1; RISE = "True"}
				else {UTSET = Sun_Hour + zero1; SET = "True"}; break;
	
			case 2:	if (YE < 0.0) {UTRISE = Sun_Hour + zero2; UTSET = Sun_Hour + zero1}
				else {UTRISE = Sun_Hour + zero1; UTSET = Sun_Hour + zero2}
				RISE = "True"; SET = "True"; break}
	if (Event == "Rise" && RISE == "True") {return}
	if (Event == "Set" && SET == "True") {return}	
	Sun_Hour++}}

//=========================================================	 
//  ** sun_rise_set subfunctions **
//=========================================================
function QUAD(yMinus, yPlus) {
	NZ = 0
	A = 0.5 * (yMinus + yPlus) - Y0
	B = 0.5 * (yPlus - yMinus)
	C = Y0
	XE = -B / (2 * A)
	YE = (A * XE + B)* XE + C
	DIS = B * B - 4.0 * A * C
	if (DIS>=0) 
		{DX = 0.5 * Math.sqrt(DIS) / Math.abs(A)
		zero1 = XE - DX
		zero2 = XE + DX
		if (Math.abs(zero1) <= 1.0) NZ = NZ + 1
		if (Math.abs(zero2) <= 1.0) NZ = NZ + 1
		if (zero1 < -1.0) zero1 = zero2}}

//=========================================================	 
function computeHeight(de, la, lo, gh) {
	lat_K = la * (degrees_to_radians)
	dec_K = de * (degrees_to_radians)
	x = Number(gh) + Number(lo)
	sinHeight = Math.sin(dec_K) * Math.sin(lat_K) + Math.cos(dec_K) * Math.cos(lat_K) * Math.cos(K*x)	
	return Math.asin(sinHeight) / K}

//=========================================================	 
function computeGHA(T, M, J, STD) {
	N = 365 * J + T + 31 * M - 46	 
	if (M < 3) {N = N + Math.floor((J-1)/4)} else {N = N - Math.floor(0.4 * M + 2.3) + Math.floor(J/4)}		 
	P = STD / 24.0
	X = (P + N - 7.22449E5) * 0.98564734 + 279.306
	X = X * (degrees_to_radians)
	XX = -104.55 * Math.sin(X) - 429.266 * Math.cos(X) + 595.63 * Math.sin(2.0 * X) - 2.283 * Math.cos(2.0 * X)
	XX = XX + 4.6 * Math.sin(3.0 * X) + 18.7333 * Math.cos(3.0 * X)
	XX = XX - 13.2 * Math.sin(4.0 * X) - Math.cos(5.0 * X) - Math.sin(5.0 * X) / 3.0 + 0.5 * Math.sin(6.0 * X) + 0.231
	XX = XX / 240.0 + 360.0 * (P + 0.5)
	if (XX > 360.0) XX = XX - 360.0
	return XX}

//=========================================================	 
function JulDay(date, month, year, UT) {
	if (month <= 2) {month = month + 12, year = year - 1}
	B = Math.floor(year / 400.0) - Math.floor(year / 100.0)  + Math.floor(year / 4.0)
	A = 365.0 * year - 679004.0
	jd = A + B + Math.floor(30.6001 * (month + 1)) + date + UT / 24.0
	jd = jd + 2400000.5
	return jd}

//=========================================================	 
function sunDec(jd) {
	cos_eps = 0.917482
	sin_eps = 0.397778				
	T = (jd - 2451545.0) / 36525.0
	M = PI_2 * frac(0.993133 + 99.997361 * T)
	DL = 6893.0 * Math.sin(M) + 72.0 * Math.sin(2.0 * M)
	L = PI_2 * frac(0.7859453 + M / PI_2 + (6191.2 * T + DL) / 1296000)
	SL = Math.sin(L)
	X = Math.cos(L)
	Y = cos_eps * SL
	Z = sin_eps * SL
	R = Math.sqrt(1.0 - Z * Z)
	dec = (360.0 / PI_2) * Math.atan(Z / R)
	if (dec >= 23.435) {dec = 23.45}
	if (dec <= -23.435) {dec = -23.45}
	return dec}

//=========================================================	 
function frac(X) {
	X = X - Math.floor(X)
	if (X < 0) X = X + 1.0
	return X}
//=========================================================	 
//  ** End sun_rise_set subfunctions **
//=========================================================	 	 

//=========================================================
//
// Calculate solstice and equinox dates for years (1900 - 3000)
//
//=========================================================
function calcSeasons(event, astro_year) { 
	if (event == "VE") {i = 1}
	if (event == "SS") {i = 2} 
	if (event == "AE") {i = 3}
	if (event == "WS") {i = 4}
	astro_date_time = calcEquiSol( i, astro_year)
	return astro_date_time}

//=========================================================	 
//  ** calcSeasons subfunctions **
//=========================================================
function INT (n) {return Math.floor(n)}
function POW2(n) {return Math.pow(n,2)}
function POW3(n) {return Math.pow(n,3)}
function POW4(n) {return Math.pow(n,4)}			
function COS(deg) {return Math.cos(deg * Math.PI/180)}

//=========================================================
function calcEquiSol(i, astro_year) {
	k = i - 1
	calcInitial( k, astro_year)
	T = (JDE0 - 2451545.0) / 36525
	W = 35999.373*T - 2.47
	dL = 1 + 0.0334*COS(W) + 0.0007*COS(2*W)
	S = periodic24(T)
	JDE = JDE0 + ((0.00001*S) / dL)
	TDT = fromJDtoUTC(JDE)
	UTC = fromTDTtoUTC(TDT).toUTCString()
	astro_day = parseFloat(UTC.substring(5,7))
	astro_hour = parseFloat(UTC.substring(17,19))
	astro_min = parseFloat(UTC.substring(20,22))
	astro_sec = parseFloat(UTC.substring(23,25))
	if (astro_sec >= 30) {astro_min = astro_min + 1}
	if (astro_min == 60) {astro_hour = astro_hour + 1}
	if (astro_hour == 24) {astro_day = astro_day + 1}
	if (astro_hour < 10) {astro_hour = "0" + astro_hour}
	if (astro_min < 10) {astro_min = "0" + astro_min}
	astro_time = astro_day + " " + astro_hour + " " + astro_min
	return astro_time}

//=========================================================
function calcInitial( k, astro_year ) { // Valid for years 1000 to 3000
	JDE0=0
	Y=(astro_year-2000)/1000
	switch(k) {
		case 0: JDE0 = 2451623.80984 + 365242.37404*Y + 0.05169*POW2(Y) - 0.00411*POW3(Y) - 0.00057*POW4(Y); break
		case 1: JDE0 = 2451716.56767 + 365241.62603*Y + 0.00325*POW2(Y) + 0.00888*POW3(Y) - 0.00030*POW4(Y); break
		case 2: JDE0 = 2451810.21715 + 365242.01767*Y - 0.11575*POW2(Y) + 0.00337*POW3(Y) + 0.00078*POW4(Y); break
		case 3: JDE0 = 2451900.05952 + 365242.74049*Y - 0.06223*POW2(Y) - 0.00823*POW3(Y) + 0.00032*POW4(Y); break}
	return JDE0} 

//=========================================================
function periodic24(T) {
	A = new Array(485,203,199,182,156,136,77,74,70,58,52,50,45,44,29,18,17,16,14,12,12,12,9,8)
	B = new Array(324.96,337.23,342.08,27.85,73.14,171.52,222.54,296.72,243.58,119.81,297.17,21.02,
			247.54,325.15,60.93,155.12,288.79,198.04,199.76,95.39,287.11,320.81,227.73,15.45)
	C = new Array(1934.136,32964.467,20.186,445267.112,45036.886,22518.443,
			65928.934,3034.906,9037.513,33718.147,150.678,2281.226,
			29929.562,31555.956,4443.417,67555.328,4562.452,62894.029,
			31436.921,14577.848,31931.756,34777.259,1222.114,16859.074)
	S = 0
	for (i=0; i < 24; i++) {S += A[i]*COS( B[i] + (C[i]*T))}
	return S} 

//=========================================================
function fromJDtoUTC(JD) {
    	Z = INT(JD + 0.5)
    	F = (JD + 0.5) - Z
    	if (Z < 2299161) {A = Z}
	else {alpha = INT( (Z-1867216.25) / 36524.25), A = Z + 1 + alpha - INT(alpha / 4)}
    	B = A + 1524
    	C = INT((B-122.1) / 365.25)
    	D = INT(365.25*C)
    	E = INT((B-D )/30.6001)
    	DT = B - D - INT(30.6001*E) + F
    	Mon = E - (E<13.5?1:13)
    	Yr  = C - (Mon>2.5?4716:4715)  
    	Day = INT(DT)
    	H = 24*(DT - Day)
    	Hr = INT(H) 
    	M = 60*(H - Hr)
    	Min = INT(M)
    	Sec = INT(60*(M-Min))		
    	theDate = new Date(0)
    	theDate.setUTCFullYear(Yr, Mon-1, Day)
    	theDate.setUTCHours(Hr, Min, Sec)
    	return(theDate)}

//=========================================================
function fromTDTtoUTC(tobj) {
	TBLfirst = 1620, TBLlast = 2002
	TBL = new Array(
		/*1620*/ 121,112,103, 95, 88, 82, 77, 72, 68, 63, 60, 56, 53, 51,  48, 46, 44, 42, 40, 38,
		/*1660*/  35, 33, 31, 29, 26, 24, 22, 20, 18, 16, 14, 12, 11, 10,   9,  8,  7,  7,  7,  7,
		/*1700*/   7,  7, 8,   8,  9,  9,  9,  9,  9, 10, 10, 10, 10, 10,  10, 10, 10, 11, 11, 11,
		/*1740*/  11, 11, 12, 12, 12, 12, 13, 13, 13, 14,  14, 14, 14, 15, 15, 15, 15, 15, 16, 16,
		/*1780*/  16, 16, 16, 16, 16, 16, 15, 15, 14, 13,  
		/*1800*/ 13.1, 12.5, 12.2, 12.0, 12.0,  12.0, 12.0, 12.0, 12.0, 11.9,  11.6, 11.0, 10.2,  9.2,  8.2,
		/*1830*/  7.1,  6.2,  5.6,  5.4,  5.3,   5.4,  5.6,  5.9,  6.2,  6.5,   6.8,  7.1,  7.3,  7.5,  7.6,
		/*1860*/  7.7,  7.3,  6.2,  5.2,  2.7,   1.4, -1.2, -2.8, -3.8, -4.8,  -5.5, -5.3, -5.6, -5.7, -5.9,
		/*1890*/ -6.0, -6.3, -6.5, -6.2, -4.7,  -2.8, -0.1,  2.6,  5.3,  7.7,  10.4, 13.3, 16.0, 18.2, 20.2,
		/*1920*/ 21.1, 22.4, 23.5, 23.8, 24.3,  24.0, 23.9, 23.9, 23.7, 24.0,  24.3, 25.3, 26.2, 27.3, 28.2,
		/*1950*/ 29.1, 30.0, 30.7, 31.4, 32.2,  33.1, 34.0, 35.0, 36.5, 38.3,  40.2, 42.2, 44.5, 46.5, 48.5,
		/*1980*/ 50.5, 52.5, 53.8, 54.9, 55.8,  56.9, 58.3, 60.0, 61.6, 63.0,  63.8, 64.3)
		// Delta T for 2000 thru 2002 from NASA
	deltaT = 0
	Year = tobj.getUTCFullYear()
	t = (Year - 2000) / 100	// Centuries from the epoch 2000.0
	if (Year >= TBLfirst && Year <= TBLlast ) { // Find correction in table
		if (Year%2) {deltaT = ( TBL[(Year-TBLfirst-1)/2] + TBL[(Year-TBLfirst+1)/2] ) / 2}
		else {deltaT = TBL[(Year-TBLfirst)/2]}} 
	else if (Year < 948) {deltaT = 2177 + 497*t + 44.1*POW2(t)}
	else if (Year >= 948) {deltaT =  102 + 102*t + 25.3*POW2(t)
		if (Year >= 2000 && Year <= 2100) {deltaT += 0.37 * (Year - 2100)}}
	return(new Date(tobj.getTime() - (deltaT * 1000)))}

//=========================================================	 
//  ** End calcSeasons subfunctions **
//=========================================================

//=========================================================
//
// Computes earth-moon distance (in kilometers)
//
//=========================================================
function MoonDistance(moon_days) {
	T = moon_days / 36525
	T2 = T * T
	T3 = T2 * T
	T4 = T3 * T
	D = angle (297.8502042
	     + 445267.1115168 * T
	     - 0.0016300 * T2
	     + T3 / 545868
	     + T4 / 113065000)
	M = angle (357.5291092
	     + 35999.0502909 * T
	     - 0.0001536 * T2
	     + T3 / 24490000)
	Mprime = angle (134.9634114
	     + 477198.8676313 * T
	     + 0.0089970 * T2
	     - T3 / 3536000
	     + T4 / 14712000)
	F = angle(93.2720993 + 483202.0175273 * T - .0034029 * T2 - T3 / 3526000 + T4 / 863310000)

        DcA = new Array(0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0,
                        4, 0, 4, 2, 2, 1, 1, 2, 2, 4, 2, 0, 2, 2, 1, 2,
                        0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2, 2, 2, 4, 0,
                        4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2, 2)

        McA = new Array(0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0,
                        0, 0, 0, 1, 1, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2,
                        1, 2, -2, 0, 0, -1, 0, 0, 1, -1, 2, 2, 1, -1, 0,
                        0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0, 0)

        MpcA = new Array(1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1,
                         1, -1, 3, -2, -1, 0, -1, 0, 1, 2, 0, -3, -2,
                        -1, -2, 1, 0, 2, 0, -1, 1, 0, -1, 2, -1, 1, -2, -1,
			-1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3, 2, 1, -1, 3, -1)

        FcA = new Array(0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,
                        0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0, 0, 0, -2, 0, 0,
                        0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, -2)

	ScAR = new Array(-20905355, -3699111, -2955968, -569925,
			48888, -3149, 246158, -152138, -170733, -204586, -129620,
			108743, 104755, 10321, 0, 79661, -34782, -23210, -21636,
			24208, 30824, -8379, -16675, -12831, -10445, -11650, 14403,
			-7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0,
			-3958, 0, 3258, 2616, -1897, -2117, 2354, 0, 0, -1423,
			-1117, -1571, -1739, 0, -4421, 0, 0, 0, 0, 1165, 0, 0, 8752)

	E = 1 - .002516 * T - .0000074 * T2
	E2 = E * E
	sigma_r = 0
        for (var i=0; i < DcA.length; ++i) {
		ME = 1
		if (MpcA[i] == 1 || MpcA[i] == -1) {ME = E}
            	else if (MpcA[i] == 2 || MpcA[i] == -2) {ME = E2}
		sigma_r += ScAR[i]  * Math.cos(DcA[i] * D + McA[i] * ME * M + MpcA[i] * Mprime + FcA[i] * F)}
	distance = sigma_r / 1000 + 385000.56
	return distance}

//=========================================================	 
//  ** MoonDistance subfunctions **
//=========================================================
   function day2000(y, m, d, h) {
	greg = y*10000 + m*100 + d
	if (m == 1 || m == 2) {y = y - 1, m = m + 12}
	a = Math.floor(y/100)
	b = 2 - a  + Math.floor(a/4) 
	c = Math.floor(365.25 * y)
	d1 = Math.floor(30.6001 * (m + 1));
	return (b + c + d1 -730550.5 + d + h/24)} 

//=========================================================	 
	function angle(deg) {
 	return range(deg) * degrees_to_radians}

//=========================================================	 
function ipart(x) {
	if (x> 0) {a = Math.floor(x)}
	else {a = Math.ceil(x)}	
	return a}

//=========================================================	 
function range(x) {
	b = x / 360
	a = 360 * (b - ipart(b))
	if (a  < 0) {a = a + 360}
	return a}

//=========================================================	 
//  ** End MoonDistance subfunctions **
//=========================================================

//=========================================================
//
// End WEATHER_FUNCTIONS.JS 
//
//=========================================================