%% Connect to the existing scenario
% Get reference to running STK instance
uiApplication = actxGetRunningServer('STK12.application');
% Get our IAgStkObjectRoot interface
root = uiApplication.Personality2;

fprintf('Connected to the scenario.\nExtracting objects and preparing to calculate...\n')

%% Set query time
root.UnitPreferences.Item('DateFormat').SetCurrentUnit('UTCG')
sc = root.CurrentScenario;
dt = 600;    % cadence in seconds
%% Define FOM paths and get references to them
% path format: "*/CoverageDefinition/<CoverageName>/FigureOfMerit/<FOM_name>"

% FOM GridPoint2SatRange (gives distance from grid to satellite; 0 if no access)
FOM1Path = '*/CoverageDefinition/CoverageDefinition1/FigureOfMerit/GridPoint2SatRange';
oFom1 = root.GetObjectFromPath(FOM1Path);
oFom1.InstanceName;

% FOM GridPoint2SatAngle
FOM2Path = '*/CoverageDefinition/CoverageDefinition1/FigureOfMerit/GridPoint2SatAngle';
oFom2 = root.GetObjectFromPath(FOM2Path);
oFom2.InstanceName;

% FOM SZA
FOM3Path = '*/CoverageDefinition/CoverageDefinition1/FigureOfMerit/SZA';
oFom3 = root.GetObjectFromPath(FOM3Path);
oFom3.InstanceName;

% Get Coverage grid properties (needed for grid area)
Cov = root.GetObjectFromPath('*/CoverageDefinition/CoverageDefinition1');
elements = {'Point Number';'Latitude';'Longitude';'Area'};
grid_data = Cov.DataProviders.Item('Grid Point Locations').ExecElements(elements);
point_numbers = cell2mat(grid_data.DataSets.GetDataSetByName('Point Number').GetValues());
lats = cell2mat(grid_data.DataSets.GetDataSetByName('Latitude').GetValues());
lons = cell2mat(grid_data.DataSets.GetDataSetByName('Longitude').GetValues());
areas = cell2mat(grid_data.DataSets.GetDataSetByName('Area').GetValues());
GridInfo = [point_numbers lats lons areas]; %gives grid information

% CovDP = Cov.DataProviders.Item('Grid Point Locations');
% CovResult = CovDP.Exec();  % this takes a while
% Covlat = CovResult.DataSets.GetDataSetByName('Latitude').GetValues; % not
% really needed
% Covlon = CovResult.DataSets.GetDataSetByName('Longitude').GetValues;
%ArrayCellArea = cell2mat(CovResult.DataSets.GetDataSetByName('Area').GetValues);   %km^2

%% Time loop (extract different FOMs and do computation)
root.UnitPreferences.Item('DateFormat').SetCurrentUnit('EpSec') % easy to work in EpSec
tstart = sc.StartTime;
tend = sc.StopTime;
tarray = tstart:dt:tend;
Earthshine_intensity = zeros(size(tarray)); % init
i = 1;

for t = tstart:dt:tend

    t_current = root.ConversionUtility.ConvertDate('EpSec','UTCG',num2str(t)); % back to UTCG
    fprintf('Processing for time = %s \n',t_current)
    dataPrv1 = oFom1.DataProviders.Item('Time Value By Point');
    dataPrv1.PreData =  t_current;

    dataPrv2 = oFom2.DataProviders.Item('Time Value By Point');
    dataPrv2.PreData = t_current;

    dataPrv3 = oFom3.DataProviders.Item('Time Value By Point');
    dataPrv3.PreData = t_current;

    % Extract FOM data for the timestep 

    RangeResult = ExecElements(dataPrv1,{'Latitude'; 'Longitude'; 'FOM Value'});
    %Pull out the individual data sets
    Lats = RangeResult.DataSets(1).GetDataSetByName('Latitude');
    Lons = RangeResult.DataSets(1).GetDataSetByName('Longitude');
    CovRange = RangeResult.DataSets(1).GetDataSetByName('FOM Value');
    %Get the values of the datasets
    ArrayCovRange = [cell2mat(Lats.GetValues()), cell2mat(Lons.GetValues()), cell2mat(CovRange.GetValues())];

    GridPt2SatResult = ExecElements(dataPrv2,{'Latitude'; 'Longitude'; 'FOM Value'});
    %Pull out the individual data sets
    Lats = GridPt2SatResult.DataSets(1).GetDataSetByName('Latitude');
    Lons = GridPt2SatResult.DataSets(1).GetDataSetByName('Longitude');
    GridPt2SatAngle = GridPt2SatResult.DataSets(1).GetDataSetByName('FOM Value');
    %Get the values of the datasets
    ArrayGridPt2SatAngle = [cell2mat(Lats.GetValues()), cell2mat(Lons.GetValues()), cell2mat(GridPt2SatAngle.GetValues())];

    SZAResult = ExecElements(dataPrv3,{'Latitude'; 'Longitude'; 'FOM Value'});
    %Pull out the individual data sets
    Lats = SZAResult.DataSets(1).GetDataSetByName('Latitude');
    Lons = SZAResult.DataSets(1).GetDataSetByName('Longitude');
    SZA = SZAResult.DataSets(1).GetDataSetByName('FOM Value');
    %Get the values of the datasets
    ArraySZA = [cell2mat(Lats.GetValues()), cell2mat(Lons.GetValues()), cell2mat(SZA.GetValues())];

    % computation of Earthshine intensity (W/m^2)
    Earthshine_intensity(i) = ComputeEarthshine(1367,0.3,ArrayCovRange(:,3),ArraySZA(:,3),ArrayGridPt2SatAngle(:,3),GridInfo(:,4));
    i = i+1;
end

% plot the computed result
figure()
plot(tarray,Earthshine_intensity,'LineWidth',1)
xlabel('EpSec')
ylabel('Earthshine (W/m^2)')


function[Earthshine] = ComputeEarthshine(SolarConstant,albedo,Range,SZA,Grid2SatAngle,GridArea)
% inputs: SolarConstant (1367 W/m^2 nominally) and albedo (0.3 nominally) constants
%         Range: distance from Gridpoint to target in km (nx1 vector)
%         SZA: solar zenith angle in deg (nx1 vector)
%         Grid2SatAngle: angle between grid normal to target in deg (nx1 vector)
%         GridArea: Area of grid cells in km^2 (nx1 vector)
% output: Earthshine: computed Earthshine intensity at current timestep in W/m^2


    EarthshinePerCell = zeros(size(Range));
    SZA(abs(SZA-999) < 1e-6) = NaN; % writing invalid (999.0) values as NaN
    Grid2SatAngle(abs(Grid2SatAngle-999) < 1e-6) = NaN;
    
    SZA(SZA > 90) = 90; % binding SZA values to +- 90 deg
    SZA(SZA < -90) = -90;
            
    EarthshinePerCell = SolarConstant.*cosd(SZA).*cosd(Grid2SatAngle).*albedo.*double(GridArea)./(Range.^2);   %W/m2
    EarthshinePerCell(isnan(EarthshinePerCell)) = 0; % set all nan values to zero
    EarthshinePerCell(EarthshinePerCell < 0) = 0;   % set negative values as zero
    Earthshine = sum(EarthshinePerCell); % Earthshine at this time instant (nansum not needed since NaN values already removed)
end