//

// Version - Nov 12, 2002

//

//*************************************************************************************************

//Assign a few common parameters

//*************************************************************************************************

	radianfactor = 3.14159265/180

	todegrees = 1/radianfactor



function setDecimalPoint(value, places)

//*************************************************************************************************

//Does rounding of numbers to desired number of places after decimal point:

//*************************************************************************************************

 {	var factor = Math.pow(10, places);

	var valInt = Math.round(value * factor);

	var valStr = valInt.toString(10); // convert valInt to a string (base 10)

    if (places == 0) { return valStr; }	// no further adjustments are needed

	var len = valStr.length;

	var radix = len - places - 1; // array index location of digit LEFT of future decimal pt.

    var radShift = 0;

    if (radix < 0) {   // will need to add leading zeroes to avoid something like ". 0".

    	radShift = -radix;

    	radix = 0;

    }

	for (var i = 0; i < radShift; i++) {

		valStr = "0" + valStr;	// prepend zeroes to fill missing spaces

	}

 	var intStr = valStr.substring(0, radix + 1); // integer portion of result

	rounding_cutoff = 5 / Math.pow(10,places + 1)

	if (value < 0.0) {if (value > -1 + rounding_cutoff) {intStr = "-0"}}

	var decStr = valStr.substring(radix + 1, len + radShift);

	return intStr + "." + decStr}



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, temp, dewpoint, elevation)

//*************************************************************************************************

//Computes mixing ratio (g/kg) from pressure (mb), temp (Far), dewPt (Far) and station elev (Ft)

//*************************************************************************************************

{	t = Far2Cel(temp)

	td = Far2Cel(dewpoint)

	altitude = (elevation/3.28)*.001 //station elevation in km

	station_pressure = pressure * Math.pow(2.71828, -altitude/8.5) //reduce pressure for station elevation

	vp_sat = 6.112 * Math.pow(10, (7.5 * t) / (237.7 + t))

	vp = 6.112 * Math.pow(10, (7.5 * td) / (237.7 + td))

	sat_mix_ratio = 621.97 * (vp_sat/(station_pressure - vp_sat))

	mix_ratio = (vp/vp_sat) * sat_mix_ratio

	return mix_ratio}





function calcWetbulb(press, tFar, dpFar)

//*************************************************************************************************

//Computes wet bulb temperature (Far) from pressure (mb), temp (Far), and dewpoint (Far)

//*************************************************************************************************

{	t = Far2Cel(tFar) //need temperature in Celsius

	dp = Far2Cel(dpFar) //need dewpoint in Celsius

	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} //disregard windchill at very light wind speeds

	feelslike = tFar //between 50 and 80 apparent temp set equal to actual temp

	if (windchill<50) {feelslike = windchill} //apparent temp below 50 assigned as wind chill

	if (heatIndx>80) {feelslike = heatIndx} //apparent temp above 80 assigned as heat index

	return feelslike}



function monthlyDepart(monthprecip, dayofmonth, monthofyear) 

//*************************************************************************************************

//Calculates departure from monthly normal precip for any given day of the month

//*************************************************************************************************

{  	if (monthofyear == 1) {average = 3.45,days = 31} //Assign normal monthly precip and days per month

	if (monthofyear == 2) {average = 2.76,days = 28}

	if (monthofyear == 3) {average = 3.86,days = 31}

	if (monthofyear == 4) {average = 3.39,days = 30}

	if (monthofyear == 5) {average = 4.30,days = 31}	

	if (monthofyear == 6) {average = 3.42,days = 30}

	if (monthofyear == 7) {average = 3.51,days = 31}

	if (monthofyear == 8) {average = 3.70,days = 31}

	if (monthofyear == 9) {average = 3.69,days = 30}

	if (monthofyear == 10){average = 3.37,days = 31}

   	if (monthofyear == 11){average = 3.60,days = 30}

  	if (monthofyear == 12){average = 2.89,days = 31}    

   	normaltodate = (average/days)*dayofmonth //Compute normal monthly precip to date

  	departmonth = monthprecip-normaltodate //Compute departure from normal

   	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}		//Assign prior month normal precip values

	if (monthofyear == 2) {totalprior = 3.45}

   	if (monthofyear == 3) {totalprior = 6.21}

	if (monthofyear == 4) {totalprior = 10.07}

	if (monthofyear == 5) {totalprior = 13.46}	

	if (monthofyear == 6) {totalprior = 17.76}

	if (monthofyear == 7) {totalprior = 21.18}

	if (monthofyear == 8) {totalprior = 24.69}

	if (monthofyear == 9) {totalprior = 28.39}

	if (monthofyear == 10){totalprior = 32.08}

   	if (monthofyear == 11){totalprior = 35.45}

   	if (monthofyear == 12){totalprior = 39.05}  

   	normaltotal = totalprior+normmonth //compute total normal precip to date

   	departyear = yearprecip-normaltotal //compute departure from normal

   	return departyear} 



	

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

//*************************************************************************************************

{	solarnoon = (sunset + sunrise)/2 

	sun_time = (solarnoon - time)/60

	sun_time = 12 - sun_time

	return sun_time}

	

function solarAltitude(Jday, suntime, latitude)

//*************************************************************************************************

//Calculates solar elevation angle for given hour, minute, and day of the year

//*************************************************************************************************

{	hour_r = (15 * (12 - suntime)) * radianfactor

	decl = 23.45 * Math.sin((Jday + 284) * (360/365) * radianfactor)

	decl_r = decl * radianfactor

	lat_r = latitude * radianfactor

	sin_alt_r = (Math.cos(lat_r) * Math.cos(decl_r) * Math.cos(hour_r)) + (Math.sin(lat_r) * Math.sin(decl_r))

	altitude = Math.asin(sin_alt_r) * todegrees

	return altitude}

	

function solarAzimuth(Jday, suntime, latitude)

//*************************************************************************************************

//Calculates solar elevation angle for given hour, minute, and day of the year

//*************************************************************************************************

{	hour_r = (15 * (12 - suntime)) * radianfactor

	decl = 23.45 * Math.sin((Jday + 284) * (360/365) * radianfactor)

	decl_r = decl * radianfactor

	lat_r = latitude * radianfactor

	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 = Math.atan(x_azm/y_azm) * todegrees

	if (x_azm > 0) {if (y_azm > 0) {azimuth = azimuth}}

	if (x_azm >= 0) {if (y_azm <= 0) {azimuth = 180 + azimuth}}	

	if (x_azm < 0) {if (y_azm < 0) {azimuth = 180 + azimuth}}

	if (x_azm <= 0) {if (y_azm >= 0) {azimuth = 360 + azimuth}}

	return azimuth}



function maxSolar(solar_elevation)

//*************************************************************************************************

//Calculates maximum expected solar radiation for given hour, minute, and day of the year

//*************************************************************************************************

{	solarconstant = 1366 //Top of atmosphere, horizontal surface 

	maxradpossible = solarconstant * Math.sin(solar_elevation * radianfactor) //Assumes no atmospheric attenuation

	attenuation = maxradpossible * .5 * (1 - Math.sin(solar_elevation * radianfactor)) //Increase attenuation as sun angle decreases

	maxradpossible = maxradpossible - attenuation

	if (maxradpossible < 0) {maxradpossible = 0}

	return maxradpossible}