"""
Calculates solar flux based upon location (lat/long), topography (slope/aspect),
and clearness (1.0 = clear sky). It's pretty amateur, and currently assumes the
time zone is pacific. Descriptions for each of the functions is provided. Runs on
QGIS on point shapefiles that have Slope and Aspect fields and is in wgs 84 projection
(EPSG:4326). Spits out a .csv file with the total insolation for each month and the year
as well as the Lat/Long for each point in the shapefile. This can then be loaded 
into your GIS platform of choice. Have fun messing with it.

This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see .
    
    If you have questions or comments regarding this code
    please e-mail me at info@solamaps.ca
    
"""

def declin(n):
	""" returns the solar angle of declination in radians
		n = day of the year; i.e. n=1 for Jan 1st
	"""
	temp_dec = 23.45*math.sin((math.pi/180)*(360.0/365.0)*(284.0+n))
	#print "declination =",temp_dec
	return math.radians(temp_dec)


def LST(time,long,n):
	""" returns the local solar time in HH.HH
		time = Clock Time, hr
		long = Local Longitude, degree
		n	 = day of the year; i.e. n = 1 for Jan 1st
	"""
	if (n<300 and n>60): #determine Daylight Saving Time Correction
		DT = 1			 #DT = 1 for daylight savings
	else:
		DT = 0			 #DT = 0 for no daylight saving time
	#print "DT =",DT
	#print 'standard meridian is', time_zone
	
	if long < 0:
		long =-long
	
	
	E = LST_E(n)
	time_zone=LST_SM(long)
	lst = (time - DT)+(1/15.0)*(time_zone - long) + E 	
	
	#print "Local Solar Time is:", lst, "Hr"
	return lst
	
def LST_E(n):
	"""Equation of time, E for use in local solar time calculation
	n	 = day of the year; i.e. n = 1 for Jan 1st
	"""
	B = LST_B(n)
	temp_E = 0.165*math.sin(2*B)-0.126*math.cos(B) - 0.025*math.sin(B)
	
	#print "E =", temp_E
	
	return temp_E
	
def LST_B(n):
	"""calculates the B in equation of time and returns it in radians
		n = day of the year; i.e. n = 1 for Jan 1st
	"""
	temp_B = (math.pi/180.0)*(360.0/364.0)*(n-81)
	#print "B=", math.degrees(temp_B)
	return temp_B
	
def LST_SM(long):
	"""returns the correct Standard Meridian based on longitude
		alaska: 135
		pacific: 120
		mountain: 105
		central: 90
		eastern: 75
		atlantic: 60
	"""
	if long >= 165.0 and long < 180.0:
		temp_SM = 165.0
	elif long >= 150.0 and long <165.0:
		temp_SM = 150.0	
	elif long >= 135.0 and long <150.0:
		temp_SM = 135.0
	elif long >=120.0 and long < 135.0:
		temp_SM = 120.0
	elif long >= 105.0 and long < 120.0:
		temp_SM = 105.0
	elif long >= 90.0 and long < 105.0:
		temp_SM = 90.0
	elif long >=75.0 and long < 90.0:
		temp_SM = 75.0
	elif long >= 60.0 and long < 75.0:
		temp_SM = 60.0
	elif long >= 45.0 and long < 60.0:
		temp_SM = 45.0
	elif long >= 30.0 and long < 45.0:
		temp_SM = -30.0
	elif long >= 15.0 and long < 30.0:
		temp_SM = 15.0
	else:
		temp_SM = 0.0
	
	#print "Standard Meridian=", temp_SM
	return temp_SM
		
def hour_angle(time,long,n):
	"""returns the hour angle in radians
	"""
	lst = LST(time,long,n)
	H = (math.pi/180.0)*15.0*(lst- 12.0)
	#print "lst =", lst
	#print "Hour Angle =", H, 'radians'
	return H

def beta(time, lat, long, n):
	"""returns the solar altitude above the horizontal in radians
	"""
	L = math.radians(lat)
	delta = declin(n)
	H = hour_angle(time, long, n)
		
	first_term = math.cos(L)*math.cos(delta)*math.cos(H)
	second_term = math.sin(L)*math.sin(delta)
	temp_beta = math.asin(first_term + second_term)
	
	#print 'Beta =', temp_beta, 'radians'
	return temp_beta
	
def theta(time, lat, long, n, slope, aspect):
	"""returns the solar incident angle in radians
		slope is in degrees
		aspect is in degrees measured from north
	"""
	temp_beta = beta(time, lat, long,n)
	gamma = math.radians(aspect-180)
	sigma = math.radians(slope)
	
	first_term = math.cos(temp_beta)*math.cos(gamma)*math.sin(sigma)
	second_term = math.sin(temp_beta)*math.cos(sigma)
	theta = math.acos(first_term + second_term)
	
	#print 'Theta =', theta, 'radians'
	return theta
	
def get_A(n):
	"""returns the correct value of A used in normal direct
		irradiation calculation
	"""
	if n < 30:
		A = 1202
	elif n >= 30 and n< 60:
		A = 1187.0
	elif n >= 60 and n< 90:
		A = 1164.0
	elif n >= 90 and n< 120:
		A = 1130.0
	elif n >= 120 and n< 150:
		A = 1106.0
	elif n >= 150 and n< 180:
		A = 1092.0
	elif n >= 180 and n< 210:
		A = 1093.0
	elif n >= 210 and n< 240:
		A = 1107.0
	elif n >= 240 and n< 270:
		A = 1136.0
	elif n >= 270 and n< 300:
		A = 1166.0
	elif n >= 300 and n< 330:
		A = 1190.0
	else:
		A = 1204.0
	#print "A =", A
	return A

def get_B(n):
	"""returns the correct value of B used in normal direct
		irradiation calculation
	"""
	if n < 30:
		B = 0.141
	elif n >= 30 and n< 60:
		B = 0.142
	elif n >= 60 and n< 90:
		B = 0.149
	elif n >= 90 and n< 120:
		B = 0.164
	elif n >= 120 and n< 150:
		B = 0.177
	elif n >= 150 and n< 180:
		B = 0.185
	elif n >= 180 and n< 210:
		B = 0.186
	elif n >= 210 and n< 240:
		B = 0.182
	elif n >= 240 and n< 270:
		B = 0.165
	elif n >= 270 and n< 300:
		B = 0.152
	elif n >= 300 and n< 330:
		B = 0.144
	else:
		B = 0.141
	#print "B =", B
	return B

def get_C(n):
	"""returns the correct value of C used in normal direct
		irradiation calculation
	"""
	if n < 30:
		C = 0.103
	elif n >= 30 and n< 60:
		C = 0.104
	elif n >= 60 and n< 90:
		C = 0.109
	elif n >= 90 and n< 120:
		C = 0.12
	elif n >= 120 and n< 150:
		C = 0.13
	elif n >= 150 and n< 180:
		C = 0.137
	elif n >= 180 and n< 210:
		C = 0.138
	elif n >= 210 and n< 240:
		C = 0.134
	elif n >= 240 and n< 270:
		C = 0.121
	elif n >= 270 and n< 300:
		C = 0.111
	elif n >= 300 and n< 330:
		C = 0.106
	else:
		C = 0.103
	#print "C =", C
	return C

