# Set prec_no and block_no of the weather station.
# They are obtained by the website of the Japan Meteorological Agency. (http://www.data.jma.go.jp/obd/stats/etrn/index.php)
prec_no = 44
block_no = 47662

# Set the latitude and longitude of the weather station.
# They are obtained by the website of the Japan Meteorological Agency. (http://www.data.jma.go.jp/obd/stats/etrn/index.php)
lat_deg = 35
lat_min = 41.4
lat_sec = 0

lon_deg = 139
lon_min = 45.6
lon_sec = 0

# Set the tilt angle of the solar panel in degrees.
#   ex. 6   sun -> 31.0
#       5.5 sun -> 28.8
#       5   sun -> 26.6
#       4.5 sun -> 24.2
#       4   sun -> 21.8
beta = 31.0

# Set the azimuth angle of the solar panel in degrees.
# Zero degrees is due south. Eastbound will be a negative value and westbound will be a positive value.
#   ex. Southeast -> -45.0
#       West      ->  90.0
gamma = 0.0

# Set a start and end date of the calculation.
start_year = 2012
start_month = 9
start_day = 1

end_year = 2012
end_month = 9
end_day = 30




################################################################################
import urllib2
import HTMLParser
import datetime
import math

class DataExtractor(HTMLParser.HTMLParser):
	def __init__(self):
		HTMLParser.HTMLParser.__init__(self)
		self.clear()

	def clear(self):
		self.dataNumber = -1
		self.isTargetData = False
		self.trCounter = 0
		self.tdCounter = 0
		self.results = []
		for i in range(0, 24):
			self.results.append(0.0)

	def classIs(self, attributes, classAttr):
		for attr in attributes:
			if attr[0] == 'class' and attr[1] == classAttr:
				return True
		return False

	def handle_starttag(self, tagname, attributes):
		self.isTargetData = False
		if tagname.lower() == 'tr' and self.classIs(attributes, 'mtx'):
			self.trCounter = self.trCounter + 1
			self.tdCounter = 0
		if tagname.lower() == 'td' and self.classIs(attributes, 'data_0_0'):
			self.tdCounter = self.tdCounter + 1
			self.isTargetData = self.tdCounter == self.dataNumber

	def handle_data(self, data):
		if self.isTargetData:
			n = self.trCounter-3
			str = data
			str = str.replace(']', '')
			str = str.replace(')', '')
			str = str.replace('/', '')
			str = str.replace('#', '')
			str = str.replace('-', '')
			str = str.strip()
			try:
				self.results[n] = float(str)
			except:
				print "Error: ", str.decode('utf-8')
				self.results[n] = 0.0

	def extract(self, htmlsrc, dataNumber):
		self.clear()
		self.dataNumber = dataNumber
		self.feed(htmlsrc)


