-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmonthly.m
146 lines (119 loc) · 5.41 KB
/
monthly.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
% Constants
solar_constant = 1353; % W/m²
Cn = 0.75; % PV system efficiency
latitude = 42.984333; % Syracuse, NY
longitude = -76.142167;
tilt_angle = 30; % Assumed tilt angle (degrees)
surface_azimuth = 10; % Assumed collector azimuth angle (degrees)
albedo = 0.2; % Ground reflectance
household_consumption = 10000; % Yearly household electricity consumption in kWh
% Dates for simulation (one day per month)
dates = ["2024-01-15", "2024-02-15", "2024-03-15", "2024-04-15", "2024-05-15", "2024-06-15", "2024-07-15", "2024-08-15", "2024-09-15", "2024-10-15", "2024-11-15", "2024-12-15"];
labels = ["January 15", "February 15", "March 15", "April 15", "May 15", "June 15", "July 15", "August 15", "September 15", "October 15", "November 15", "December 15"];
% Times from 4 AM to 9 PM (in 1-hour increments)
times = arrayfun(@(x) sprintf('%02d:00', x), 4:21, 'UniformOutput', false);
% Prepare figure
figure;
hold on;
colors = lines(length(dates));
% Variables to store monthly energy and total yearly energy
monthly_energy = zeros(1, length(dates));
total_yearly_energy = 0;
% Variables to store sunrise and sunset times
sunrise_times = strings(1, length(dates));
sunset_times = strings(1, length(dates));
for i = 1:length(dates)
date = dates{i};
label = labels{i};
radiation_values = zeros(1, length(times));
non_zero_times = {};
for j = 1:length(times)
time = times{j};
Ic = solar_radiation(latitude, longitude, date, time, tilt_angle, surface_azimuth, albedo, Cn, 0.144, 0.06, solar_constant);
radiation_values(j) = Ic;
if Ic > 0
non_zero_times{end+1} = time;
end
end
% Calculate daily energy in kWh/m²
daily_energy = sum(radiation_values) / 1000; % Convert Wh/m² to kWh/m²
monthly_energy(i) = daily_energy * 30; % Approximate monthly energy (30 days)
% Plot the radiation values
plot(1:length(times), radiation_values, 'DisplayName', label, 'Color', colors(i,:));
% Calculate and store sunrise and sunset times
if ~isempty(non_zero_times)
first_non_zero = datetime(non_zero_times{1}, 'InputFormat', 'HH:mm');
last_non_zero = datetime(non_zero_times{end}, 'InputFormat', 'HH:mm');
sunrise_times(i) = datestr(first_non_zero, 'HH:MM AM');
sunset_times(i) = datestr(last_non_zero, 'HH:MM PM');
else
sunrise_times(i) = 'N/A';
sunset_times(i) = 'N/A';
end
end
% Sum up total annual energy
total_yearly_energy = sum(monthly_energy);
% Customize the plot
xlabel('Time (Hour)');
ylabel('Solar Radiation (W/m²)');
title('Total Solar Radiation on the Panel vs. Local Time');
legend show;
grid on;
hold off;
% Display sunrise and sunset times
for i = 1:length(dates)
fprintf('%s:\n', labels{i});
fprintf(' Sunrise: %s\n', sunrise_times(i));
fprintf(' Sunset: %s\n', sunset_times(i));
end
% Display monthly energy and total annual energy
fprintf('\nMonthly Energy (kWh/m²):\n');
for i = 1:length(dates)
fprintf('%s: %.2f kWh/m²\n', labels{i}, monthly_energy(i));
end
fprintf('\nTotal Annual Energy: %.2f kWh/m²\n', total_yearly_energy);
% Define the function to calculate solar radiation
function Ic = solar_radiation(lat, lon, date, time, tilt_angle, surface_azimuth, albedo, Cn, k, C, solar_constant)
% Convert inputs to radians
lat = deg2rad(lat);
tilt_angle = deg2rad(tilt_angle);
surface_azimuth = deg2rad(surface_azimuth);
% Parse the date and time
datetime_obj = datetime(date + " " + time, 'InputFormat', 'yyyy-MM-dd HH:mm');
% Calculate day of the year (n)
day_of_year = day(datetime_obj, 'dayofyear');
% Equation of Time (ET)
B = deg2rad(360 * (day_of_year - 81) / 364); % Radians
ET = 9.87 * sin(2 * B) - 7.53 * cos(B) - 1.5 * sin(B); % Minutes
% Convert EST to Solar Time
local_standard_meridian = -75; % For EST (degrees)
solar_time_offset = ET + 4 * (lon - local_standard_meridian);
solar_time = datetime_obj + minutes(solar_time_offset);
solar_hour_angle = deg2rad(15 * (hour(solar_time) + minute(solar_time) / 60 - 12)); % Radians
% Solar Declination (delta)
declination = deg2rad(23.45 * sin(deg2rad(360 / 365 * (day_of_year - 81))));
% Solar Altitude Angle (alpha)
sin_alpha = sin(lat) * sin(declination) + cos(lat) * cos(declination) * cos(solar_hour_angle);
alpha = asin(sin_alpha);
% Direct Normal Irradiance (Ib,n)
if rad2deg(alpha) > 0
Ib_n = Cn * solar_constant * (1 + 0.034 * cos(deg2rad(360 * day_of_year / 365.25))) * exp(-k / sin(alpha));
else
Ib_n = 0;
end
% Diffuse Radiation (Id)
Id = C * Ib_n * cos(tilt_angle / 2) * cos(tilt_angle / 2);
% Ground Reflected Radiation (Ir)
Ir = albedo * Ib_n * sin(alpha + C) * sin(tilt_angle / 2) * sin(tilt_angle / 2);
% Angle of Incidence (theta)
cos_theta = (sin(declination) * sin(lat) * cos(tilt_angle)) - ...
(sin(declination) * cos(lat) * sin(tilt_angle) * cos(surface_azimuth)) + ...
(cos(declination) * cos(lat) * cos(tilt_angle) * cos(solar_hour_angle)) + ...
(cos(declination) * sin(lat) * sin(tilt_angle) * cos(surface_azimuth) * cos(solar_hour_angle)) + ...
(cos(declination) * sin(tilt_angle) * sin(surface_azimuth) * sin(solar_hour_angle));
cos_theta = max(cos_theta, 0);
% Beam Radiation on the Collector (Ib_n * cos(theta))
Ib_cos_theta = Ib_n * cos_theta;
% Total Radiation on the Collector (Ic)
Ic = Ib_cos_theta + Id + Ir;
end