def solar_flux(time, lat, long, n, slope, aspect, C_N, albedo):
	"""calculates the instantaneous direct incidence of solar radiation in 
	Watts per meter squared
	time: time of day in hours
	lat: latitude degrees
	longitude: in degrees
	n: day of the year, n=1 for January 1st
	time_zone: time zone (Pacific = 'pst', Mountain = 'mst', Central = 'cst'
				Eastern = 'est', Atlantic = 'ast', Alaskan = 'ast')
	slope: slope of surface in degrees
	aspect: aspect of surface in degrees, north = 0
	C_N: sky clearness number, C_N = 1 for clear sky
	albedo: reflectiveness of the surrounding surfaces
	"""
	A = get_A(n)
	B = get_B(n)
	C = get_C(n)
	
	temp_beta = beta(time, lat, long, n)
	#print temp_beta
	if temp_beta >0:
		
		if abs(temp_beta) < 0.001:
			temp_beta = 0.001
	
		
		temp_theta = theta(time, lat, long, n, slope, aspect)
		#print math.sin(temp_beta)
	
		G_ND = A / (math.exp(B/math.sin(temp_beta)))
		temp_G_D = C_N * G_ND * math.cos(temp_theta)
		temp_G_d = C * G_ND * ((1+math.cos(math.radians(slope)))/2)
		temp_G_R = G_ND * (math.sin(temp_beta)+C)*albedo*((1-math.cos(math.radians(slope)))/2)
	
		if temp_G_D > 0:
			G_D = temp_G_D
		else:
			G_D = 0.0
		
		if temp_G_d > 0:
			G_d = temp_G_d
		else:
			G_d = 0.0
	
		if temp_G_R > 0:
			G_R = temp_G_R
		else:
			G_R = 0.0	
		
		G_t = G_D + G_d + G_R
	
		#print "G_D =", G_D, "G_d =", G_d, "G_R =", G_R
		#print "G_t =", G_t
		return G_t
	else:
		G_t = 0.0
		#print "G_t =", G_t
		return G_t
	
	
def daily_total(lat, long, n, slope, aspect, C_N, albedo):
	"""adds the total daily flux to calculate the total daily insolation
		and returns the result.
	"""
	t = 1
	sum = 0
	photo_hr = 0.0
	while t <= 24:
		#print t
		term = solar_flux(t, lat, long, n, slope, aspect, C_N, albedo)
		if term > 0:
			photo_hr +=1
		sum += term
		#print sum/1000
		t += 1
	daily_tot = [sum/1000,photo_hr]
	#print "total for the", n, 'th day is', sum/1000,"kW*hr/m2"
	#print daily_tot
	return daily_tot
	
def annual_total(lat, long, slope, aspect, C_N, albedo):
	"""calculates the daily total insolation for each day of the year
		and adds each value to a list which is the return
	"""
	n = 1
	sum = 0
	pow_table = []
	time_table = []
	while n <= 365:
		term = daily_total(lat, long, n, slope, aspect, C_N, albedo)
		pow_table += [term[0]]
		time_table += [term[1]]
		sum += term[0]
		n += 1
		#print n
	
	#print table
	#print "The total solar radiation for the year is"
	#print sum,"kW*hr/m2"
	#print pow_table, time_table, len(pow_table), len(time_table)
	solar_data = [pow_table,time_table]
	return solar_data
	
def test_totals(lat, long, time_zone, slope, aspect, C_N, albedo):
	"""test the daily totals one function iteratively"""
	n=1
	while n<=365:
		print  "n=", n, "tot=" ,daily_total(lat, long, n, time_zone, slope, aspect, C_N, albedo)
		n +=1
	return
	
def solar_user():
	"""Welcome to Solar Calculator:1.0
		calculates the annual insolation
		in kW*Hr/(yr*m2)squared.
	 	time: time of day in hours
	 	lat: latitude degrees
	 	longitude: in degrees
	 	time_zone: time zone (Pacific = 120, Mountain = 105, Central = 90
	 				Eastern = 75, Atlantic = 60, Alaskan = 135)
	 	slope: slope of surface in degrees
	 	aspect: aspect of surface in degrees, north = 0
	 	C_N: sky clearness number, C_N = 1 for clear sky
	 	albedo: reflectiveness of the surrounding surfaces
	"""
	print "Lat(Degrees Vancouver = 49.20):"
	lat = float(raw_input())
	
	print "Long(Degrees Vancouver = -123.17):"
	temp_long = float(raw_input())
	long = abs(temp_long)*-1
	
	#print "Time Zone(Pacific = 120, Mountain = 105, Central = 90 Eastern = 75, Atlantic = 60, Alaskan = 135):" 
	#temp_zone = float(raw_input())
	time_zone = float(120)
	
	print "Slope (Degrees):"
	slope = float(raw_input())
	
	print "Aspect (Degrees North = 0):"
	aspect = float(raw_input())
	
	print "Clearness Factor (0-1, Clear Sky = 1):"
	C_N = float(raw_input())
	
	print "Surrounding Surface Albedo (0-1):"
	albedo = float(raw_input())
	
	solar_stats(lat, long, time_zone, slope, aspect, C_N, albedo)
	
	return 
	
