clear; close all; clc;


%%
Latitude = 38; place = 'Athen'; luxmax = 120;

%%
dateaxis = 0;
for ii = 1:12
    dateaxis = [dateaxis,datenum(2000,ii,1),datenum(2000,ii,15)];
end
dateaxis = [dateaxis(2:end),datenum(2000,12,31)];
dateaxis = dateaxis - dateaxis(1) + 1;

timeaxis = linspace(datenum(2000,1,1,2,0,0),datenum(2000,1,1,22,0,0),11);

%%
hour = linspace(0,24,251).';
day  = 1:366;
Elevation = zeros(length(hour),length(day));

for ii = 1:length(day)
    UTC = datenum([2001,1,1,0,0,0])-1+day(ii);
    UTC = UTC + hour/24;
    Lat = Latitude * ones(size(UTC));
    Lon = zeros(size(UTC));
    Alt = zeros(size(UTC));

    [Az,El] = SolarAzEl(UTC,Lat,Lon,Alt);
    
    Elevation(:,ii) = El;
end

figure();

imagesc(day,hour,Elevation)
hold on;
contour(day,hour,Elevation,[0,0],'-k');

ylabel('time of day');
set(gca,'XTick',dateaxis);
set(gca,'XTickLabel',datestr(dateaxis,'dd.mmm'));
xtickangle(90)
colorbar;
title(sprintf('Solar Altitude at %d ° latitude (e.g. %s)',Latitude,place))


%%
Elevatation2 = Elevation;
Elevatation2(Elevatation2<0)=0;
Lux = luxmax/exp(-0.5) * exp(-0.5 * csc( Elevatation2/180*pi ) );

f = figure('Position',[50,50,600,500]); 
axes('Position',[0.1,0.5,0.78,0.4]); 
imagesc(day,hour,Lux)
hold on;
contour(day,hour,Elevation,[0,0],'-w');

set(gca,'XTick',dateaxis);
set(gca,'XTickLabel',datestr(dateaxis,'dd.mmm'));
xtickangle(90)

set(gca,'YTick',2:2:24);
set(gca,'YTickLabel',datestr(timeaxis,'HH:MM'));
h_cb  = colorbar;
set(get(h_cb,'ylabel'),'String','Beleuchtungsstärke (klx)');
set(h_cb,'Position',[0.9,0.5,0.02,0.4])
title(sprintf('Beleuchtungsstärke bei %d° Breitengrad (zB %s)',Latitude,place))
ylim([4,20])

axes('Position',[0.1,0.1,0.78,0.3]); 
plot(day,Lux((end+1)/2,:),'k-')
xlim([1,366]);
set(gca,'XTick',[]);
ylim([0,1.1*luxmax])


axes('Position',[0.1,0.1,0.78,0.3],'Color','none');
[~,id1] = min(Lux((end+1)/2,:));
[~,id2] = max(Lux((end+1)/2,:));
plot(hour,Lux(:,id2),'k-');
hold on;
plot(hour,Lux(:,id2),'r-')
plot(hour,Lux(:,id1),'r--');
xlim([4,22]);
ylim([0,1.1*luxmax])
set(gca,'Color','none');
set(gca,'XTick',2:2:24);
set(gca,'XTickLabel',datestr(timeaxis,'HH:MM'));
legend('Maximum übers Jahr','Tagesverlauf Sommer','Tagesverlauf Winter')
ylabel('Beleuchtungsstärke (klx)')

print(f,'lux_year.png','-dpng');