%% Use this script to perform all bacis computations 
% This script contains all simple computation mentioned in section
% 'Sensitivity analyis'. Use this script prior creating figure showing the
% effect of intergration radius and different layers.
% 
addpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
addpath(fullfile('..','MatlabLibrary','aux_fce'));
clear
close all
clc

%% Main settings valid for all computations
% Site names. Will be used for looping the computation
sites = {'sg030','igrav'}; % 'sg029',

% Set GausKrueger coordinates + height interpolated using local DEM.
coor.sg029 = [4564213.48 5445521.70 613.287]; % x, y, Height at the floor! (m)
coor.sg030 = [4564041.87 5445662.88 606.471]; % x, y, Height at the floor! (m)
coor.igrav = [4564082.00 5445669.70 609.755]; % centre of pillar, Height at surface (m)

% Set local coordinate system, i.e., x and y along building
% and z above the ground
coor_relative.sg029 = [3.78 2.44 0.0]; % x, y, z includes base plate!!
coor_relative.sg030 = [8.65 6.32 0.0]; % x, y, z
coor_relative.igrav = [0 0 0]; % centre of pillar

% Set sensor heights relative to floor/surface
sensor.sg029 = 0.25; % m
sensor.sg030 = 0.25; % m
sensor.igrav = 1.05; % m
% Set sensor height for sensitivity test
sensor.igrav_test = 0.25:0.1:2; % m

% Set building dimensions
build.sg029 = [5.99 10.04, 0, 2.5, 30]; % dx, dy (m), z, height to roof, roof angle in degrees
build.sg030 = [10.99 7.99, 0, 2.5, 30]; 
pillar.sg029 = [1.1 1.1 1.7]; % dx, dy, depth below surface (m)
pillar.sg030 = [1.1 1.1 1.7]; % dx, dy, depth below surface (m)
pillar.igrav = [0.5 1.2]; % radius of the pillar/cylinder, depth below surface

% Set distance to iBox and its dimensions
ibox_dim = [0.6 1.6]; % dx, dy
ibox_dist = 12; % m
ibox_water = 0.1 ; %water below ibox in m

% Set section to be computed/carried out. 1 = on, 0 = off. Settings for
% individual sections can be found below. SOME sections require DEM!
section.area = 0;           % Area covered by building/pillar
section.umbrella = 0;       % Gravity effect of area bellow building
section.integ_initial = 0;  % Gravity effect as function of integration radius = inputs.
section.integ_sum = 0;      % Read + Sum the results for integration radius analysis
section.integ_plot = 0;     % Plot the results
section.sensor_height = 0;  % Sensitivity to sensor height
section.ibox_effect = 0;    % compute effect of water below ibox

% Set if Matlab should be closed after finishing the compuation
close_after_computing = 0;

font_size = 11; % for plots

