Skip to content

GitLab

Commit eb0c3e35 authored by Angioni's avatar Angioni
Browse files

inseretd functions and readme

parent f943dde1
Expand all Hide whitespace changes
Inline Side-by-side
DHI_calculation.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate direct sun irradiation
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DHI] = DHI_calculation(Extraterr_radiation, theta_zenith)
% Tlk = 4.5; % value of Linke's turbidity factor when the sky is fully clean and dry (cit. Comparison Between Atmospheric Turbidity Coefficients of Desert and Temperate Climates)
Tlk = 3;
h0 = (pi/2)-theta_zenith;
Tn = -0.015843 + 0.030543*Tlk + 0.0003797*(Tlk^2); % this and folllowing values are collected from: http://re.jrc.ec.europa.eu/pvgis/solres/solmod3.htm
A1 = 0.26463 - 0.061581*Tlk + 0.0031408*(Tlk^2);
A2 = 2.04020 + 0.018945*Tlk - 0.011161*(Tlk^2);
A3 = -1.3025 + 0.039231*Tlk + 0.0085079*(Tlk^2);
Fd = A1 + A2.*sin(h0) + A3.*(sin(h0).^2);
DHI = Extraterr_radiation*Tn*Fd;
\ No newline at end of file
DNI_calculation.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate direct irradiation
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DNI] = DNI_calculation(Extraterr_radiation, theta_zenith)
% h0 = pi/2 - theta_zenith;
% Dh0ref = 0.061359.*(0.1594 + 1.123.*h0 + 0.065656*(h0.^2))./(1 + 28.9344.*h0 + 277.3971*(h0.^2));
% h0ref = h0 + Dh0ref;
% Tlk = 4.5; % value of Linke's turbidity factor when the sky is fully clean and dry (cit. Comparison Between Atmospheric Turbidity Coefficients of Desert and Temperate Climates)
Tlk = 3;
% m = 1/(cos(theta_zenith) + 0.15.*(93.885 - theta_zenith*180/pi).^(-1.253)); % relative optical air mass (cit. Comparison Between Atmospheric Turbidity Coefficients of Desert and Temperate Climates)
% m = 1./(sin(h0ref) + 0.50572*(h0ref + 6.07995).^(-1.6364));
% Air mass calculation: reference of formula -> formula of Young (1994) on: https://en.wikipedia.org/wiki/Air_mass_%28astronomy%29#CITEREFYoung1994
m = (1.002432.*(cos(theta_zenith).^2) + 0.148386.*cos(theta_zenith) + 0.0096467) ./ ((cos(theta_zenith).^3) + 0.149864.*(cos(theta_zenith).^2) + 0.0102963.*cos(theta_zenith) + 0.000303978);
% Computation spectral optical thickness of the clean dry atmosphere dr
% (http://re.jrc.ec.europa.eu/pvgis/solres/solmod3.htm Kasten 1996)
if m <=20
dr = 1./(6.6296 + 1.7513.*m - 0.1202.*(m.^2) + 0.0065*(m.^3) - 0.00013.*(m.^4));
else
dr = 1./(10.4 + 0.718.*m);
end
DNI = Extraterr_radiation.*exp(-0.8662*Tlk.*m.*dr);
Gen_Base_Profile.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to generate base power profile
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DNI,DHI] = Gen_Base_Profile(num_day)
% GHI_matrix = zeros(2401,365);
% POA_matrix = zeros(2401,365);
% lat = 50.77;
% lon = 6.09;
% lat = 39.2167;
% lon = 9.1167;
lat = 36.06;
lon = -115.08;
GMT_diff = -8;
albedo = 0.2;
tilt_pv = 30;
azimuth_pv = 180;
%1 sunny
%2 slightly cloudy
%3 cloudy
%4 total overcast
% cloudiness = [0 10 4;
% 10 12.75, 3;
% 12.75 13.2 2;
% 13.2 24 1];
% cloudiness = [0 8.5 2;
% 8.5 12.0, 2;
% 12.0 13.0 2;
% 13.0 15.5 3;
% 15.5 24 3];
cloudiness = [ 0 24 1];
% for i =1:365
%num_day = 33; %20 79 171
[theta_zenith, ~] = sun_position(num_day, lat, lon, GMT_diff);
[~, DNI, DHI, ~] = irradiation_calc_test(num_day, theta_zenith, cloudiness);
% [PV_output] = POA_irradiance(GHI, DNI, DHI, theta_zenith, theta_azimuth, extraterr_radiation, albedo, tilt_pv, azimuth_pv);
% GHI_matrix(:,i) = GHI;
% POA_matrix(:,i) = POA;
% end
end
Mean_direct.mat 0 → 100644
View file @ eb0c3e35
File added
POA_irradiance.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate PV output
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [PV_output] = POA_irradiance(GHI, DNI, DHI, theta_zenith, theta_azimuth, extraterr_radiation, albedo, tilt_pv, azimuth_pv)
tilt_pv = tilt_pv*pi/180;
azimuth_pv = azimuth_pv*pi/180;
AOI = acos(cos(theta_zenith).*cos(tilt_pv) + sin(theta_zenith)*sin(tilt_pv).*cos(theta_azimuth-azimuth_pv));
Ai = DNI/extraterr_radiation;
Eb = DNI.*cos(AOI);
Ed = DHI.*(Ai.*cos(AOI) + (1-Ai).*(1+cos(tilt_pv))/2);
% Eg = 0;
Eg = GHI*albedo*(1-cos(tilt_pv))/2;
POA = Eb + Eg + Ed;
POA(POA<0) = 0;
if mod(length(POA),2) == 1
POA1 = POA(1:length(POA)/2-0.5);
POA2 = POA(length(POA)/2+0.5:end);
elseif mod(length(POA),2) == 0
POA1 = POA(1:length(POA)/2);
POA2 = POA(length(POA)/2+1:end);
end
idx1 = find(POA1==0);
len_idx1 = length(idx1);
idx2 = idx1(len_idx1,1);
POA1(1:idx2,1) = 0;
idx1 = find(POA2==0);
idx2 = idx1(1,1);
POA2(idx2:end,1) = 0;
PV_output = [POA1; POA2]/1000;
% PV_output = [POA1; POA2];
PV_generator.m 0 → 100644
View file @ eb0c3e35
This diff is collapsed.
README.md 0 → 100644
View file @ eb0c3e35
# Software tool to generate PV power profiles
This repository contains the following material:
- Matlab code to generate realistic PV power injection profiles for 1 day.
## Instructions
1) Run in the command window "PVgenerator.m"
2) Define, in the GUI, the following parameters:
- Location inputs. Latitude, longitude
- Time inputs Time of the year
- Cloudiness settings: Cloudiness level in different parts of the day
- PV settings: Geometrical placement and power settings of PV panels
3) Press "Generate profiles".
4) This will generate a figure with the profile in the GUi and create a file "PV_Power_profiles.mat" with all the power profiles, with time resolution of 1 second, in the current directory.
## Copyright
2019, Institute for Automation of Complex Power Systems, EONERC
## License
This project is released under the terms of the [GPL version 3](COPYING.md).
```
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
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 <http://www.gnu.org/licenses/>.
```
For other licensing options please consult [Prof. Antonello Monti](mailto:amonti@eonerc.rwth-aachen.de).
## Contact
[![EONERC ACS Logo](https://git.rwth-aachen.de/acs/public/villas/VILLASnode/raw/develop/doc/pictures/eonerc_logo.png)](http://www.acs.eonerc.rwth-aachen.de)
- Marco Pau <mpau@eonerc.rwth-aachen.de>
- Andrea Angioni <aangioni@eonerc.rwth-aachen.de>
[Institute for Automation of Complex Power Systems (ACS)](http://www.acs.eonerc.rwth-aachen.de)
[EON Energy Research Center (EONERC)](http://www.eonerc.rwth-aachen.de)
[RWTH University Aachen, Germany](http://www.rwth-aachen.de)
cloud_generator.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate impact of clouds
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DNI, DHI] = cloud_generator(DNI, DHI, cloudiness)
DNI_start = DNI;
DHI_start = DHI;
load('Mean_diffuse.mat');
load('Mean_direct.mat');
mean_dir_tot(1,1) = 0;
mean_dif_tot(1,1) = 0;
for i = 1:size(cloudiness,1)
start_time = cloudiness(i,1)*3600;
end_time = cloudiness(i,2)*3600;
[d_dir, d_dif] = reduction_factor_creation_new(start_time, end_time, cloudiness(i,3));
DNI(start_time+1:end_time,1) = (d_dir' + mean_dir_tot(cloudiness(i,3),1)).*DNI_start(start_time+1:end_time,1) + DNI_start(start_time+1:end_time,1);
DHI(start_time+1:end_time,1) = (d_dif' + mean_dif_tot(cloudiness(i,3),1)).*DNI_start(start_time+1:end_time,1) + DHI_start(start_time+1:end_time,1);
for x = start_time+1 : end_time
if DNI(x) > DNI_start(x)
delta_dir = DNI(x) - DNI_start(x);
DNI(x) = DNI_start(x) - delta_dir;
end
if DNI(x) < 0
DNI(x) = 0;
end
if DHI(x) < 0
DHI(x) = - DHI(x);
end
end
var =300;
if start_time>var
for j = 1:var
DNI(start_time+j,1) = ((var-j)/var)*DNI(start_time,1) + (j/var)*DNI(start_time+j,1);
DHI(start_time+j,1) = ((var-j)/var)*DHI(start_time,1) + (j/var)*DHI(start_time+j,1);
end
end
end
% DHI_start(DHI_start<0) = 0;
% DHI(DHI<0) = 0;
% DNI_start(DNI_start<0) = 0;
% DNI(DNI<0) = 0;
% figure(1)
% hold on
% plot(DNI,'b')
% plot(DHI,'g')
% axis([20000 70000 0 900])
% a = DHI;
irradiation_calc.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate total irradiation
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [GHI, DNI, DHI, Extraterr_radiation] = irradiation_calc(num_day, theta_zenith, cloudiness)
% Air mass calculation: reference of formula -> formula of Young (1994) on: https://en.wikipedia.org/wiki/Air_mass_%28astronomy%29#CITEREFYoung1994
% air_mass = (1.002432.*(cos(theta_zenith).^2) + 0.148386.*cos(theta_zenith) + 0.0096467) ./ ((cos(theta_zenith).^3) + 0.149864.*(cos(theta_zenith).^2) + 0.0102963.*cos(theta_zenith) + 0.000303978);
% air_mass(theta_zenith>pi/2) = -1;
Esolar_constant = 1367; % Energy solar constant in Watt per squared meter
b = 2*pi*num_day/365;
R = 1.00011 + 0.034221*cos(b) + 0.00128*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b); % Multiplicative factor taken from: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/irradiance-and-insolation-2/extraterrestrial-radiation/
Extraterr_radiation = Esolar_constant*R; % in Watts per squared meters, taken from: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/irradiance-and-insolation-2/extraterrestrial-radiation/
[DNI] = DNI_calculation(Extraterr_radiation, theta_zenith); % Direct Normal Irradiance (normal to the sun beam)
[DHI] = DHI_calculation(Extraterr_radiation, theta_zenith); % Diffuse Horizontal Irradiance (on a horizontal surface)
DNI(DNI<0.01) = 0;
idx = find(DNI(1:43200,1)==0);
len = length(idx);
idx_start = idx(len);
DNI(1:idx_start,1) = zeros(idx_start,1);
DHI(1:idx_start,1) = zeros(idx_start,1);
idx = find(DNI(43201:end,1)==0);
idx_end = idx(1)+43200;
len = 86400-idx_end;
DNI(idx_end+1:end,1) = zeros(len,1);
DHI(idx_end+1:end,1) = zeros(len,1);
[DNI, DHI] = cloud_generator(DNI, DHI, cloudiness);
GHI = DHI + DNI.*cos(theta_zenith);
% DNIh = DNI.*cos(theta_zenith);
reduction_factor_creation_new.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate reduction factor
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [d_dir, d_dif] = reduction_factor_creation_new(start_time, end_time, cloudiness)
len = end_time - start_time;
if cloudiness == 1
d_dir = zeros(1,len);
d_dif = zeros(1,len);
else
load('Transition_matrix.mat');
load('White_covariance_matrix.mat');
F_matrix = reshape(F_matrix_out_tot(cloudiness,:,:),[2 2]);
cov_w1 = reshape(cov_w1_out_tot(cloudiness,:,:),[2 2]);
cov_w1 = (cov_w1 + cov_w1.')/2;
d(1:2,1) = mvnrnd([0;0],cov_w1);
% F_matrix
% cov_w1
for x = 2 : len
d(1:2,x) = F_matrix * d(1:2,x-1) + mvnrnd([0;0],cov_w1)';
end
d_dir = d(1,:);
d_dif = d(2,:);
end
sun_position.m 0 → 100644
View file @ eb0c3e35
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate sun position
%
% @author Marco Pau <mpau@eonerc.rwth-aachen.de>
% @copyright 2019, Institute for Automation of Complex Power Systems, EONERC
% @license GNU General Public License (version 3)
%
% PVgenerator
%
% 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
% 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 <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [theta_zenith, theta_azimuth] = sun_position(num_day, lat, lon, GMT_diff)
% Inputs: num_day = number of the day of the year
% lat = latitude of the considered location in decimal degrees (northern emisphere has positive values)
% lon = longitude of the considered location in decimal degrees (east with repect to Greenwich has positive values)
% GMT_diff = Greenwich Mean Time, time offset with respect to UTC
% Outputs: theta_zenith = zenith angle of the sun (in degrees)
% theta_azimuth = azimuth angle of the sun (in degrees)
% air mass = air mass value
% Tlocal=[0:0.01:24]'; % time of the day with resolution of 1 second
Tlocal=[24/86400:24/86400:24]';
theta_d= (23.45*pi/180)*sin(2*pi*(284+num_day)/365); % approximated calculation of the declination angle of the sun (in radians) according to the number of the day (PV Sandia model at: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/sun-position/simple-models/)
if(num_day >= 1 && num_day <= 106) % approximated calculation of the equation of time according to the number of the day (PV Sandia model at: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/sun-position/simple-models/)
eq_time = -14.2*sin(pi*(num_day+7)/111);
elseif(num_day > 106 && num_day <= 166)
eq_time = 4.0*sin(pi*(num_day-106)/59);
elseif(num_day > 166 && num_day <= 246)
eq_time = -6.5*sin(pi*(num_day-166)/80);
elseif(num_day > 246 && num_day <= 366) %andrea --> modified from 365 to 366
eq_time = 16.4*sin(pi*(num_day-247)/113);
end
lon_sm = 15*(GMT_diff); % longitude of the standard meridian of the observer’s time zone (degrees)
Tsolar = Tlocal+eq_time/60+(lon-lon_sm)/15; % Calculation of the solar time based on the position of the sun and not on the local time (wrong sign in Sandia website)
theta_hr = pi*(Tsolar-12)/12; % Calculation of the hour angle (in radians) (wrong sign in Sandia website)
theta_zenith = acos(sin(lat*pi/180)*sin(theta_d)+cos(lat*pi/180)*cos(theta_d).*cos(theta_hr)); % Solar zenith angle (in radians) accordin to PV Sandia model (https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/sun-position/simple-models/)
% Possible validation of zenith angles at: https://midcdmz.nrel.gov/solpos/spa.html
theta_azimuth1 = acos((sin(theta_d)*cos(lat*pi/180)-cos(theta_d)*sin(lat*pi/180).*cos(theta_hr))./sin(theta_zenith)); % Calculation of the solar azimuth angle (Wikipedia)
theta_azimuth2 = asin(-sin(theta_hr).*cos(theta_d)./sin(theta_zenith));
theta_azimuth = theta_azimuth1; % just for initialization purposes
% If procedure to put the azimuth between 0 and pi when the hour angle is negative and between -pi and 0 when the hour angle is positive
for i = 1:length(theta_azimuth)
if theta_azimuth2(i,1) < 0 && theta_azimuth1(i,1) < pi/2
theta_azimuth(i,1) = 2*pi + theta_azimuth2(i,1);
elseif theta_azimuth2(i,1) < 0 && theta_azimuth1(i,1) > pi/2
theta_azimuth(i,1) = 2*pi - theta_azimuth1(i,1);
end
end
% theta_zenith=theta_zenith *180/pi;
% theta_azimuth=theta_azimuth *180/pi;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment