//
// Version - Nov 15, 2002
//
//*************************************************************************************************
//Assign a few common parameters
//*************************************************************************************************
	degrees_to_radians = 3.14159265/180
	radians_to_degrees = 1/degrees_to_radians
	solar_constant = 1366

function setDecimal(value, places)
//*************************************************************************************************
//Does rounding of numbers to desired number of places after decimal point:
//*************************************************************************************************
{	if (value < 0) {sign = "-"} else {sign = ""}
	value = Math.abs(value) 
	factor = Math.pow(10, places)
	valInt = Math.round(value * factor)
	if (places == 0) {if (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 (var 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}

function Far2Cel(tFar)
//*************************************************************************************************
//Converts degrees F to degrees C
//*************************************************************************************************
{	tCel = (tFar-32) * (5/9)
	return tCel}

function Cel2Far (tCel)
//*************************************************************************************************
//Converts degrees C to degrees F
//*************************************************************************************************
{	tFar = tCel * (9/5) + 32
	return tFar}

function mb2inches(pressMB)
//*************************************************************************************************
//Converts pressure from millibars to inches of mercury
//*************************************************************************************************
{	pressIN = pressMB/33.864
	return pressIN}

function calcMixRatio(pressure, dewpoint, elevation)
//*************************************************************************************************
//Computes mixing ratio (g/kg) from pressure (mb), temp (Far), dewPt (Far) and station elev (Ft)
//***********************************************************************************************
{	td = Far2Cel(dewpoint)
	altitude = (elevation/3.28)*.001 //station elevation in km
	station_pressure = pressure * Math.pow(2.71828, -altitude/8.5)
	vp = 6.112 * Math.pow(10, (7.5 * td) / (237.7 + td))
	mix_ratio = 621.97 * (vp/(station_pressure - vp))
	return mix_ratio}


function calcWetbulb(press, tFar, dpFar)
//*************************************************************************************************
//Computes wet bulb temperature (Far) from pressure (mb), temp (Far), and dewpoint (Far)
//*************************************************************************************************
{	t = Far2Cel(tFar)
	dp = Far2Cel(dpFar)
	tmin = Math.min(dp,t)
    	tmax = Math.max(dp,t)
    	e = 6.112 * Math.pow(10, (7.5 * dp) / (237.7 + dp))
    	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) * press * (t - tcur)
        	diff = peq - vpcur + e
        	if (Math.abs(diff) < 0.01) break
        	if (diff < 0) tmax = tcur
        else tmin = tcur}
    wetbulbf = Cel2Far(tcur)    
    return wetbulbf}

function apparentTemp(heatIndx, tFar, windMPH)
//*************************************************************************************************
//Computes windchill (2001 formula) and heat index and return appropriate value for display 
//*************************************************************************************************
{	windchill = (35.74 + 0.6215 * tFar - 35.75 * Math.pow(windMPH,0.16) + 0.4275 * tFar * Math.pow(windMPH,0.16))
	if (windMPH<3) {windchill = tFar}
	feelslike = tFar
	if (windchill<50) {feelslike = windchill}
	if (heatIndx>80) {feelslike = heatIndx}
	return feelslike}

function monthlyDepart(monthprecip, dayofmonth, monthofyear) 
//*************************************************************************************************
//Calculates departure from monthly normal precip for any given day of the month
//*************************************************************************************************
{  	if (monthofyear == 1) {average = 0.77,days = 31}
	if (monthofyear == 2) {average = 0.80,days = 28}
	if (monthofyear == 3) {average = 2.13,days = 31}
	if (monthofyear == 4) {average = 2.94,days = 30}
	if (monthofyear == 5) {average = 4.44,days = 31}	
	if (monthofyear == 6) {average = 3.95,days = 30}
	if (monthofyear == 7) {average = 3.86,days = 31}
	if (monthofyear == 8) {average = 3.21,days = 31}
	if (monthofyear == 9) {average = 3.17,days = 30}
	if (monthofyear == 10){average = 2.21,days = 31}
   	if (monthofyear == 11){average = 1.82,days = 30}
  	if (monthofyear == 12){average = 0.92,days = 31}    
   	normaltodate = (average / days) * dayofmonth
  	departmonth = monthprecip - normaltodate
   	return departmonth}

function yearlyDepart(monthofyear, yearprecip, normmonth)
//*************************************************************************************************
//Calculates departure from yearly normal precip for any given day of the year
//*************************************************************************************************
{	if (monthofyear == 1) {totalprior = 0}
	if (monthofyear == 2) {totalprior = 0.77}
   	if (monthofyear == 3) {totalprior = 1.57}
	if (monthofyear == 4) {totalprior = 3.70}
	if (monthofyear == 5) {totalprior = 6.64}	
	if (monthofyear == 6) {totalprior = 11.08}
	if (monthofyear == 7) {totalprior = 15.03}
	if (monthofyear == 8) {totalprior = 18.89}
	if (monthofyear == 9) {totalprior = 22.10}
	if (monthofyear == 10){totalprior = 25.27}
   	if (monthofyear == 11){totalprior = 27.48}
   	if (monthofyear == 12){totalprior = 29.30}  
   	normaltotal = totalprior + normmonth
   	departyear = yearprecip - normaltotal
   	return departyear} 

function monthlySnowDepart(monthsnow, dayofmonth, monthofyear) 
//*************************************************************************************************
//Calculates departure from monthly normal precip for any given day of the month
//*************************************************************************************************
{  	if (monthofyear == 1) {averagesnow = 8.0,days = 31}
	if (monthofyear == 2) {averagesnow = 6.5,days = 28}
	if (monthofyear == 3) {averagesnow = 6.9,days = 31}
	if (monthofyear == 4) {averagesnow = 1.5,days = 30}
	if (monthofyear == 5) {averagesnow = 0.0,days = 31}	
	if (monthofyear == 6) {averagesnow = 0.0,days = 30}
	if (monthofyear == 7) {averagesnow = 0.0,days = 31}
	if (monthofyear == 8) {averagesnow = 0.0,days = 31}
	if (monthofyear == 9) {averagesnow = 0.0,days = 30}
	if (monthofyear == 10){averagesnow = 0.5,days = 31}
   	if (monthofyear == 11){averagesnow = 3.5,days = 30}
  	if (monthofyear == 12){averagesnow = 5.5,days = 31}    
   	normalsnowtodate = (averagesnow / days) * dayofmonth
  	departsnow = monthsnow-normalsnowtodate
   	return departsnow}

function seasonSnowDepart(monthofyear, seasonsnow, normmonth)
//*************************************************************************************************
//Calculates departure from yearly normal precip for any given day of the year
//*************************************************************************************************
{	if (monthofyear == 1) {totalpriorsnow = 9.5}
	if (monthofyear == 2) {totalpriorsnow = 17.5}
   	if (monthofyear == 3) {totalpriorsnow = 24.0}
	if (monthofyear == 4) {totalpriorsnow = 29.9}
	if (monthofyear == 5) {totalpriorsnow = 32.4}	
	if (monthofyear == 6) {totalpriorsnow = 32.4}
	if (monthofyear == 7) {totalpriorsnow = 0.0}
	if (monthofyear == 8) {totalpriorsnow = 0.0}
	if (monthofyear == 9) {totalpriorsnow = 0.0}
	if (monthofyear == 10){totalpriorsnow = 0.0}
   	if (monthofyear == 11){totalpriorsnow = 0.5}
   	if (monthofyear == 12){totalpriorsnow = 4.0}  
   	normalsnow = totalpriorsnow + normmonth
   	departseason = seasonsnow - normalsnow
   	return departseason} 


	
function calcJday(current_month, current_day)
//*************************************************************************************************
//Calculates maximum expected solar radiation for given hour, minute, and day of the year
//*************************************************************************************************
{	if (current_month == 1) {priordays = 0}
	if (current_month == 2) {priordays = 31}
	if (current_month == 3) {priordays = 60}
	if (current_month == 4) {priordays = 91}
	if (current_month == 5) {priordays = 120}
	if (current_month == 6) {priordays = 151}
	if (current_month == 7) {priordays = 181}
	if (current_month == 8) {priordays = 212}
	if (current_month == 9) {priordays = 243}
	if (current_month == 10){priordays = 273}
 	if (current_month == 11){priordays = 303}
	if (current_month == 12){priordays = 334}
	Jday = priordays + current_day
	return Jday}

function sunTime(sunrise, sunset, time)
//*************************************************************************************************
//Calculates time relative to solar noon
//*************************************************************************************************
{	solar_noon = (sunrise + sunset)/2 
	suntime = (solar_noon - time)/60
	suntime = 12 - suntime
	return suntime}
	
function solarAltitude(Julian_day, suntime, latitude)
//*************************************************************************************************
//Calculates solar elevation angle for given hour, minute, and day of the year
//*************************************************************************************************
{	hour_r = (15 * (12 - suntime)) * degrees_to_radians
	decl = 23.45 * Math.sin((Julian_day + 284) * (360/365) * degrees_to_radians)
	decl_r = decl * degrees_to_radians
	lat_r = latitude * degrees_to_radians
	sin_alt_r = (Math.cos(lat_r) * Math.cos(decl_r) * Math.cos(hour_r)) + (Math.sin(lat_r) * Math.sin(decl_r))
	altitude_angle = Math.asin(sin_alt_r) * radians_to_degrees
	return altitude_angle}
	
function solarAzimuth(Julian_day, suntime, latitude)
//*************************************************************************************************
//Calculates solar elevation angle for given hour, minute, and day of the year
//*************************************************************************************************
{	hour_r = (15 * (12 - suntime)) * degrees_to_radians
	decl = 23.5 * Math.sin((Julian_day + 284) * (360/365) * degrees_to_radians)
	decl_r = decl * degrees_to_radians
	lat_r = latitude * degrees_to_radians
	y_azm = (-1 * (Math.cos(hour_r)) * Math.cos(decl_r) * Math.sin(lat_r)) + (Math.cos(lat_r) * Math.sin(decl_r))
	x_azm = Math.sin(hour_r) * Math.cos(decl_r)
	azimuth_angle = Math.atan(x_azm/y_azm) * radians_to_degrees
	if (x_azm > 0) {if (y_azm > 0) {azimuth_angle = azimuth_angle}}
	if (x_azm >= 0) {if (y_azm <= 0) {azimuth_angle = 180 + azimuth_angle}}	
	if (x_azm < 0) {if (y_azm < 0) {azimuth_angle = 180 + azimuth_angle}}
	if (x_azm <= 0) {if (y_azm >= 0) {azimuth_angle = 360 + azimuth_angle}}
	return azimuth_angle}

function maxSolar(solar_elevation)
//*************************************************************************************************
//Calculates maximum expected solar radiation for given hour, minute, and day of the year
//*************************************************************************************************
{	maxradpossible = solar_constant * Math.sin(solar_elevation * degrees_to_radians)
	attenuation = maxradpossible * .5 * (1 - Math.sin(solar_elevation * degrees_to_radians))
	maxradpossible = maxradpossible - attenuation
	if (maxradpossible < 0) {maxradpossible = 0}
	return maxradpossible} 