In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
il y a une librairie : PySolar !
http://pysolar.org/
http://docs.pysolar.org/en/latest/
pip install pysolar
Rq: voir aussi une app js en ligne: http://suncalc.net/
In [2]:
map_coords = (45.1973288, 5.7103223) #( 45.166672, 5.71667 )
In [3]:
import pysolar.solar as solar
import datetime as dt
In [16]:
d = dt.datetime.now()
#d = dt.datetime(2017, 6, 20, 13, 30, 0, 130320)
In [18]:
solar.get_altitude( *map_coords, d)
Out[18]:
In [19]:
solar.get_azimuth(*map_coords, d)
Out[19]:
In [20]:
Alt = [ solar.get_altitude(*map_coords, dt.datetime(2017, 12, 21, h, 0, 0, 0)) for h in range(0, 24) ]
Az = [ solar.get_azimuth(*map_coords, dt.datetime(2017, 12, 21, h, 0, 0, 0)) for h in range(0, 24) ]
In [13]:
Az = np.array( Az )
Az[ Az < -180 ] = Az[ Az < -180 ]+360
In [14]:
plt.plot( Alt )
plt.plot( Az )
plt.plot([0, 24], [0, 0], ':'); plt.ylim([-120, 120]); plt.xlabel('hour of the day');
In [21]:
import pysolar.radiation as radiation
In [23]:
radiation.get_radiation_direct( d, 15 ) # W/m2
Out[23]:
Remarque: Le flux solaire au dessus de l'atmosphère est de F = 1 360,8 W/m2
In [12]:
from numpy import genfromtxt
In [13]:
horizon_data = genfromtxt('horizon.csv', delimiter=',').T
In [14]:
horizon_data
Out[14]:
In [53]:
def isUpperHorizon( azimuth, altitude_deg ):
h = np.interp(-azimuth, horizon_data[0, :], horizon_data[1, :])
if h > altitude_deg:
return 0
else:
return 1
In [54]:
isUpperHorizon( 20, 2 )
Out[54]:
In [55]:
horizon_data[1, :].max()
Out[55]:
http://www.a-ghadimi.com/files/Courses/Renewable%20Energy/REN_Book.pdf
page 414
In [56]:
import math
import pysolar.radiation as radiation
import pysolar.solar as solar
import datetime as dt
In [57]:
def get_radiation_direct(d, alt):
if alt>0:
return radiation.get_radiation_direct( d, alt ) # W/m2
else:
return 0
def get_flux_surface( coords, date, sigma, phi_C ):
# Surface orientation :
# sigma : deg, vertical angle of the surface, ref. to the horizontal
# phi_C : deg, azimuth, relative to south, with positive values in the southeast direction and negative values in
# the southwest
# Sun position
phi_S_deg = solar.get_azimuth( *coords, date ) # deg, azimuth of the sun,relative to south
beta_deg = solar.get_altitude( *coords, date ) # deg, altitude angle of the sun
I0 = get_radiation_direct( d, beta_deg ) # W/m2
I0 = I0* isUpperHorizon( phi_S_deg, beta_deg )
beta = beta_deg*math.pi/180 # rad
phi_S = phi_S_deg*math.pi/180 #rad
sigma = sigma*math.pi/180
phi_C = phi_C*math.pi/180
cosTheta = math.cos(beta)*math.cos( phi_S - phi_C )*math.sin( sigma ) + math.cos( sigma )*math.sin( beta )
if cosTheta >0 :
Isurf = I0*cosTheta # flux projeté, W/m2
else:
Isurf = 0 # mais diffuse...
return Isurf
def get_flux_total( coords, date ):
# Sun position
beta_deg = solar.get_altitude( *coords, date ) # deg, altitude angle of the sun
I0 = get_radiation_direct( d, beta_deg ) # W/m2
return I0
In [58]:
get_radiation_direct( d, -4 )
Out[58]:
In [59]:
d = dt.datetime(2017, 6, 22, 11, 0, 0, 0)
sigma = 37
phi_C = 50
F = get_flux_surface( map_coords, d, sigma, phi_C )
print( F )
In [66]:
d
Out[66]:
In [60]:
import pandas as pd
In [61]:
start = dt.datetime(2017, 6, 22, 0, 0, 0, 0)
end = dt.datetime(2017, 6, 22, 23, 59, 0, 0)
d_range = pd.date_range( start=start, end=end, freq='5min' )
In [62]:
F_tot = [ get_flux_total(map_coords, d ) for d in d_range ]
F_est = [ get_flux_surface(map_coords, d, sigma, phi_C ) for d in d_range ]
F_ouest = [ get_flux_surface(map_coords, d, sigma, phi_C+180 ) for d in d_range ]
F_sud = [ get_flux_surface(map_coords, d, 90, phi_C-90 ) for d in d_range ]
In [63]:
x = d_range.hour + d_range.minute/60
plt.figure(figsize=(12, 5))
plt.plot( x, F_est )
plt.plot( x, F_ouest )
plt.plot( x, F_sud )
plt.plot( x, F_tot, 'k:' )
plt.xlabel('hour of the day');
plt.ylabel('flux solaire projeté');
In [64]:
d_range.hour + d_range.minute/60
Out[64]:
In [ ]:
# Sun position
phi_S = solar.get_azimuth( *map_coords, d ) # deg, azimuth of the sun,relative to south
beta = solar.get_altitude( *map_coords, d ) # deg, altitude angle of the sun
I0 = radiation.get_radiation_direct( d, 65 ) # W/m2
cosTheta = math.cos(beta)*math.cos( phi_S - phi_C )*math.sin( sigma ) + math.cos( sigma )*math.sin( beta )
Isurf = I0*cosTheta # flux projeté, W/m2
In [ ]:
cosTheta
In [51]:
Azi = np.array( [ solar.get_azimuth( *map_coords, d ) for d in d_range ] )
Azi[ Azi < -180 ] = Azi[ Azi < -180 ]+360
Alt = [ solar.get_altitude( *map_coords, d ) for d in d_range ]
Hor = [ np.interp(-a, horizon_data[0, :], horizon_data[1, :]) for a in Azi ]
In [65]:
plt.plot( Azi, Hor )
plt.plot( Azi, Alt )
plt.ylim([0, 80]);
In [40]:
Azi
Out[40]:
In [ ]: