%% Compute gravity effect of n layers
% Use this script to compute the gravity response to soil moisture change
% of 1Vol%. The output file can be then used for computation of gravity
% effect time series.
% Be patient, the computation will take some time!!! Reduce vertical
% resolution of sm_layers to accelerate the computation
% 
addpath(fullfile('..','MatlabLibrary','aux_fce'));
close all
clear
clc
%% Main settings
% output_file_name = fullfile('..','..','Data','Output','ForwardModel','iGrav_gravity_response_1VolPerc_10m_deep_VertResolution_0p01m.txt');
output_file_name = fullfile('..','..','Data','ForwardModel','SG030_gravity_response_1VolPerc_10m_deep_VertResolution_0p01m.txt');

% Input coordinates of the gravimeter. Height at surface. Sensor height
% will be added (m)
% coordinates = [4564082.00 5445669.70 609.755]; % iGrav
coordinates = [4564041.87 5445662.88 606.471]; % SG030
% sensor_height = 1.05; % m (iGrav)
sensor_height = 0.25; % m (SG030)

% Digital elevation models
% resolution of  the {first,second,third} zone
dem_interp = {0.2, 2, 20};     
% max.integration radius of  the {first,second,third} zone
radius = {50,300,11900};  
% Input file names
dem{1} = fullfile('..','..','Data','DEM','Sample_DEM.asc');  % DEM for the first zone (fines)
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 (lowest approximation)

% Umbrella space
exclusion_depth = 2;                                                        % depth of the umbrella effect
% Exclusion polygon or radius of the cylinder that will be subtracted
% exclusion_polygon = 0.5 ;% iGrav (double or string) 
exclusion_polygon = fullfile('..','..','Data','DEM','Sample_Polygon.txt'); % SG030

% Soil layers settings
soil_moisture_volPerc = 1;                                                  % 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.

% Set soil moisture layers
sm_layer.start = 0.00:0.01:10.09;                                           % soil moiture layer starting depth
sm_layer.stop  = 0.01:0.01:10.10;  

%% Prepare inputs for mLocal toolbox.
% Generate inputs for soil layers
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
    % Compute umbrella effect only up to given depth and only if string =
    % exclusion polygon is set. Otherwise, compute the cylinder effect
    % outside mLocal
    if (sm_layer.start(i) < exclusion_depth) && ischar(exclusion_polygon)                                 
        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. Use soilMoisture*10 and then devide by 100 to overcome the
% precission issues
dlmwrite(sm{i,5},[datenum(clock),soil_moisture_volPerc*100],'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

% Create output file name
output_file{1} = fullfile('TEMP','Temporary_SM_output_VolPerc.txt'); 
output_file{2} = 1;                                                         % txt output on
output_file{5} = 0;output_file{3} = 0;output_file{4} = 0;                   % matrix, tsf and xls off



%% Call mLocal
% Compute (call mLocal toolbox)
mLocal_response([coordinates,sensor_height],...         % input coordinates and sensor height
    interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file);% Use settings declared above

% Remove mLocal toolbox from search path                           
% Remove temporary file
delete(sm{i,5});
close 
    
%% Read the results
data = load(output_file{1});
% Devity by 100 to get back 1 VolPerc. See row 67
data(:,4:end) = data(:,4:end)./100;
% Remove the used file
delete(output_file{1});
% Declare variable
g_sm = zeros(length(sm_layer.start),4);
for i = 1:length(sm_layer.start)
    g_sm(i,2:4) = data(7+i:length(sm_layer.start):end);
    % If iGrav, subtract effect of the pillar/cylinder
    if (sm_layer.start(i) < exclusion_depth) && ~ischar(exclusion_polygon)  
        temp_g = cylinderEffect(sensor_height+sm_layer.start(i),...       height above the cyliner/ubrella zone = sensor height + distance to the layer
                            exclusion_polygon,... %     radius of the cylinder
                            abs(sm_layer.start(i)-sm_layer.stop(i)),... total length of the cylinder = i-th layer
                            10*soil_moisture_volPerc)*10; % density *10 => kg/m^3, *10 => nm/s^2
    else
        temp_g = 0;
    end
    % Correct for pillar/ubrella effect
    g_sm(i,2) = g_sm(i,2) - temp_g;
end
% Sum the total effect
g_sm(:,1) = sum(g_sm(:,2:end),2);

%% Save the results
fid = fopen(output_file_name,'w');
fprintf(fid,'%% Input coordinates: %.3f, %.3f, %.3f + %.3f m\n',coordinates(1),coordinates(2),coordinates(3),sensor_height);
fprintf(fid,'%% mLocal settings:\n');
fprintf(fid,'%% DEM interpolation zone 1, 2, 3: %.2f, %.2f, %.2f m\n',dem_interp{1},dem_interp{2},dem_interp{3});
fprintf(fid,'%% DEM radius zone 1, 2, 3: %.2f, %.2f, %.2f m\n',radius{1},radius{2},radius{3});
for i = 1:3
    fprintf(fid,'%% DEM zone %02d input: %s\n',i,dem{i});
end
fprintf(fid,'%% Interpolate DEM (01=yes): %02d\n',interp_height);
fprintf(fid,'%% Umbrella depth: %.2f m\n',exclusion_depth);
if ischar(exclusion_polygon)
	fprintf(fid,'%% Umbrella polygon: %s m\n',exclusion_polygon);
else
	fprintf(fid,'%% Umbrella cylinder radius: %.2f m\n',exclusion_polygon);
end
fprintf(fid,'%% Input soil moisture: %.2f VolPercent\n',soil_moisture_volPerc);
fprintf(fid,'%% Results: depths in m and gravity effect in nm/s^2:\n');
fprintf(fid,'%% StartDepth EndDepth Sum Zone1 Zone2 Zone3\n');
for i = 1:length(sm_layer.start)
    fprintf(fid,'%6.3f %6.3f %8.5f %8.5f %8.5f %8.5f\n',sm_layer.start(i),sm_layer.stop(i),g_sm(i,1),g_sm(i,2),g_sm(i,3),g_sm(i,4));
end
fclose(fid);