class Erbs:
	def __init__(self, lat, lon, beta, gamma):
		self.lat = lat
		self.lon = lon
		self.beta = beta
		self.gamma = gamma
		self.results = []

	def compute(self, date, hradiations):
		self.results = []

		lat   = math.radians(self.lat)
		beta  = math.radians(self.beta)
		gamma = math.radians(self.gamma)

		# Number of days of year
		n     = (datetime.date(date.year, date.month, date.day) - datetime.date(date.year, 1, 1) + datetime.timedelta(1)).days
		coeff = math.radians((n-1) * 360.0 / 365.0)

		# Declination
		delta = math.radians(360.0 / 2.0 / math.pi * (0.006918 - 0.399912*math.cos(coeff) + 0.070257*math.sin(coeff) - 0.006758*math.cos(2*coeff) + 0.000907*math.sin(2*coeff) - 0.002697*math.cos(3*coeff) + 0.00148*math.sin(2*coeff)))

		# Equation of time
		ET = (0.0172 + 0.4281*math.cos(coeff) - 7.3515*math.sin(coeff) - 3.3495*math.cos(2*coeff) - 9.3619*math.sin(2*coeff)) / 60

		# Albedo		
		rho = 0.2

		# Solar constant		
		Isc = 1.366

		# Geocentric distance of the sun
		rr = 1.0 / (1.00011 + 0.034221*math.cos(coeff) + 0.00128*math.sin(coeff) + 0.000719*math.cos(2*coeff) + 0.000077*math.sin(2*coeff))**0.5

		for i in range(0, 24):
			# Japan Standard Time
			jst = float(i) + 0.5

			# Hour angle
			omega = math.radians(15.0 * (jst + self.lon/15.0 - 9.0 + ET - 12.0))

			# Solar altitude
			sin_h  = math.sin(delta)*math.sin(lat) + math.cos(delta)*math.cos(lat)*math.cos(omega)
			sin_h2 = sin_h
			if sin_h2 <= 0.05:
				sin_h2 = 0.0

			# Angle of incidence of direct solar radiation to the surface of the solar cell
			cos_theta  = math.sin(delta)*math.sin(lat)*math.cos(beta) - math.sin(delta)*math.cos(lat)*math.sin(beta)*math.cos(gamma) + math.cos(delta)*math.cos(lat)*math.cos(beta)*math.cos(omega) + math.cos(delta)*math.sin(lat)*math.sin(beta)*math.cos(gamma)*math.cos(omega) + math.cos(delta)*math.sin(beta)*math.sin(gamma)*math.sin(omega)
			cos_theta2 = cos_theta
			if cos_theta2 <= 0.0:
				cos_theta2 = 0.0

			# Total solar radiation on a horizontal surface outside the atmosphere
			H0 = Isc * sin_h2 / rr / rr

			# Global solar radiation on a horizontal surface (measured value)
			H = hradiations[i] / 3.6

			HH0 = 0.0
			if H0 != 0.0:
				HH0 = H / H0

			# The amount of diffuse solar radiation on a horizontal surface
			Hd = 0.0
			if   HH0 <= 0.0:  Hd = H
			elif HH0 <  0.22: Hd = (1.0 - 0.09*HH0) * H
			elif HH0 <= 0.8:  Hd = (0.9511 - 0.1604*HH0 + 4.388*HH0*HH0 - 16.638*HH0*HH0*HH0 + 12.366*HH0*HH0*HH0*HH0) * H
			else:             Hd = 0.165 * H

			# The amount of direct solar radiation on a horizontal surface
			Hb = H - Hd

			# The amount of solar radiation on the inclined surface
			hb = Hb * cos_theta2 / sin_h					# Direct
			hr = H  * rho * (1.0 - math.cos(beta)) / 2.0	# Reflection
			hd = Hd * (1.0 + math.cos(beta)) / 2.0			# Diffuse

			self.results.append(hb + hr + hd)


if __name__ == "__main__":
	# Output files
	global_solar_radiation_output_filename     = "%04d%02d%02d-%04d%02d%02d-g.csv" % (start_year, start_month, start_day, end_year, end_month, end_day)
	inclined_surface_radiation_output_filename = "%04d%02d%02d-%04d%02d%02d-i.csv" % (start_year, start_month, start_day, end_year, end_month, end_day)

	fileG = open(global_solar_radiation_output_filename, "w")
	fileI = open(inclined_surface_radiation_output_filename, "w")

	d = datetime.date(start_year, start_month, start_day)

	while True:
		# Download the hourly data from the website of the Japan Meteorological Agency.
		url = 'http://www.data.jma.go.jp/obd/stats/etrn/view/hourly_s1.php?prec_no=%d&block_no=%d&year=%d&month=%d&day=%d&view=' % (prec_no, block_no, d.year, d.month, d.day)
		print '%04d/%02d/%02d %s' % (d.year, d.month, d.day, url)
		source = urllib2.urlopen(url)
		htmlsrc = source.read()
		source.close()

		# Extracct the amount of global solar radiation from the data table
		extractor = DataExtractor()
		extractor.extract(htmlsrc, 11)

		# Output the amount of global solar radiation to the file
		fileG.write("%d/%d/%d" % (d.year, d.month, d.day))
		sum = 0.0
		for res in extractor.results:
			fileG.write(",%.2lf" % res)
			sum = sum + res
		fileG.write(",%.2lf" % sum)
		fileG.write("\n")

		# Calculate the amount of solar radiation on inclined plane
		lat = lat_deg + lat_min/60.0 + lat_sec/3600.0
		lon = lon_deg + lon_min/60.0 + lon_sec/3600.0
		erbs = Erbs(lat, lon, beta, gamma)
		erbs.compute(d, extractor.results)

		# Output the amount of solar radiation on inclined plane to the file
		fileI.write("%d/%d/%d" % (d.year, d.month, d.day))
		sum = 0.0
		for res in erbs.results:
			fileI.write(",%lf" % res)
			sum = sum + res
		fileI.write(",%lf" % sum)
		fileI.write("\n")

		# Break when d will be the end date.
		d = d + datetime.timedelta(1)
		if d > datetime.date(end_year, end_month, end_day):
			break

	fileG.close()
	fileI.close()