Another approach —
n = 32;
S = 23.45*sind(360*(284 + n)/365);
phi = 38.9072;
tsol = 0:1/60:24;
w = 15*(tsol-12);
Gsc = 1353;
Go = Gsc*(1 + 0.033*cosd((360*n)/365))*(cosd(phi)*cosd(S)*cosd(w) + sind(phi)*sind(S));
plot(tsol,Go)
hold on
Y = 30;
B = 50;
theta1 = sind(S)*sind(phi)*cosd(B) - sind(S)*cosd(phi)*sind(B)*cosd(Y) + cosd(S)*cosd(phi)*cosd(B)*cosd(w)+ cosd(S)*sind(phi)*sind(B)*cosd(Y)*cosd(w) + cosd(S)*sind(B)*sind(Y)*sind(w);
Got = Gsc*(1 + 0.033*cosd((360*n)/365))*theta1;
plot(tsol,Got)
Go_idx = find(diff(sign(Go))); % Approximate Indices
for k = 1:numel(Go_idx)
idxrng = max([1,Go_idx(k)-1]) : min([numel(Go),Go_idx(k)+1]); % Index Range For Interpolation
Go_t(k) = interp1(Go(idxrng), tsol(idxrng), 0); % Precise Values
end
Got_idx = find(diff(sign(Got))); % Approximate Indices
for k = 1:numel(Got_idx)
idxrng = max([1,Got_idx(k)-1]) : min([numel(Got),Got_idx(k)+1]); % Index Range For Interpolation
Got_t(k) = interp1(Got(idxrng), tsol(idxrng), 0); % Precise Values
end
plot([Go_t Got_t], zeros(1,numel(Go_t)+numel(Got_t)), '+k','MarkerSize',10)
hold off
grid on
title('Homework #3, Question 1b')
xlabel('Time of day (Hours)')
ylabel('Solar Radiation on horizontal plane')
xlim([0 24])
legend('Horrizontal','Angled')
syms tsol Go
tsol = solve([0 == Gsc*(1 + 0.033*cosd((360*n)/365))*(cosd(phi)*cosd(S)*cosd(w) + sind(phi)*sind(S))-Go], tsol)
.