%% Area covered by building/pillar
% Umbrella effect for this section (different value may be used in other
% sections!!!
umbrella_depth = 2; % depth in m
umbrella_vol = 10;  % mm. For this value will be computed the effect. 

if section.area == 1
    fprintf('Umbrella effect: area sealed by building or pillar\n');
    for i = 1:length(sites)
        if strcmp(char(sites(i)),'igrav')
            area.(char(sites(i))) = pi*pillar.(char(sites(i)))(1)^2; % area = pi * radius^2
        else
            area.(char(sites(i))) = build.(char(sites(i)))(1)*build.(char(sites(i)))(2); % x * y
            area.(char(sites(i))) = build.(char(sites(i)))(1)*build.(char(sites(i)))(2); % x * y
        end
        % Write output to command line
        fprintf('Area covered: %s = %4.1f m^2\n',...
                char(sites(i)),area.(char(sites(i))));
    end
end

%% Umbrella effect of iBox
if section.ibox_effect == 1
    dg_ibox = sorokin([0,0,sensor.igrav],0,ibox_dist,0,1000,ibox_water,ibox_dim)*10;
    fprintf('\nUmbrella effect iBox: %.3f m water below iBox = %.3f nm/s^2\n',...
        ibox_water,dg_ibox);
end
%% Gravity effect of area bellow building
if section.umbrella == 1
    fprintf('\nUmbrella effect building: gravity effect for depth of %3.1f m and %3.1f mm\n',...
            umbrella_depth,umbrella_vol);
    for i = 1:length(sites)
        if strcmp(char(sites(i)),'igrav')
             g_umbrella.(char(sites(i))) = cylinderEffect(sensor.(char(sites(i))),...height above the cyliner/ubrella zone = sensor height
                            pillar.(char(sites(i)))(1),... %  radius of the cylinder
                            umbrella_depth,umbrella_vol/umbrella_depth)*10; % depth and density (10 = 1Vol%) of the umbrella effect. *10 = uGal >> nm/s^2
        else
            g_umbrella.(char(sites(i))) = ...
                    sorokin([coor_relative.(char(sites(i)))(1),coor_relative.(char(sites(i)))(2),coor_relative.(char(sites(i)))(3) + sensor.(char(sites(i)))],... % position of the sensor
                             build.(char(sites(i)))(1)/2,build.(char(sites(i)))(2)/2,build.(char(sites(i)))(3),... position==centre of the building
                             umbrella_vol/umbrella_depth,umbrella_depth,... density (10 = 1Vol%) and depth of the umbrella effect
                             [build.(char(sites(i)))(1),build.(char(sites(i)))(2)])*10;% - ... % building dimennsions. *10 = uGal >> nm/s^2. Minus pillar as there cannot by water!
        end
        % Write output to command line
        fprintf('Gravity effect: %s = %4.2f nm/s^2\n',...
                char(sites(i)),g_umbrella.(char(sites(i))));
    end
end
%% Gravity effect as function of integration radius
% Set integration radius values (set it outside 'if section.integ_initial' 
% because these values will/could be used in following section)
% Warning: DO NOT SET minimum integration radius below the iGrav pillar!!
integration_radius = sort([1:5,10:10:11900]); % m
% integration_radius = sort([logspace(1,4,40),1:5,50,300:100:900,11120,11900]); % m
% mLocal settings
soil_moisture_volPerc = 10;                                                 % Compute the gravity response to soil moisture layers = X Vol%
interp_height = 1;                                                          % use interpolated height (+sensor height) (1=on/0=off)
coor_system = 1;                                                            % 1 = x, y.  2 = lon, lat.
exclusion_depth = 2;                                                        % depth of the umbrella effect
exclusion_polygon = fullfile('..','..','Data','DEM','Sample_DEM.asc');      % file name of exclusion polygon. The cylinder for iGrav will be subtracted separately
dem{1} = fullfile('..','..','Data','DEM','Sample_DEM.asc');                 % DEM for the first zone (on hygra server)
dem{2} = fullfile('..','..','Data','DEM','Sample_DEM.asc');                 % DEM for the second zone
dem{3} = fullfile('..','..','Data','DEM','Sample_DEM.asc');                 % DEM for the third zone
dem_interp = {0.1, 1, 10};                                                  % resolution of  the {first,second,third} zone
radius = {50,300,11900};                                                   % max.integration radius of  the {first,second,third} zone
% Set soil moisture layers
sm_layer.start = 0:1:6;                                                     % soil moiture layer starting depth
sm_layer.stop  = 1:1:7;                                                     % soil moisture layer ending depth

if section.integ_initial == 1
    fprintf('\nIntegration radius analysis...\n');
    % Open a waitbar to see the progress
    h = waitbar(0,'Preparing integation radius analysis...');

    % Prepare inputs for mLocal toolbox. The is common step for all three
    % gravimeters
    for i = 1:length(sm_layer.start)                                            % prepare input data for each layer
        sm{i,1} = 1;                                                            % switch 1 = layer on, 0 = layer off
        sm{i,2} = sm_layer.start(i);                                            % starting depth
        sm{i,3} = sm_layer.stop(i);                                             % end depth
        if sm_layer.start(i) < exclusion_depth                                  % compute umbrella effect only up to given depth
            sm{i,4} = exclusion_polygon;                                        % full file name to exclusion polygon
        else
            sm{i,4} = [];                                                       % no exclusion
        end
        sm{i,5} = fullfile('TEMP','Temporary_SM_input_VolPerc.txt'); % Input file name. This file will be created automatically for given soil moisture (see soil_moisture_volPerc variable)!
        sm{i,6} = 2;                                                            % data column (in file sm{i,5})
    end
    % Create file with soil moisture in mLocal file format: timeInMatlabFormat soilMoisture
    dlmwrite(sm{i,5},[datenum(clock),soil_moisture_volPerc],'delimiter',' ','precision','%14.4f'); % write SM to a txt file
    % Groundwater and snow effect
    gw = {0,[],0,0,[],0,0};                                                     % no Groundwater effect
    snow = {0,[],0,0,0,[],0,0};                                                 % no snow effect
    output_file{2} = 1;output_file{3} = 0;output_file{4} = 0;                   % txt output on
    output_file{5} = 1;                                                         % matrix ON - will be used for further computation

    % Run loop for all gravimeters
    for i = 1:length(sites)
        waitbar(i/length(sites),h,sprintf('Computing gravity response for integration radius (%s)',char(sites(i))));
        % Create output file name
        output_file{1} = fullfile('TEMP',sprintf('%s_forIntegRadius_%04dlayers_radius%04dm_umbella%04dcm.txt',...
                                        char(sites(i)),length(sm_layer.start),length(integration_radius),round(exclusion_depth*100))); 
        % Compute (call mLocal toolbox)
        mLocal_response([coor.(char(sites(i))),sensor.(char(sites(i)))],...         % input coordinates and sensor height
            interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file);% Use settings declared above
    end                   
    % Remove used variables
    clear interp_height coor_system dem dem_interp radius sm gw snow i
    close(h);
end

%% Read + Sum the results for integration radius analysis
if section.integ_sum == 1
    % The results for all sites, all layers and each integration radius will
    % be saved in following file:
    output_integ_radius_all = fullfile('..','..','Data','ForwardModel',sprintf('All_IntegRadius_%04dlayers_radius%04dm_umbella%04dcm.mat',...
                                        length(sm_layer.start),length(integration_radius),round(exclusion_depth*100)));
    
    % Open a waitbar to see the progress
    h = waitbar(0,'Reading integration radius input data...');
    c = 1; % just to show the progress (in waitbar)

    for i = 1:length(sites)
        % Declare output variable
        dg_integRadius.(char(sites(i))) = 0;

        % create file name for loading. Use the same code as in mLocal loop but
        % with a different file suffix (mat)
        file_name_mat = fullfile('TEMP',sprintf('%s_forIntegRadius_%04dlayers_radius%04dm_umbella%04dcm.mat',...
                                        char(sites(i)),length(sm_layer.start),length(integration_radius),round(exclusion_depth*100))); 
        % Read the matlab structure aray (see mLocal for output structure)
        all_layers = importdata(file_name_mat);

        % Run loop for all SM layers
        for l = 1:length(sm_layer.start)
            % Get gravity response for current layer. This is only done to
            % avoid using 3 dimensions (= 3x indices).
            zone1 = all_layers.sm_response_zone1.dg(:,:,l);
            zone2 = all_layers.sm_response_zone2.dg(:,:,l);
            zone3 = all_layers.sm_response_zone3.dg(:,:,l);
            % Run loop for all integration radius values
            for r = 1:length(integration_radius)
                waitbar((c)/(length(sites)*length(sm_layer.start)*length(integration_radius)),h,sprintf('Integrating gravity effect... (%s)',char(sites(i))));
                % Sum all grid cells within the current radius values and convert
                % the result to nm/s^2 
                dg_integ_zone1 = sum(sum(zone1(all_layers.sm_response_zone1.dist<=integration_radius(r))))*10^9;
                dg_integ_zone2 = sum(sum(zone2(all_layers.sm_response_zone2.dist<=integration_radius(r))))*10^9;
                dg_integ_zone3 = sum(sum(zone3(all_layers.sm_response_zone3.dist<=integration_radius(r))))*10^9;
                % Subtract the pillar effect only for iGrav and only if current
                % layer if affected by umbrella effect (=exlusion_depth)
                if strcmp(char(sites(i)),'igrav') && sm_layer.start(l) < exclusion_depth 
                    dg_integ_zone1 = dg_integ_zone1 - cylinderEffect(sensor.(char(sites(i)))+sm_layer.start(l),...height above the cyliner/ubrella zone = sensor height + distance to layer. Do not add pillar length like in 'g_umbrella'!; *10 => nm/s^2
                            pillar.(char(sites(i)))(1),... %  radius of the cylinder
                            sm_layer.stop(l)-sm_layer.start(l),10)*10; % cylinder/layer thickness, density (10 = 1 Vol%), *10 => nm/s^2
                end
                % Sum all zones and multiply with given soil moiture value
                % as the 'mat' results do not count the 'Temporary_SM_input_VolPerc.txt'
                % values into account (only for 'txt' outputs).
                % Columns store layers and rows radius values
                dg_integRadius.(char(sites(i)))(r,l) = (dg_integ_zone1 + dg_integ_zone2 + dg_integ_zone3)*soil_moisture_volPerc;
                c = c + 1;
            end
        end
        clear all_layers
    end
    % Save results
    dg_integRadius.radius = integration_radius';
    dg_integRadius.start = sm_layer.start;
    dg_integRadius.stop = sm_layer.stop;
    dg_integRadius.volPerc = soil_moisture_volPerc;
    save(output_integ_radius_all,'dg_integRadius');
    clear c;
    close(h);
end

%% Plot the result: integration radius
if section.integ_plot == 1
    % Load results from previous section
    output_integ_radius_all = fullfile('..','..','Data','ForwardModel',sprintf('All_IntegRadius_%04dlayers_radius%04dm_umbella%04dcm.mat',...
                                        length(sm_layer.start),length(integration_radius),round(exclusion_depth*100)));
    load(output_integ_radius_all);   
    % Create new figure/window. Use fixed Position (location+size) so identical
    % plots are created regardles of computed resolution.
    f1 = figure('Units','pixels','Position',[200 400 800, 280],'PaperPositionMode','auto');
    % Set plotting options
    xlimits = [1 dg_integRadius.radius(end)];             % limits on X axis
    xtick = [1 10 100 1000 10000];  % ticks on X axis
    ylimits = [0 0.55];             % limits on Y axis
    ytick = 0:0.1:0.6;              % ticks on X axis
    % Temporary variables for axes position
    if length(sites) == 3
        xstart = 0.065;                 % defines the staring (left) point of axes
        xwidth = 0.29;
        xstep = 0.315;
    else
        xstart = 0.08;
        xwidth = 0.42;
        xstep = 0.47;
    end
    % Run loop for all sites
    for i = 1:length(sites)
        % Convert results (in nm/s^2/Vol%) to nm/s^2/mm
        for j = 1:length(dg_integRadius.start)
            % mm = volume * density
            dg_integRadius.(sites{i})(:,j) = dg_integRadius.(sites{i})(:,j)./((dg_integRadius.stop(j) - dg_integRadius.start(j))*(soil_moisture_volPerc*10));
        end
        clear j

        % Open new axes for each site
        ax = axes('Units','normalized','Position',[xstart,0.2,xwidth,0.72]);
        xstart = xstart + xstep;    % add to move to next axes
        % Plot the data in log (X) axis
        semilogx(dg_integRadius.radius,dg_integRadius.(char(sites(i)))(:,1:length(dg_integRadius.start)))
        % Depending on the plot/site, add Y label and legend
        if i == 1
            ylab = 'gravity effect per 1 mm (nm.s^{-2})';
            ytl = ytick;
        else
            ylab = '';
            ytl = 0;                % do not show ticks for other than first site
        end
        xlab = 'integration radius (m)';
        % Create legend and show in in the last plot
        if i == length(sites)
            % Remove 'leg' variable first so it is possible to use indices
            clear leg
            for j = 1:length(dg_integRadius.start)
                leg(j) = {sprintf('%2.1f-%2.1f',dg_integRadius.start(j),dg_integRadius.stop(j))};
            end
        else
            leg = [];
        end
        % Add X% threshold
        threshold = 99;
        if strcmp(sites{i},'igrav');
            % Compute the threshold for the lowest layer = highest effect
            temp = dg_integRadius.(char(sites(i)))(:,length(dg_integRadius.start));
            temp_perc = temp./max(temp);
            % Find where Threshold value is reached
            r = find(temp_perc*100 >= threshold);
            % Write to command promtp
            fprintf('Integration radius Plot iGrav: %.1f%% threshold at %4.2f km\n',...
                threshold,dg_integRadius.radius(r(1))/1000);
            % Write [min, max] gravity effect range
            temp = dg_integRadius.(char(sites(i)))(end,1:length(dg_integRadius.start));
            fprintf('Integration radius Plot iGrav: min = %.2f nm/s^2 max = %4.2f nm/s^2\n',...
                min(temp),max(temp));
            % Add to plot
            hold on
            plot([dg_integRadius.radius(r(1)) dg_integRadius.radius(r(1))],ylimits,'k-')
            text(dg_integRadius.radius(r(1)),.10,sprintf('%.0f%% threshold',threshold),'VerticalAlignment','top','Rotation',90,'FontSize',font_size-1);
            hold off
        end
        % Use mm_setaxes to set all required parameters
        [~,~,~,leg_o] = mm_setaxes(ax,'font_size',font_size,'xlabel',xlab,'ylabel',ylab,...
                       'legend',leg,'xlim',xlimits,'ylim',ylimits,'xticklabel',xtick,'xtick',xtick,'ytick',ytick);
        % Place the legend in suitable place
        if ~isempty(leg_o)
            set(leg_o,'Location','southeast'); 
        end
        % Add gravimeters/site name
        text(60,ylimits(end)-5,char(sites(i)),'FontSize',font_size)
    end
end

%% Sensitivity to sensor height (iGrav only)
if section.sensor_height == 1
    fprintf('\nSensitivity to sensor height...\n');
    % Open a waitbar to see the progress
    h = waitbar(0,'Preparing snow cover analysis...');
    c = 1; % just to show the progress (in waitbar)

    soil_moisture_volPerc = 1;                                              % Compute the gravity response to soil moisture layers = X Vol%. 1 => 10 mm
    interp_height = 1;                                                      % use interpolated height (+sensor height) (1=on/0=off)
    coor_system = 1;                                                        % 1 = x, y.  2 = lon, lat.
    exclusion_depth = 2;                                                    % depth of the umbrella effect
    exclusion_polygon = fullfile('..','..','Data','DEM','Sample_Polygon.txt'); % file name of exclusion polygon. The cylinder for iGrav will be subtracted separately
    dem{1} = fullfile('..','..','Data','DEM','Sample_DEM.asc');    % DEM for the first zone (on hygra server)
    dem{2} = fullfile('..','..','Data','DEM','Sample_DEM.asc');           % DEM for the second zone
    dem{3} = fullfile('..','..','Data','DEM','Sample_DEM.asc');           % DEM for the third zone
    dem_interp = {0.2, 2, 20};                                              % resolution of  the {first,second,third} zone
    radius = {50,300,11900};                                                % max.integration radius of  the {first,second,third} zone
    % Set soil moisture layers
    sm_layer.start = [0];                                                 % soil moiture layer starting depth
    sm_layer.stop  = [1];   

    % Prepare inputs for mLocal toolbox. The is common step for all three
    % gravimeters
    for i = 1:length(sm_layer.start)                                        % prepare input data for each layer
        sm{i,1} = 1;                                                        % switch 1 = layer on, 0 = layer off
        sm{i,2} = sm_layer.start(i);                                        % starting depth
        sm{i,3} = sm_layer.stop(i);                                         % end depth
        if sm_layer.start(i) < exclusion_depth                              % compute umbrella effect only up to given depth
            sm{i,4} = exclusion_polygon;                                    % full file name to exclusion polygon
        else
            sm{i,4} = [];                                                   % no exclusion
        end
        sm{i,5} = fullfile('TEMP','Temporary_SM_input_VolPerc.txt'); % Input file name. This file will be created automatically for given soil moisture (see soil_moisture_volPerc variable)!
        sm{i,6} = 2;                                                            % data column (in file sm{i,5})
    end
    % Create file with soil moisture in mLocal file format: timeInMatlabFormat soilMoisture
    dlmwrite(sm{i,5},[datenum(clock),soil_moisture_volPerc],'delimiter',' ','precision','%14.4f'); % write SM to a txt file
    % Groundwater and snow effect
    gw = {0,[],0,0,[],0,0};                                                 % no Groundwater effect
    snow = {0,[],0,0,0,[],0,0};                                             % no snow effect
    output_file{2} = 1;output_file{3} = 0;output_file{4} = 0;               % txt output on
    output_file{5} = 0;                                                     % matrix ON - will be used for further computation

    % Run loop for all sensor height
    for i = 1:length(sensor.igrav_test)
        waitbar(c/length(sensor.igrav_test),h,sprintf('Sensitivity to sensor height (%3.0f cm)',100*sensor.igrav_test(i)));
        % Create output file name
        output_file{1} = fullfile('TEMP',sprintf('SensitivityToSensHeight_%3.0fcm_%04dlayers_radius%04dm_umbella%04dcm.txt',...
                                        sensor.igrav_test(i)*100,length(sm_layer.start),length(integration_radius),round(exclusion_depth*100))); 
        % Compute (call mLocal toolbox)
        mLocal_response([coor.(char(sites(end))),sensor.igrav_test(i)],...         % input coordinates and sensor height
            interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file);% Use settings declared above
        % Read results
        temp = load(output_file{1});
        % Subtract cylnder/pillar
        sensor_effect(i) = temp(4) - cylinderEffect(sensor.igrav_test(i)+sm_layer.start(1),...height above the cyliner/ubrella zone = sensor height + distance to layer. Do not add pillar length like in 'g_umbrella'!; *10 => nm/s^2
                            pillar.(char(sites(end)))(1),... %  radius of the cylinder
                            sm_layer.stop(1)-sm_layer.start(1),10)*10*soil_moisture_volPerc; % cylinder/layer thickness, density (10 = 1 Vol% = 10 mm), *10 => nm/s^2;
        c = c + 1;
    end                          
    % Remove used variables
    clear interp_height coor_system dem dem_interp radius sm gw snow i c
    close(h);
    % Output file name
    output_file = fullfile('..','..','Manuscript','Figures','Figure7.jpg');
    % Figure position
    figure_position = [200 400 390, 260];
    % Plot settings
    font_size = 11;
    % X and Y labels
    y_label = 'sensor height (m)';
    x_label = 'gravity effect per 10 mm (nm.s^{-2})';
    y_limits = [0 2];y_ticks = 0:0.5:2;
    x_limits = [3.9 4.8];x_ticks = 3.8:0.2:4.8;
    figure('Position',figure_position,'PaperPositionMode','auto');
    plot(sensor_effect,sensor.igrav_test,'linewidth',1)
    [~,~,~,leg_o] = mm_setaxes(gca,'fontsize',[font_size,font_size,font_size,font_size,font_size-1],...
                'xlabel',x_label,'ylabel',y_label,'xlim',x_limits,'ylim',y_limits,...
                'xtick',x_ticks,'ytick',y_ticks);
    print(gcf,output_file,'-r300','-djpeg');
end

%% Close matlab if required
rmpath(fullfile('..','MatlabLibrary','hydroGravityLib'));
rmpath(fullfile('..','MatlabLibrary','aux_fce'));
if close_after_computing == 1
    quit
end