//
// Version - Dec 7, 2002
//
//*************************************************************************************************
//Assign a few common parameters
//*************************************************************************************************
	degrees_to_radians = 3.14159265/180
	radians_to_degrees = 1/degrees_to_radians
	solar_constant = 1368

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 Inches2mm (inches)
//*************************************************************************************************
//Converts inches of precipitation to mm
//*************************************************************************************************
{	mm = inches * 25.4
	return mm}

function mph2kts (windinmph)
//*************************************************************************************************
//Converts windspeed from mph to knots
//*************************************************************************************************
{	windinknots = windinmph / 1.152
	return windinknots}

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), dewPt (Cel) and station elev (m)
//***********************************************************************************************
{	td = dewpoint
	altitude = elevation *.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 monthlyDepart(monthprecip, dayofmonth, monthofyear) 
//*************************************************************************************************
//Calculates departure from monthly normal precip for any given day of the month
//*************************************************************************************************
{  	if (monthofyear == 1) {average = 72,days = 31}
	if (monthofyear == 2) {average = 45,days = 28}
	if (monthofyear == 3) {average = 55,days = 31}
	if (monthofyear == 4) {average = 52,days = 30}
	if (monthofyear == 5) {average = 53,days = 31}	
	if (monthofyear == 6) {average = 75,days = 30}
	if (monthofyear == 7) {average = 80,days = 31}
	if (monthofyear == 8) {average = 79,days = 31}
	if (monthofyear == 9) {average = 83,days = 30}
	if (monthofyear == 10){average = 83,days = 31}
   	if (monthofyear == 11){average = 90,days = 30}
  	if (monthofyear == 12){average = 78,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 = 117}
   	if (monthofyear == 3) {totalprior = 172}
	if (monthofyear == 4) {totalprior = 224}
	if (monthofyear == 5) {totalprior = 277}	
	if (monthofyear == 6) {totalprior = 352}
	if (monthofyear == 7) {totalprior = 432}
	if (monthofyear == 8) {totalprior = 511}
	if (monthofyear == 9) {totalprior = 594}
	if (monthofyear == 10){totalprior = 677}
   	if (monthofyear == 11){totalprior = 767}
   	if (monthofyear == 12){totalprior = 845}  
   	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 * .3 * (1 - Math.sin(solar_elevation * degrees_to_radians))
	maxradpossible = maxradpossible - attenuation
	if (maxradpossible < 0) {maxradpossible = 0}
	return maxradpossible}