def solar_stats(lat, long, slope, aspect, C_N, albedo):
	
	dataset = annual_total(lat, long, slope, aspect, C_N, albedo)
	pow_table = dataset[0]
	time_table = dataset[1]
	
	pow_jan = pow_table[:30]
	pow_feb = pow_table[31:59]
	pow_mar = pow_table[60:90]
	pow_apr = pow_table[91:121]
	pow_may = pow_table[122:153]
	pow_jun = pow_table[154:184]
	pow_jul = pow_table[185:216]
	pow_aug = pow_table[217:247]
	pow_sep = pow_table[248:279]
	pow_oct = pow_table[280:310]
	pow_nov = pow_table[311:341]
	pow_dec = pow_table[341:]
	
	
	ann_tot = math.fsum(pow_table)
	ann_ave = ann_tot/len(pow_table)
	
	jan_tot = math.fsum(pow_jan)
	jan_ave = jan_tot/len(pow_jan)
	
	feb_tot = math.fsum(pow_feb)
	feb_ave = feb_tot/len(pow_feb)
	
	mar_tot = math.fsum(pow_mar)
	mar_ave = mar_tot/len(pow_mar)
	
	apr_tot = math.fsum(pow_apr)
	apr_ave = apr_tot/len(pow_apr)
	
	may_tot = math.fsum(pow_may)
	may_ave = may_tot/len(pow_may)
	
	jun_tot = math.fsum(pow_jun)
	jun_ave = jun_tot/len(pow_jun)
	
	jul_tot = math.fsum(pow_jul)
	jul_ave = jul_tot/len(pow_jul)
	
	aug_tot = math.fsum(pow_aug)
	aug_ave = aug_tot/len(pow_aug)
	
	sep_tot = math.fsum(pow_aug)
	sep_ave = sep_tot/len(pow_sep)
	
	oct_tot = math.fsum(pow_oct)
	oct_ave = oct_tot/len(pow_oct)
	
	nov_tot = math.fsum(pow_nov)
	nov_ave = nov_tot/len(pow_nov)
	
	dec_tot = math.fsum(pow_dec)
	dec_ave = dec_tot/len(pow_dec)
	
	month_avs = [jan_ave,feb_ave,mar_ave,apr_ave,may_ave,jun_ave,jul_ave,aug_ave,sep_ave,
	oct_ave,nov_ave,dec_ave]
	
	monthly_tots = [jan_tot, feb_tot, mar_tot, apr_tot, may_tot, jun_tot, jul_tot, aug_tot, sep_tot, oct_tot, nov_tot, dec_tot]
	
	time_jan = time_table[:30]
	time_feb = time_table[31:59]
	time_mar = time_table[60:90]
	time_apr = time_table[91:121]
	time_may = time_table[122:153]
	time_jun = time_table[154:184]
	time_jul = time_table[185:216]
	time_aug = time_table[217:247]
	time_sep = time_table[248:279]
	time_oct = time_table[280:310]
	time_nov = time_table[311:341]
	time_dec = time_table[341:]
	
	time_ann_tot = math.fsum(time_table)
	time_ann_ave = time_ann_tot/len(time_table)
	
	time_jan_tot = math.fsum(time_jan)
	time_jan_ave = time_jan_tot/len(time_jan)
	
	time_feb_tot = math.fsum(time_feb)
	time_feb_ave = time_feb_tot/len(time_feb)
	
	time_mar_tot = math.fsum(time_mar)
	time_mar_ave = time_mar_tot/len(time_mar)
	
	time_apr_tot = math.fsum(time_apr)
	time_apr_ave = time_apr_tot/len(time_apr)
	
	time_may_tot = math.fsum(time_may)
	time_may_ave = time_may_tot/len(time_may)
	
	time_jun_tot = math.fsum(time_jun)
	time_jun_ave = time_jun_tot/len(time_jun)
	
	time_jul_tot = math.fsum(time_jul)
	time_jul_ave = time_jul_tot/len(time_jul)
	
	time_aug_tot = math.fsum(time_aug)
	time_aug_ave = time_aug_tot/len(time_aug)
	
	time_sep_tot = math.fsum(time_aug)
	time_sep_ave = time_sep_tot/len(time_sep)
	
	time_oct_tot = math.fsum(time_oct)
	time_oct_ave = time_oct_tot/len(time_oct)
	
	time_nov_tot = math.fsum(time_nov)
	time_nov_ave = time_nov_tot/len(time_nov)
	
	time_dec_tot = math.fsum(time_dec)
	time_dec_ave = time_dec_tot/len(time_dec)
	
	time_month_avs = [time_jan_ave,time_feb_ave,time_mar_ave,time_apr_ave,time_may_ave,time_jun_ave,time_jul_ave,time_aug_ave,time_sep_ave,time_oct_ave,time_nov_ave,time_dec_ave]
	return [monthly_tots, time_month_avs]
	
    
    
    
