%%% this script computes the daily lighting and change in daily lighting
%%% for sites over a year. It's a little more complex than it
%%% probably needs to be so it can handle times of complete darkness/light,
%%% solar eclipses, ...
%%%
%%% author: jens ramrath, agi
%%% date: 3 jan 2020


%% build scenario
% start new instance of STK
uiApp = actxserver('STK11.application');
uiApp.Visible = 1;
root = uiApp.Personality2;

% set up scenario
root.NewScenario('DaylightSingleSite');
sc = root.CurrentScenario;
sc.Animation.AnimStepValue = 86400;
root.Rewind;
root.ExecuteCommand('Terrain * TerrainServer UseTerrainForAnalysis No');

% facilities
fac1 = sc.Children.New('eFacility', 'AGI_HQ');
fac1.Position.AssignGeodetic(40.0386, -75.5966, 0.0);
startTime(1) = {'31 Dec 2019 19:00:00.0'};
stopTime(1) = {'31 Dec 2020 19:00:00.0'};

fac2 = sc.Children.New('eFacility', 'Singapore');
fac2.Position.AssignGeodetic(1.28677, 103.85451, 0.0);
startTime(2) = {'1 Jan 2020 08:00:00.0'};
stopTime(2) = {'1 Jan 2021 08:00:00.0'};

fac3 = sc.Children.New('eFacility', 'Svalbard');
fac3.Position.AssignGeodetic(78.22831, 15.60047, 0.0);
startTime(3) = {'31 Dec 2019 23:00:00.0'};
stopTime(3) = {'31 Dec 2020 23:00:00.0'};

%% loop through facilities
facs = sc.Children.GetElements('eFacility');


% loop through facilities
for f = 0:facs.Count - 1
    root.UnitPreferences.SetCurrentUnit('DateFormat', 'UTCG');
    sc.SetTimePeriod(cell2mat(startTime(f + 1)), cell2mat(stopTime(f + 1)));
    root.UnitPreferences.SetCurrentUnit('DateFormat', 'EpDay');
    
    t1 = datetime(2020,1,1) + caldays(0:str2num(sc.AnalysisInterval.FindStopTime) - 1);
    t2 = datetime(2020,1,1) + caldays(0:str2num(sc.AnalysisInterval.FindStopTime) - 2);
    
    % loop through days
    for day = 0:str2num(sc.AnalysisInterval.FindStopTime) - 1
        fac = facs.Item(cast(f, 'int32'));
        dp = fac.DataProviders.GetDataPrvIntervalFromPath('Lighting Times/Sunlight');
        results = dp.Exec(day, day + 1);
        % if eclipse, there may be multiple lighting intervals per day
        % if no sunlight at all, no data is returned
        if results.DataSets.Count > 0
            lightingDuration(day + 1) = sum(cell2mat(results.DataSets.GetDataSetByName('Duration').GetValues));
        else
            lightingDuration(day + 1) = 0;
        end
        
        if day > 0
            lightingChange(day) = lightingDuration(day+1) - lightingDuration(day);
        end


    end
    
    ff = figure(f+1);
    
    plot(t1, lightingDuration/60.0/60.0, 'k');
    title(['Daylight Durations - ' strrep(fac.InstanceName, '_', ' ') ' (latitude ' num2str(fac.Position.QueryPlanetodetic) ')']);
    ylabel('Daily Sunlight (hours)');
	if f == 1
        ylim([11.8, 12.2]);
    end
    yyaxis right
    plot(t2, lightingChange);
    ylabel('Daily Sunlight Change (sec)');
    xticks('manual')
    xticks(datetime(2020,1,1) + calmonths(0:11));
    if f == 1
        ylim([-6, 6]);
    end
    
    set(gca,'LooseInset',get(gca,'TightInset'));
    set(gcf, 'OuterPosition', [100, 100, 800, 400]);
end


%% svalbard Sun elevation
sunDP = fac3.DataProviders.GetDataPrvTimeVarFromPath('Lighting AER');

% last day of no Sun - 1
sunResults = sunDP.Exec(47.0, 48.0, 60.0);
day47 = cell2mat(sunResults.DataSets.GetDataSetByName('Elevation').GetValues);

% last day of no Sun
sunResults = sunDP.Exec(48.0, 49.0, 60.0);
day48 = cell2mat(sunResults.DataSets.GetDataSetByName('Elevation').GetValues);

% first day of Sun
sunResults = sunDP.Exec(49.0, 50.0, 60.0);
day49 = cell2mat(sunResults.DataSets.GetDataSetByName('Elevation').GetValues);

% next day
sunResults = sunDP.Exec(50.0, 51.0, 60.0);
day50 = cell2mat(sunResults.DataSets.GetDataSetByName('Elevation').GetValues);

figure(4)
time = hours(0:1/60:24);
plot(time, day47, time, day48, time, day49, time, day50);
patch([1 25 25 1],[-0.275 -0.275 -10 -10],[0.4660 0.6740 0.1880], 'EdgeColor', 'none');
patch([1 25 25 1],[0.275 0.275 -0.275 -0.275],[0.9290 0.6940 0.1250], 'EdgeColor', 'none');
patch([1 25 25 1],[0.275 0.275 10 10],'y', 'EdgeColor', 'none');
alpha(0.3);
legend('17 Feb','18 Feb','19 Feb','20 Feb');
title('Svalbard Sun elevation');
ylabel('Sun elevation angle (deg)');
xlim([hours(10.25), hours(14.25)]);
ylim([-1.0, 1.0]);
set(gca,'LooseInset',get(gca,'TightInset'));
set(gcf, 'OuterPosition', [100, 100, 800, 400]);