-
Notifications
You must be signed in to change notification settings - Fork 0
/
test3.m
131 lines (104 loc) · 4.95 KB
/
test3.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
% Solar Radiation Calculation and Plotting
% Created on Mon Nov 18 23:21:19 2024
% Author: Kalhaku Mclester
% Constants
solar_constant = 1353; % W/m²
% Parameters
latitude = 42.984333; % Syracuse, NY
longitude = -76.142167;
tilt_angle = 30; % Collector tilt angle (degrees)
surface_azimuth = 10; % Collector azimuth angle (degrees)
albedo = 0.2; % Ground reflectance
% Dates for different seasons
dates = ["2024-03-21", "2024-06-21", "2024-09-21", "2024-12-21"];
labels = ["March 21", "June 21", "September 21", "December 21"];
% Times from 4 AM to 9 PM (in 1-minute increments)
times = arrayfun(@(x) sprintf('%02d:%02d', floor(x/60), mod(x, 60)), 4*60:21*60-1, 'UniformOutput', false);
% Sunrise and Sunset Errors (in HH:MM format)
actual_times = containers.Map({'March 21', 'June 21', 'September 21', 'December 21'}, ...
{{'7:05 AM', '7:18 PM'}, {'5:25 AM', '8:47 PM'}, {'6:51 AM', '7:03 PM'}, {'7:33 AM', '4:33 PM'}});
% Plotting
figure('Position', [100, 100, 1200, 700]);
hold on;
% Calculate and plot solar radiation
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, 1, 0.144, 0.06);
radiation_values(j) = Ic;
if Ic > 0
non_zero_times{end+1} = time;
end
end
% Plot solar radiation
plot(1:length(times), radiation_values, 'DisplayName', label);
% Calculate sunrise and sunset
first_non_zero = datetime(non_zero_times{1}, 'InputFormat', 'HH:mm') + hours(1);
last_non_zero = datetime(non_zero_times{end}, 'InputFormat', 'HH:mm') + hours(1);
predicted_sunrise = datestr(first_non_zero, 'HH:MM PM');
predicted_sunset = datestr(last_non_zero, 'HH:MM PM');
% Print the times and errors
actual_sunrise = datetime(actual_times(label){1}, 'InputFormat', 'hh:mm a');
actual_sunset = datetime(actual_times(label){2}, 'InputFormat', 'hh:mm a');
sunrise_error = minutes(first_non_zero - actual_sunrise);
sunset_error = minutes(last_non_zero - actual_sunset);
fprintf('%s:\n', label);
fprintf(' Predicted Sunrise: %s, Actual Sunrise: %s, Error: %.2f minutes\n', predicted_sunrise, actual_times(label){1}, sunrise_error);
fprintf(' Predicted Sunset: %s, Actual Sunset: %s, Error: %.2f minutes\n', predicted_sunset, actual_times(label){2}, sunset_error);
end
% Plot settings
xlabel('Time');
ylabel('Solar Radiation (W/m²)');
title('Solar Radiation throughout the Day for Different Seasons');
legend('show');
grid on;
hold off;
% Function to calculate solar radiation
function Ic = solar_radiation(lat, lon, date, time, tilt_angle, surface_azimuth, albedo, Cn, k, C)
% 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