def print_stats(): #this is scaffolding, could be functional again by making it call solar_stats
	print
	print"======================================================================="
	print "Solar Insolation Data for", lat, 'latitude,',long, "longitude"	print "slope =", slope, "aspect=", aspect
	print"======================================================================="
	print 
	print "Total Annual Insolation =",ann_tot,"(kW*hr/m2*year)"
	print
	print "Average Daily Vales (by month)"
	print
	print '          Insolation in (kW*Hr/m2*day) | Photoperiod in (Hrs) | Power Output in (kW/m2)'
	for i in range(len(month_avs)):
		insol = month_avs[i]
		insol_space = ' ' * (27-len(str(insol)))
		time = time_month_avs[i]
		time_space = ' ' * (19-len(str(time)))
		print "Month",i+1,':', insol, insol_space,"|", time, time_space,"|", insol/time		
	print"========================================================================="
	print
	return
	
	
	
def qgis_singlept_func():
    layer = iface.activeLayer()
    layer_name = layer.name()
    print layer_name
    output_file = open('/Your/File/Path/{}_monthly.csv'.format(layer_name) ,'a')
    header ='%s,%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s\n' % ('Latitude', 'Longitude', 'January' , 'February', 'March', 'April' , 'May' , 'June','July','August','September','October','November','December', 'Total')
    unicode_header = header.encode('utf-8')
    output_file.write(unicode_header)
    for f in layer.getFeatures():
        geom = f.geometry()
        #name = f['GGRPHCLNM']
        lat = geom.asPoint().y()
        long = geom.asPoint().x()
        print lat,long
        #time_zone = 120
        slope = f['Slope']
        aspect = f['Aspect']
        C_N = 1.0
        albedo = .33
        
        monthly_data = solar_stats(lat, long, slope, aspect, C_N, albedo)
        monthly_insol = monthly_data[0]
        monthly_time = monthly_data[1]
        
        line = '%f, %f, %f, %f, %f, %f, %f, %f, %f, %f,%f, %f, %f, %f, %f\n' % (lat, 
        	long, monthly_insol[0],
            monthly_insol[1],monthly_insol[2],monthly_insol[3],monthly_insol[4],monthly_insol[5],
            monthly_insol[6],monthly_insol[7],monthly_insol[8],monthly_insol[9],monthly_insol[10],
            monthly_insol[11], math.fsum(monthly_insol))
        unicode_line = line.encode('utf-8')
        output_file.write(unicode_line)
        
    print "done"
    output_file.close()
    
    





import math


qgis_singlept_func() #Calls the qgis_singlept_func