function out = mLocal_response(Input,interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file)
%MLOCAL_RESPONSE Function computes the gravity response to SM, GW and snow
% Input variables:
%   Input        ...    Gravimeter coordinates [x,y,height,sensor height]
%                       or [longitude(deg), latitude(deg),height,sensor 
%                       height]. The sensor height will be added to height.
%                       Only one point per computation allowed (vector).
%   interp_height...    If == 1 than the gravimeter height will be
%                       interpolated using first digital elevation model +
%                       the sensor height (scalar, 0|1).
%   coor_system  ...    Coordinate system switch
%                       == 1 => x,y,height (m)
%   dem          ...    Full file names of DEMs stored in Matlab cell 
%                       format, i.e. dem{1} = DEM for zone 1
%                                    dem{2} = DEM for zone 2
%                                    dem{3} = DEM for zone 3
%   dem_interp   ...    Information related to re-interpolation of loaded
%                       DEM (matlab cell).
%                       dem_interp{1} = new resolution of zone 1 (m|deg)
%                       dem_interp{2} = new resolution of zone 2 (m|deg)
%                       dem_interp{3} = new resolution of zone 3 (m|deg)
%                       if dem_interp{i} == [], than no interpolation of
%                       i-th zone will be performed. The digital elevation
%                       model can be stored in *.mat, *.asc or *.grd
%                       format.
%   radius       ...    Radius of computation zones (matlab cell)
%                       radius{1} = max.integration radius of zone 1(m|deg)
%                       radius{2} = max.integration radius of zone 1(m|deg)
%                       radius{3} = max.integration radius of zone 1(m|deg)
%   sm           ...    All information related to soil moisture effect
%                       computation (matlab cell). The rows define the
%                       number of sm layers, whereas the columns define
%                       their parameters. The maximal number of layers is
%                       not limited.
%                       sm{1,1} = layer 1 on/off == 1|0
%                       sm{1,2} = layer 1 starting depth (m)(0 for surface)
%                       sm{1,3} = layer 1 ending depth (m)
%                       sm{1,4} = full file name of the exclusion polygon.
%                        All grid cells (layer1) inside of this polygon will 
%                        be removed from the integration. Polygons can be
%                        stored either in *.shp or *.txt format. The *.shp
%                        file format allows the loading of various polygons
%                        stored in one file (requires *.X and *.Y layer for
%                        both lon/lat or x/y coordinate system). If the
%                        *.txt file contains only one number, this number
%                        will be used as an exclusion radius (m|deg). To
%                        store more polygons in one *.txt file, use NaNs to
%                        separate them.
%                       sm{1,5} = full file name of the input file. The
%                        time series can be stored in *.txt, *.xls or 
%                        *.tsf. The time vector (matlab time = datenum) 
%                        must be located in the first column except of the 
%                        *.tsf format where the time is stored in first six
%                        columns. [] if not selected.
%                       sm{1,6} = data column. This value points to the
%                        column with soil moisture (Vol%). In case of *.tsf,
%                        this value determines the channel not the column.
%                       sm{2,1} = layer 2 on/off == 1|0
%                       ...
%   gw           ...    All information related to the groundwater effect
%                       computation (matlab cell). One row cell.
%                       gw{1} = groundwater on/off == 1|0
%                       gw{2} = full file name of the exclusion polygon.
%                        All grid cells  inside of this polygon will be
%                        removed from the integration. Polygons can be
%                        stored either in *.shp or *.txt format. The *.shp
%                        file format allows the loading of various polygons
%                        stored in one file (requires *.X and *.Y layer for
%                        both lon/lat or x/y coordinate system). If the
%                        *.txt file contains only one number, this number
%                        will be used as a exclusion radius (m|deg).
%                       gw{3} = specific yield/effective porosity in %.
%                       gw{4} = computation method switch (1|2).
%                        == 1 the groundwater altitude will be loaded and
%                        the max. groundwater level will be determine
%                        using these altitude values.
%                        == 2 the groundwater depth will be loaded and the
%                        max. groundwater is computed using DEMs and given
%                        gw. depths, i.e. the groundwater copies the
%                        terrain.
%                       gw{5} = used only if gw{4} == 1. In such case, this
%                        value will be used to compute the minimal distance
%                        (in m) between the groundwater and terrain (max.
%                        groundwater is given by loaded values). This zone
%                        represents the unsaturated layer.
%                       gw{6} = full file name of the input file. The
%                        time series can be stored in *.txt, *.xls or 
%                        *.tsf. The time vector (matlab time = datenum) 
%                        must be located in the first column except of the 
%                        *.tsf format where the time is stored in first six
%                        columns.
%                       gw{7} = data column. This value points to the
%                        column with groundwater depth or altitude (m). 
%                        In case of *.tsf, this value determines the channel 
%                        not the column.
%   snow         ...    All information related to the snow effect
%                       computation (matlab cell). One row cell.
%                       snow{1} = snow effect on/off == 1|0
%                       snow{2} = full file name of the exclusion polygon.
%                        All grid cells  inside of this polygon will be
%                        removed from the integration. Polygons can be
%                        stored either in *.shp or *.txt format. The *.shp
%                        file format allows the loading of various polygons
%                        stored in one file (requires *.X and *.Y layer for
%                        both lon/lat or x/y coordinate system). If the
%                        *.txt file contains only one number, this number
%                        will be used as a exclusion radius (m|deg). To
%                        include the snow on the roof, modify the DEM.
%                       snow{3} = constant term used for the computation of
%                        snow density (kg/m^3). 
%                       snow{4} = snow height depended term used for the 
%                        computation of snow density (kg/m^3/m). The snow{4}
%                        and snow{5} values will be used only if now snow
%                        water equivalent is loaded (see snow{8}).
%                       snow{5} = snow gradient, i.e. snow height
%                        decrease/increase per 100 m (mm/100 m). Set this
%                        value to zero to neglect the altitudal snow
%                        dependency. Always positive.
%                       snow{6} = full file name of the input file. The
%                        time series can be stored in *.txt, *.xls or 
%                        *.tsf. The time vector (matlab time = datenum) 
%                        must be located in the first column except of the 
%                        *.tsf format where the time is stored in first six
%                        columns.
%                       snow{7} = data column. This value points to the
%                        column with snow height/depth (mm). 
%                        In case of *.tsf, this value determines the channel 
%                        not the column.
%                       snow{8} = data column. This value points to the
%                        column with snow water equivalent (mm). 
%                        In case of *.tsf, this value determines the channel 
%                        not the column.
%   output_file  ...    Information related to the output files (matlab 
%                       cell).
%                       output_file{1} = full file name (including
%                        extension)
%                       output_file{2} = 1 if *.txt file should be written
%                                      = 0 if *.txt file should not be 
%                                          written 
%                       output_file{3} = 0 if *.xls file should be written
%                       output_file{4} = 0 if *.tsf file should be written
%                       output_file{5} = 0 if *.mat file should be saved.
%                        This file contains response matrices (m/s^2).
% 
% Output:
%   out         ...     output matlab array containing all results.
% 
% Example 1. Compute the gravity response to two soil moisture layer,
% groundwater and snow + write to txt file:
%   interp_height = 1;
%   coor_system = 1;
%   dem{1} = 'DEMzone1.mat';dem{2} = 'DEMzone2.mat';dem{3} = 'DEMzone3.mat';
%   dem_interp = {0.1, 1, 10};
%   radius = {100,300,10000};
%   sm = {1,0,1,'exclusion_polyg.shp','input_sm1.txt',2};
%   sm{2,1} = 1;sm{2,2} = 0; sm{2,3} = 1; sm{2,4} = []; sm{2,5} = {'input_sm2.txt'}; sm{2,6} = 2;
%   gw = {1,[],10,1,'input_gw_altitude.txt',5,2};
%   snow = {1,'exclusion_building_roof.txt',148,108,0,'input_snow_depth.txt',2,0};
%   output_file{1} = 'Output_file.txt';output_file{2} = 1;output_file{3} = 0;output_file{4} = 0;output_file{5} = 0;
%   out = mLocal_response([4564214.0,5445520.0,610.2,0.25],...
%                          interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file)
% 
% Example 2. Compute the gravity response to one soil moisture layer  
%   within the radius of 300 m.
%   interp_height = 1;
%   coor_system = 2;
%   dem{1} = 'DEMzone1_WGS84.mat';dem{2} = 'DEMzone2_WGS84.mat';dem{3} = [];
%   dem_interp = {0.000002, 0.00002, 0};
%   radius = {100,300,10000};
%   sm = {1,0,1,'exclusion_polyg.shp','input_sm1.txt',2};
%   gw = {0,[],0,0,[],0,0};
%   snow = {0,[],0,0,0,[],0,0};
%   output_file{1} = 'Output_file.txt';output_file{2} = 1;output_file{3} = 0;output_file{4} = 0;output_file{5} = 0;
%   out = mLocal_response([12.87866,49.14354,610.2,0.25],...
%                          interp_height,coor_system,dem,dem_interp,radius,sm,gw,snow,output_file)
% 
% 
%                                                   M. Mikolaj, 30.09.2014
%                                                   mikolaj@gfz-potsdam.de

%% Load DEM + interpolate if required
set(findobj('Tag','text_status_mLocal'),'String','Gravity response: loading DEMs...');drawnow; % send message to GUI
% DEM1
dem_name = dem{1};                                                          % get the DEM name
switcher = dem_name(end-2:end);                                             % get the DEM file format (extension)
switch switcher                                                             % switch between file formats
    case 'mat'
        demt = importdata(dem_name);                                        % load *.mat file format
        if coor_system == 2                                                 % switch lon/lat layers to x/y if geodetic coor. are used
            demt.x = demt.lon;
            demt.y = demt.lat;
            demt.lon = [];demt.lat = [];
        end
    case 'asc'
        demt = mm_ascii2mat(dem_name);                               % load arc ascii format
end
if ~isempty(dem_interp{1})                                                  % interpolate only if required
    if coor_system == 1                                                     % cartesian coordinates
        [dem1.x,dem1.y] = meshgrid((Input(1)-radius{1}-dem_interp{1}*2):dem_interp{1}:(Input(1)+radius{1}+dem_interp{1}*2),...% max. distance given by integration radius
                                   (Input(2)-radius{1}-dem_interp{1}*2):dem_interp{1}:(Input(2)+radius{1}+dem_interp{1}*2));
    else                                                                    % gedetic coordinates
        [dem1.x,dem1.y] = meshgrid(min(min(demt.x)):dem_interp{1}:max(max(demt.x)),... % must be DEM_min to DEM_max because radius (radius{i}) is not not the same as latitude/longitude difference!
                                   min(min(demt.y)):dem_interp{1}:max(max(demt.y)));
    end
    dem1.height = interp2(demt.x,demt.y,demt.height,dem1.x,dem1.y);clear demt dem_name switcher % interpolate to new grid and remove used DEM
else
    dem1 = demt;clear demt dem_name switcher                                % copy if the re-sampling is not required
end
% Interpolate gravimeter height if required
InputUse = Input(1:3);                                                      % declare the InputUse variable (will be used for compuation)
if interp_height == 1
    InputUse(3) = interp2(dem1.x,dem1.y,dem1.height,Input(1),Input(2)) + Input(4); % interpolate height and add the sensor height
else
    InputUse(3) = Input(3)+Input(4);                                        % add the sensor height only
end
% Check if Gravimeter + integration radius is located within the DEM
idt = inpolygon(Input(1),Input(2),[min(min(dem1.x)),max(max(dem1.x)),max(max(dem1.x)),min(min(dem1.x))],...
                                  [min(min(dem1.y)), min(min(dem1.y)), max(max(dem1.y)), max(max(dem1.y))]);
if ~idt
	warning_message = 'Warning: the computation point is not within the DEM of zone 1!';
    fprintf('%s\n',warning_message);
end
if coor_system == 1                                                         % Secondly compute the minimal distance to DEM border
    distt(1) = sqrt(0 + (Input(2)-min(min(dem1.y)))^2);                     % distance to the lower boundary - DEM of zone1
    distt(2) = sqrt(0 + (Input(2)-max(max(dem1.y)))^2);                     % distance to the upper boundary - DEM of zone1
    distt(3) = sqrt((Input(1)-min(min(dem1.x)))^2 + 0);                     % distance to the left boundary - DEM of zone1
    distt(4) = sqrt((Input(1)-max(max(dem1.x)))^2 + 0);                     % distance to the right boundary - DEM of zone1
else
    distt(1) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,min(min(dem1.y)).*pi/180)*180/pi; % compute spherical distance (Input, not InputUse)                   % distance to the lower boundary - DEM of zone1
    distt(2) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,max(max(dem1.y)).*pi/180)*180/pi; % distance to the upper boundary - DEM of zone1
    distt(3) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,min(min(dem1.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the left boundary - DEM of zone1
    distt(4) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,max(max(dem1.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the right boundary - DEM of zone1
end
for di = 1:4
   if distt(di) < radius{1}
      warning_message = 'Warning: the integration radius of zone 1 is beyond the selected DEM!'; 
      fprintf('%s\n',warning_message);
   end
end
clear distt di idt
% Transform WGS84 latitude to sphere if required
if coor_system == 2
    dem1.ySphere = mGlobe_elip2sphere(dem1.x*pi/180,dem1.y*pi/180)*180/pi;  % latitude on sphere (input and output in radians, therefore *180/pi to deg)
    InputUse(2) = mGlobe_elip2sphere(InputUse(1)*pi/180,InputUse(2)*pi/180)*180/pi; % Computation point: latitude on sphere (POINT deg)
    % Calc new grid resolution
    resol1lat = abs(dem1.ySphere - circshift(dem1.ySphere,[1 0]));          % approx. grid resolution in latitude direction (on sphere)
    resol1lat(1,:) = resol1lat(2,:);                                        % correct the edges
    resol1lat(end,:) = resol1lat(end-1,:);
end
dem1.height(isnan(dem1.height)) = 9e+36;                                    % Remove NaNs
resol1 = [abs(dem1.x(1,1)-dem1.x(1,2)),abs(dem1.y(1,1)-dem1.y(2,1))];       % DEM resolution (not valid for latitude direction)

% DEM2
if ~isempty(dem{2})                                                         % Load second DEM only if selected
    dem_name = dem{2};                                                      % do the same as for dem1
    switcher = dem_name(end-2:end);
    switch switcher
        case 'mat'
            demt = importdata(dem_name);
            if coor_system == 2
                demt.x = demt.lon;
                demt.y = demt.lat;
                demt.lon = [];demt.lat = [];
            end
        case 'asc'
            demt = mm_ascii2mat(dem_name);
    end
    if ~isempty(dem_interp{2})
        if coor_system == 1
            [dem2.x,dem2.y] = meshgrid((Input(1)-radius{2}-dem_interp{2}*2):dem_interp{2}:(Input(1)+radius{2}+dem_interp{2}*2),...
                                       (Input(2)-radius{2}-dem_interp{2}*2):dem_interp{2}:(Input(2)+radius{2}+dem_interp{2}*2));
        else
            [dem2.x,dem2.y] = meshgrid(min(min(demt.x)):dem_interp{2}:max(max(demt.x)),...
                                       min(min(demt.y)):dem_interp{2}:max(max(demt.y)));
        end
        dem2.height = interp2(demt.x,demt.y,demt.height,dem2.x,dem2.y);clear demt dem_name switcher
    else
        dem2 = demt;clear demt dem_name switcher
    end
    % Check if Gravimeter + integration radius is located within the DEM
    idt = inpolygon(Input(1),Input(2),[min(min(dem2.x)),max(max(dem2.x)),max(max(dem2.x)),min(min(dem2.x))],...
                                      [min(min(dem2.y)), min(min(dem2.y)), max(max(dem2.y)), max(max(dem2.y))]);
    if ~idt
        warning_message = 'Warning: the computation point is not within the DEM of zone 2!';
        fprintf('%s\n',warning_message);
    end
    if coor_system == 1                                                     % Secondly compute the minimal distance to DEM border
        distt(1) = sqrt(0 + (Input(2)-min(min(dem2.y)))^2);                 % distance to the lower boundary - DEM of zone1
        distt(2) = sqrt(0 + (Input(2)-max(max(dem2.y)))^2);                 % distance to the upper boundary - DEM of zone1
        distt(3) = sqrt((Input(1)-min(min(dem2.x)))^2 + 0);                 % distance to the left boundary - DEM of zone1
        distt(4) = sqrt((Input(1)-max(max(dem2.x)))^2 + 0);                 % distance to the right boundary - DEM of zone1
    else
        distt(1) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,min(min(dem2.y)).*pi/180)*180/pi; % compute spherical distance (Input, not InputUse)                   % distance to the lower boundary - DEM of zone1
        distt(2) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,max(max(dem2.y)).*pi/180)*180/pi; % distance to the upper boundary - DEM of zone1
        distt(3) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,min(min(dem2.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the left boundary - DEM of zone1
        distt(4) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,max(max(dem2.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the right boundary - DEM of zone1
    end
    for di = 1:4
       if distt(di) < radius{2}
          warning_message = 'Warning: the integration radius of zone 2 is beyond the selected DEM!'; 
          fprintf('%s\n',warning_message);
       end
    end
    clear distt di idt
    dem2.height(isnan(dem2.height)) = 9e+36;
    resol2 = [abs(dem2.x(1,1)-dem2.x(1,2)),abs(dem2.y(1,1)-dem2.y(2,1))];
    if coor_system == 2
        dem2.ySphere = mGlobe_elip2sphere(dem2.x*pi/180,dem2.y*pi/180)*180/pi; % latitude on sphere (grid deg)
        % Calc new grid resolution
        resol2lat = abs(dem2.ySphere - circshift(dem2.ySphere,[1 0]));      % approx. grid resolution in latitude direction (on sphere)
        resol2lat(1,:) = resol2lat(2,:);
        resol2lat(end,:) = resol2lat(end-1,:);
    end
else
    resol2 = [];                                                            % set to empty if dem2 is not selected
    radius{2} = [];
end

% DEM3
if ~isempty(dem{3})                                                         % load dem3 (zone3) if required
    dem_name = dem{3};
    switcher = dem_name(end-2:end);
    switch switcher
        case 'mat'
            demt = importdata(dem_name);
            if coor_system == 2
                demt.x = demt.lon;
                demt.y = demt.lat;
                demt.lon = [];demt.lat = [];
            end
        case 'asc'
            demt = mm_ascii2mat(dem_name);
    end
    if ~isempty(dem_interp{3})
        if coor_system == 1
            [dem3.x,dem3.y] = meshgrid((Input(1)-radius{3}-dem_interp{3}*2):dem_interp{3}:(Input(1)+radius{3}+dem_interp{3}*2),...
                                       (Input(2)-radius{3}-dem_interp{3}*2):dem_interp{3}:(Input(2)+radius{3}+dem_interp{3}*2));
        else
            [dem3.x,dem3.y] = meshgrid(min(min(demt.x)):dem_interp{3}:max(max(demt.x)),...
                                       min(min(demt.y)):dem_interp{3}:max(max(demt.y)));
        end
        dem3.height = interp2(demt.x,demt.y,demt.height,dem3.x,dem3.y);clear demt dem_name switcher
    else
        dem3 = demt;clear demt dem_name switcher
    end
    % Check if Gravimeter + integration radius is located within the DEM
    idt = inpolygon(Input(1),Input(2),[min(min(dem3.x)),max(max(dem3.x)),max(max(dem3.x)),min(min(dem3.x))],...
                                      [min(min(dem3.y)), min(min(dem3.y)), max(max(dem3.y)), max(max(dem3.y))]);
    if ~idt
        warning_message = 'Warning: the computation point is not within the DEM of zone 3!';
        fprintf('%s\n',warning_message);
    end
    if coor_system == 1                                                     % Secondly compute the minimal distance to DEM border
        distt(1) = sqrt(0 + (Input(2)-min(min(dem3.y)))^2);                 % distance to the lower boundary - DEM of zone1
        distt(2) = sqrt(0 + (Input(2)-max(max(dem3.y)))^2);                 % distance to the upper boundary - DEM of zone1
        distt(3) = sqrt((Input(1)-min(min(dem3.x)))^2 + 0);                 % distance to the left boundary - DEM of zone1
        distt(4) = sqrt((Input(1)-max(max(dem3.x)))^2 + 0);                 % distance to the right boundary - DEM of zone1
    else
        distt(1) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,min(min(dem3.y)).*pi/180)*180/pi; % compute spherical distance (Input, not InputUse)                   % distance to the lower boundary - DEM of zone1
        distt(2) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,Input(1)*pi/180,max(max(dem3.y)).*pi/180)*180/pi; % distance to the upper boundary - DEM of zone1
        distt(3) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,min(min(dem3.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the left boundary - DEM of zone1
        distt(4) = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,max(max(dem3.x))*pi/180,Input(2)*pi/180)*180/pi; % distance to the right boundary - DEM of zone1
    end
    for di = 1:4
       if distt(di) < radius{3}
          warning_message = 'Warning: the integration radius of zone 3 is beyond the selected DEM!'; 
          fprintf('%s\n',warning_message);
       end
    end
    clear distt di idt
    dem3.height(isnan(dem3.height)) = 9e+36;
    resol3 = [abs(dem3.x(1,1)-dem3.x(1,2)),abs(dem3.y(1,1)-dem3.y(2,1))];
    if coor_system == 2
        dem3.ySphere = mGlobe_elip2sphere(dem3.x*pi/180,dem3.y*pi/180)*180/pi; % latitude on sphere (grid deg)
        % Calc new grid resolution
        resol3lat = abs(dem3.ySphere - circshift(dem3.ySphere,[1 0]));      % approx. grid resolution in latitude direction (on sphere)
        resol3lat(1,:) = resol3lat(2,:);
        resol3lat(end,:) = resol3lat(end-1,:);
    end
else
    resol3 = [];
    radius{3} = [];
end

%% Soil moisture (1%)
sm_exist = 0;sm_file_exist = 0;                                             % check if any SM layer is selected
for i = 1:size(sm,1)
    if sm{i,1} == 1 
        sm_exist = sm_exist + 1;                                            % add to get the final number of layers
        if ~isempty(sm{i,5})
            sm_file_exist = sm_file_exist + 1;                              % add to get the final number of SM input files                                  
        end
    end
end;clear i;

dg_sm_1perc(1:size(sm,1),1:3) = 0;                                          % prepare variable (grav.response to sm = 1Vol%: all 3 zones)
if output_file{5} == 1 && sm_exist >= 1                                     % prepare gravity response matrices (if required)
    sm_response_zone1.x = dem1.x;sm_response_zone1.y = dem1.y;              % x and y (lon/lat) coordinates
    sm_response_zone1.dg = 0;                                               % future gravity effect matrix (m/s^2)
    sm_response_zone1.input = 1;                                            % Constant input for SM (=1Vol%)
    sm_response_zone1.dg_units = 'm/s2';                                    % Output units
    sm_response_zone1.input_units = 'Vol%';
    sm_response_zone1.start = [];sm_response_zone1.stop = [];               % declare SM layer start/stop depth variables
    if coor_system == 1
        sm_response_zone1.dist = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2); % distance between gravimeter and grid cells
        sm_response_zone1.dist_units = 'm';                                 % for cartesian coordinates
    else
        sm_response_zone1.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi; % compute spherical distance (Input, not InputUse)
        sm_response_zone1.ySphere = dem1.ySphere;
        sm_response_zone1.dist_units = 'deg';                               % for geodetic coordinates
    end
   if ~isempty(dem{2})
        sm_response_zone2.x = dem2.x;sm_response_zone2.y = dem2.y;          % do the same for zone 2
        sm_response_zone2.dg = 0;
        sm_response_zone2.input = 1;
        sm_response_zone2.dg_units = 'm/s2';
        sm_response_zone2.input_units = 'Vol%';
        sm_response_zone2.start = [];sm_response_zone2.stop = [];
        if coor_system == 1
            sm_response_zone2.dist = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2); % distance in m
            sm_response_zone2.dist_units = 'm';
        else
            sm_response_zone2.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi; % compute spherical distance
            sm_response_zone2.ySphere = dem2.ySphere;
            sm_response_zone2.dist_units = 'deg';
        end
   end
   if ~isempty(dem{3})                                                      % do the same for zone 3
        sm_response_zone3.x = dem3.x;sm_response_zone3.y = dem3.y;
        sm_response_zone3.dg = 0;
        sm_response_zone3.input = 1;
        sm_response_zone3.dg_units = 'm/s2';
        sm_response_zone3.input_units = 'Vol%';
        sm_response_zone3.start = [];sm_response_zone3.stop = [];
        if coor_system == 1
            sm_response_zone3.dist = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2);
            sm_response_zone3.dist_units = 'm';
        else
            sm_response_zone3.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi; % compute spherical distance
            sm_response_zone3.ySphere = dem3.ySphere;
            sm_response_zone3.dist_units = 'deg';
        end
   end
end

for i = 1:size(sm,1)                                                        % loop for each SM layer
    if sm{i,1} == 1                                                         % compute only if current layer is selected
        set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: soil moisture layer %2.0f, zone 1...',i));drawnow; % send message to GUI
        start = sm{i,2};                                                    % current SM layer starting depth
        stop = sm{i,3};                                                     % current SM layer ending depth
        % ZONE 1
        density = dem1.height;density(:,:) = 10;                            % by default, all grid cells are included (10 == 1Vol%)
        density(dem1.height>=9e+35) = 0;                                    % set density to zero for NaN pixels
        resolution = resol1;                                                % spatial resolution of the current zone (not used for latitude)
        polyg_file = sm{i,4};                                               % exclusion polygon file name ([] if not selected). Used for all zones
        if ~isempty(polyg_file)                                             % load exclusion polygon if selected
            polyg_switcher = polyg_file(end-2:end);                         % file type switch
            switch polyg_switcher
                case 'txt'
                    polyg = load(polyg_file);                               % load TXT file
                    if numel(polyg) == 1                                    % if the TXT file contains only one number, than use this number as exclusion radius
                        switch coor_system                                  % switch between coor.systems
                            case 1
                                id = find(sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2)<=polyg); % find grid cells using plane distance (less or equal to given exclusion radius!!)
                            case 2
                                id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi<=polyg); % find grid cells using spherical distance (input in degrees!)
                        end
                    else
                        id = inpolygon(dem1.x,dem1.y,polyg(:,1),polyg(:,2)); % use the original latitude, i.e. not transpormed to spherical (~dem1.ySphere). Only one polygon in TXT file! X/lon in first row, Y/Lat in second
                    end
                    density(id) = 0;                                        % Exclude grid cells within loaded polygon
            end
            clear polyg p polyg_switcher id                                 % delete used variables (polygon will be loaded for each zone separaterly)
        end
        if coor_system == 1                                                 % For cartesian coordinates
            dist = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2);   % temp. distance between grid cells and Gravimeter
            density(dist>radius{1}) = 0;clear dist;                         % remove points/pixels out of zone 1 (zet density to zero)
            if output_file{5} == 0                                          % if output matrix is not required
                dg_sm_1perc(i,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height-start,density,abs(start-stop),resolution,1); % sum of all pixels
            else                                                            % if output matrix is required
                if size(sm_response_zone1.dg,1) == 1                        % for the first layer
                    sm_response_zone1.dg = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height-start,density,abs(start-stop),resolution,0); % matrix output
                    dg_sm_1perc(i,1) = sum(sum(sm_response_zone1.dg));      % sum all pixels
                else                                                        % for all other layers
                    sm_response_zone1.dg(:,:,end+1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height-start,density,abs(start-stop),resolution,0); % matrix output
                    dg_sm_1perc(i,1) = sum(sum(sm_response_zone1.dg(:,:,end))); % sum all pixels
                end
                sm_response_zone1.start = vertcat(sm_response_zone1.start,start); % current starting depth
                sm_response_zone1.stop = vertcat(sm_response_zone1.stop,stop); % current ending depth
            end
        else                                                                % for cartesian coordinates
            dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180); % compute spherical distance
            dist = dist.*180/pi;                                            % transform to degrees       
            density(dist>radius{1}) = 0;clear dist;                         % remove points/pixels out of zone 1 (set density to zero)
            if output_file{5} == 0                                          % do the same as for cartesian coordinates (using tesseroid function)
                dg_sm_1perc(i,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height-stop,abs(start-stop),density,resolution(1),resol1lat,6371008,1); 
            else
                if size(sm_response_zone1.dg,1) == 1
                    sm_response_zone1.dg = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height-stop,abs(start-stop),density,resolution(1),resol1lat,6371008,0);
                    dg_sm_1perc(i,1) = sum(sum(sm_response_zone1.dg));
                else
                    sm_response_zone1.dg(:,:,end+1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height-stop,abs(start-stop),density,resolution(1),resol1lat,6371008,0);
                    dg_sm_1perc(i,1) = sum(sum(sm_response_zone1.dg(:,:,end)));
                end
                sm_response_zone1.start = vertcat(sm_response_zone1.start,start);
                sm_response_zone1.stop = vertcat(sm_response_zone1.stop,stop);
            end
        end
        clear density
        
        % ZONE 2
        set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: soil moisture layer %2.0f, zone 2...',i));drawnow;
        if ~isempty(dem{2})
            density = dem2.height;density(:,:) = 10;                        % by default, all grid cells are included (10 = 1 Vol% to kg/m^3)
            density(dem2.height>=9e+35) = 0;                                % set density to zero for NaN pixels
            resolution = resol2;                                            % spatial resolution of second zone (except of latitude)
            if ~isempty(polyg_file)
                polyg_switcher = polyg_file(end-2:end);
                switch polyg_switcher
                    case 'txt'
                        polyg = load(polyg_file);
                        if numel(polyg) == 1
                            switch coor_system
                                case 1
                                    id = find(sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2)<=polyg);
                                case 2
                                    id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi<=polyg);
                            end
                        else
                            id = inpolygon(dem2.x,dem2.y,polyg(:,1),polyg(:,2));
                        end
                        density(id) = 0;	                                % Exclude grid cells within loaded polygon
                end
                clear polyg id polyg_switcher p                             % delete used variables
            end
            if coor_system == 1
                dist = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2); % temp. distance between pixels and Gravimeter
                density(dist<=radius{1}) = 0;                               % remove points/pixels out of zone 2 (set density to zero)
                density(dist>radius{2}) = 0;clear dist
                if output_file{5} == 0 
                    dg_sm_1perc(i,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height-start,density,abs(start-stop),resolution,1);
                else
                    if size(sm_response_zone2.dg,1) == 1
                        sm_response_zone2.dg = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height-start,density,abs(start-stop),resolution,0); % matrix output
                        dg_sm_1perc(i,2) = sum(sum(sm_response_zone2.dg));   % sum all pixels
                    else
                        sm_response_zone2.dg(:,:,end+1) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height-start,density,abs(start-stop),resolution,0); % matrix output
                        dg_sm_1perc(i,2) = sum(sum(sm_response_zone2.dg(:,:,end)));   % sum all pixels
                    end
                    sm_response_zone2.start = vertcat(sm_response_zone2.start,start);
                    sm_response_zone2.stop = vertcat(sm_response_zone2.stop,stop);
                end
            else
                dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180); % compute spherical distance
                dist = dist.*180/pi;                                        % transform to degrees                     
                density(dist<=radius{1}) = 0;                               % remove points/pixels out of zone 2 (set density to zero)
                density(dist>radius{2}) = 0;clear dist
                if output_file{5} == 0 
                    dg_sm_1perc(i,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height-stop,abs(start-stop),density,resolution(1),resol2lat,6371008,1);
                else
                    if size(sm_response_zone2.dg,1) == 1
                        sm_response_zone2.dg = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height-stop,abs(start-stop),density,resolution(1),resol2lat,6371008,0);
                        dg_sm_1perc(i,2) = sum(sum(sm_response_zone2.dg));
                    else
                        sm_response_zone2.dg(:,:,end+1) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height-stop,abs(start-stop),density,resolution(1),resol2lat,6371008,0);
                        dg_sm_1perc(i,2) = sum(sum(sm_response_zone2.dg(:,:,end)));
                    end
                    sm_response_zone2.start = vertcat(sm_response_zone2.start,start);
                    sm_response_zone2.stop = vertcat(sm_response_zone2.stop,stop);
                end
            end
            clear density
        end
        
        % ZONE 3
        set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: soil moisture layer %2.0f, zone 3...',i));drawnow;
        if ~isempty(dem{3})
            density = dem3.height;density(:,:) = 10;                        % by default, all grid cells are included
            density(dem3.height>=9e+35) = 0;                                % set density to zero for NaN pixels
            resolution = resol3;                                            % spatial resolution of third zone (except of latitude)
            if ~isempty(polyg_file)
                polyg_switcher = polyg_file(end-2:end);
                switch polyg_switcher
                    case 'txt'
                        polyg = load(polyg_file);
                        if numel(polyg) == 1
                            switch coor_system
                                case 1
                                    id = find(sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2)<=polyg);
                                case 2
                                    id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi<=polyg);
                            end
                        else
                            id = inpolygon(dem3.x,dem3.y,polyg(:,1),polyg(:,2));
                        end
                        density(id) = 0;                                    % Exclude grid cells within loaded polygon
                end
                clear polyg id polyg_switcher p                             % delete used variables                                               % delete used variables
            end
            if coor_system == 1
                dist = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2); % temp. plane distance between pixels and Gravimeter
                density(dist<=radius{2}) = 0;                               % remove points/pixels out of zone 3 (set density to zero)
                density(dist>radius{3}) = 0;clear dist
                if output_file{5} == 0 
                    dg_sm_1perc(i,3) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height-start,density,abs(start-stop),resolution,1);
                else
                    if size(sm_response_zone3.dg,1) == 1
                        sm_response_zone3.dg = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height-start,density,abs(start-stop),resolution,0); % matrix output
                        dg_sm_1perc(i,3) = sum(sum(sm_response_zone3.dg));   % sum all pixels
                    else
                        sm_response_zone3.dg(:,:,end+1) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height-start,density,abs(start-stop),resolution,0); % matrix output
                        dg_sm_1perc(i,3) = sum(sum(sm_response_zone3.dg(:,:,end)));   % sum all pixels
                    end
                    sm_response_zone3.start = vertcat(sm_response_zone3.start,start);
                    sm_response_zone3.stop = vertcat(sm_response_zone3.stop,stop);
                end
            else
                dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180); % compute spherical distance
                dist = dist.*180/pi;                                        % transform to degrees                     
                density(dist<=radius{2}) = 0;                               % remove points/pixels out of zone 3 (set density to zero)
                density(dist>radius{3}) = 0;clear dist
                if output_file{5} == 0 
                    dg_sm_1perc(i,3) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height-stop,abs(start-stop),density,resolution(1),resol3lat,6371008,1);
                else
                    if size(sm_response_zone3.dg,1) == 1
                        sm_response_zone3.dg = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height-stop,abs(start-stop),density,resolution(1),resol3lat,6371008,0);
                        dg_sm_1perc(i,3) = sum(sum(sm_response_zone3.dg));
                    else
                        sm_response_zone3.dg(:,:,end+1) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height-stop,abs(start-stop),density,resolution(1),resol3lat,6371008,0);
                        dg_sm_1perc(i,3) = sum(sum(sm_response_zone3.dg(:,:,end)));
                    end
                    sm_response_zone3.start = vertcat(sm_response_zone3.start,start);
                    sm_response_zone3.stop = vertcat(sm_response_zone3.stop,stop);
                end
            end
            clear density
        end
        
        % Output table (SM)
        if ~isempty(sm{i,5})                                                % if input file is selected
            data_sm = mLocal_load1d(sm{i,5},1,sm{i,6});                     % Load data. Time vector is always (except for TSF) in the first row
            if ~exist('ref_time','var')
                ref_time = data_sm(:,1);                                    % crate referenece time (from first loaded file)
                out.sm_input(1:length(ref_time),i) = data_sm(:,2);          % Input SM (Vol%). Each column represent different SM layer
            elseif length(data_sm(:,1)) == length(ref_time)                 % Interpolate only if the other file has different time sampling
               if sum(data_sm(:,1)-ref_time)==0                             % no interpolation required as both time series have the same time sampling
                   out.sm_input(1:length(ref_time),i) = data_sm(:,2);
               else                                                         % interpolation required as the time series have the same number of data points but different sampling points/epochs
                   out.sm_input(1:length(ref_time),i) = interp1(data_sm(:,1),data_sm(:,2),ref_time);
               end
            else                                                            % interpolation required as the time series have different number of data points/epochs
                out.sm_input(1:length(ref_time),i) = interp1(data_sm(:,1),data_sm(:,2),ref_time); 
            end
            out.sm_zone1_dg(1:length(ref_time),i) = dg_sm_1perc(i,1).*out.sm_input(1:length(ref_time),i)*1e+9; % final SM gravity effect (nm/s^2) for zone1
            out.sm_zone2_dg(1:length(ref_time),i) = dg_sm_1perc(i,2).*out.sm_input(1:length(ref_time),i)*1e+9; % final SM gravity effect (nm/s^2) for zone2
            out.sm_zone3_dg(1:length(ref_time),i) = dg_sm_1perc(i,3).*out.sm_input(1:length(ref_time),i)*1e+9; % final SM gravity effect (nm/s^2) for zone3
            % Create GW and Snow variables (fill with zeros)
            out.gw_zone1_dg(1:length(ref_time),1) = 0;                      % prepare these variable only, they will be filled later
            out.gw_zone2_dg(1:length(ref_time),1) = 0;
            out.gw_zone3_dg(1:length(ref_time),1) = 0;
            out.snow_zone1_dg(1:length(ref_time),1) = 0;
            out.snow_zone2_dg(1:length(ref_time),1) = 0;
            out.snow_zone3_dg(1:length(ref_time),1) = 0;
        else
            out.sm_input(1,i) = 1;                                          % If no file is selected, all computation for 1 Vol%
            out.sm_zone1_dg(1,i) = dg_sm_1perc(i,1)*1e+9;                   % nm/s^2
            out.sm_zone2_dg(1,i) = dg_sm_1perc(i,2)*1e+9;
            out.sm_zone3_dg(1,i) = dg_sm_1perc(i,3)*1e+9;
            % Prepare other variables for writing
            out.snow_zone1_dg=0;out.snow_zone2_dg = 0;out.snow_zone3_dg = 0;
            out.gw_zone1_dg=0;out.gw_zone2_dg = 0;out.gw_zone3_dg = 0;
        end
    end
end
clear i;                                                                    % clear the loop index

% Save matrices + delete
if output_file{5} == 1 && sm_exist >= 1                                     % if response matrices are required
    out_name = output_file{1};                                              % output file name
    switch exist('sm_response_zone1','var')+exist('sm_response_zone2','var')+exist('sm_response_zone3','var') % switch between number of computed zones
        case 1
            save([out_name(1:end-3),'mat'],'sm_response_zone1','InputUse'); % first zone only (this will overwrite the existing MAT file) + InputUse (shperical coor.)
            clear sm_response_zone1                                         % clear saved variable
        case 2
            save([out_name(1:end-3),'mat'],'sm_response_zone1','sm_response_zone2','InputUse'); % first and second zone
            clear sm_response_zone1 sm_response_zone2
        case 3
            save([out_name(1:end-3),'mat'],'sm_response_zone1','sm_response_zone2','sm_response_zone3','InputUse'); % all zones
            clear sm_response_zone1 sm_response_zone2 sm_response3
    end
end

%% Groundwater
if gw{1}==1 && ~isempty(gw{6})                                              % compute only if file and computation checkbox are selected
    set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: groundwater...'));drawnow; % send message to status bar
    data_gw = mLocal_load1d(gw{6},1,gw{7});                                 % Load GW data (one dataset for all zones)
    if ~exist('ref_time','var')
        ref_time = data_gw(:,1);                                            % Create ref.time if no SM has been selected
        out.gw_input = data_gw(:,2);                                        % use loaded GW data as input (no interpolation)
    elseif length(data_gw(:,1)) == length(ref_time)                         % Interpolate only if the other file has different time sampling
       if sum(data_gw(:,1)-ref_time)==0                                     % => same time
           out.gw_input = data_gw(:,2);                                     % no interpolation as both time series have the same time vector
       else
           out.gw_input = interp1(data_gw(:,1),data_gw(:,2),ref_time);      % time vectors are not equal => interpolate
       end
    else
        out.gw_input = interp1(data_gw(:,1),data_gw(:,2),ref_time);         % interpolate gw data
    end
    
    if length(ref_time) > 20                                                % Select between computing dg for each time epoch (<20) or interpolate from pre-calculated range (to spare time)
        ref_gw = linspace(min(out.gw_input),max(out.gw_input),21);          % create new range that will be used for interpolation
    else
        ref_gw = out.gw_input;                                              % use the loaded time series
    end                                                                     % from now on, time of GW observation is not important any more
    
    dg_gw(1:length(ref_gw),1:3) = 0;                                        % Prepare dg variable
    if output_file{5} == 1                                                  % Create matrices for gravity response
        if length(ref_gw) > 2
            gw_response_mat = 2;                                            % Max. 2 matrices, i.e. for min. and max. observed GW
        else
            gw_response_mat = length(ref_gw);                               % Ohterwise, only one matrix
        end
        gw_response_zone1.x = dem1.x;gw_response_zone1.y = dem1.y;
        gw_response_zone1.dg(1:size(dem1.x,1),1:size(dem1.x,2),1:gw_response_mat) = 0;
        gw_response_zone1.input_units = 'm';gw_response_zone1.dg_units = 'm/s^2';
        gw_response_zone1.specYield = gw{3};
        if coor_system == 1
            gw_response_zone1.dist = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2);
            gw_response_zone1.dist_units = 'm';
        else
            gw_response_zone1.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi; % compute spherical distance
            gw_response_zone1.ySphere = dem1.ySphere;
            gw_response_zone1.dist_units = 'deg';
        end
        if ~isempty(dem{2})
            gw_response_zone2.x = dem2.x;gw_response_zone2.y = dem2.y;
            gw_response_zone2.dg(1:size(dem2.x,1),1:size(dem2.x,2),1:gw_response_mat) = 0;
            gw_response_zone2.input_units = 'm';gw_response_zone2.dg_units = 'm/s^2';
            gw_response_zone2.specYield = gw{3};
            if coor_system == 1
                gw_response_zone2.dist = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2);
                gw_response_zone2.dist_units = 'm';
            else
                gw_response_zone2.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi; % compute spherical distance
                gw_response_zone2.ySphere = dem2.ySphere;
                gw_response_zone2.dist_units = 'deg';
            end
        end
        if ~isempty(dem{3})
            gw_response_zone3.x = dem3.x;gw_response_zone3.y = dem3.y;
            gw_response_zone3.dg(1:size(dem3.x,1),1:size(dem3.x,2),1:gw_response_mat) = 0;
            gw_response_zone3.input_units = 'm';gw_response_zone3.dg_units = 'm/s^2';
            gw_response_zone3.specYield = gw{3};
            if coor_system == 1
                gw_response_zone3.dist = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2);
                gw_response_zone3.dist_units = 'm';
            else
                gw_response_zone3.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi; % compute spherical distance
                gw_response_zone3.ySphere = dem3.ySphere;
                gw_response_zone3.dist_units = 'deg';
            end
        end
        clear gw_response_mat;
    end
    
    polyg_file = gw{2};                                                     % Selected exclusion polygon
    if ~isempty(polyg_file)                                                 % Load Exclusion polygon (same for all zones)
        polyg_switcher = polyg_file(end-2:end);
        switch polyg_switcher
            case 'txt'
                polyg = load(polyg_file);                                   % Load TXT
        end
    end
    
    for i = 1:length(ref_gw)                                                % Loop for each reference GW (max 21)
        % ZONE 1
        if gw{4} == 2
            min_gw = dem1.height-max(ref_gw)-1;                             % Lower boundary = fixed in time (depth input)
        else
            min_gw = dem1.height - gw{5} - 1;                               % Lower boundary = fixed in time (altitude input) 
            min_gw(min_gw>min(ref_gw)-1) = min(ref_gw) - 1;                 % set values greater than min. value
        end
        density = dem1.height;density(:,:) = gw{3}*10;                      % Denity in kg/m^3 (specific yield in Vol%*10)
        density(dem1.height>=9e+35) = 0;                                    % Set former NaNs to zero
        resolution = resol1;                                                % spatial resolution of the first zone
        if gw{4} == 2                                                       % Copy terrain in distance given by GW depth
            max_gw = dem1.height-ref_gw(i);                                 % ref_gw == GW Depth (not altitude)!
        else                                                                % Copy terrain in constant distance or max altitude given by model
            max_gw = dem1.height - gw{5};                                   % ref_gw == GW Altitude (not depth)!
            max_gw(max_gw>ref_gw(i)) = ref_gw(i);
        end
        if ~isempty(polyg_file)                                             % use exclusion polygon if selected
            switch polyg_switcher
                case 'shp'
                    for p = 1:length(polyg)
                        id = inpolygon(dem1.x,dem1.y,polyg(p,1).X,polyg(p,1).Y);% find grid cells for each polygon in SHP file
                        density(id) = 0;    	                            % set grid cells 
                    end
                case 'txt'
                    if numel(polyg) == 1
                        switch coor_system
                            case 1
                                id = find(sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2)<=polyg);
                            case 2
                                id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi<=polyg);
                        end
                    else
                        id = inpolygon(dem1.x,dem1.y,polyg(:,1),polyg(:,2));
                    end
                    density(id) = 0;                                        % Exclude grid cells within loaded polygon
            end
            clear id p
        end
        
        if coor_system == 1                                                 % For cartesian coordinates
            dist = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2);   % temp. (plane) distance between pixels and Gravimeter
            density(dist>radius{1}) = 0;clear dist;                         % remove points/pixels out of zone 1 (set density to zero)
            if output_file{5} == 0 	                                        % if matrix/scalar output is required
                dg_gw(i,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
            else 
                if ref_gw(i) == min(ref_gw)                                 % Compute matrix only for min. and max. GW!
                    gw_response_zone1.dg(:,:,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                    gw_response_zone1.input(1) = ref_gw(i);
                    dg_gw(i,1) = sum(sum(gw_response_zone1.dg(:,:,1)));     % sum all pixels
                elseif ref_gw(i) == max(ref_gw)                             % For max. GW value
                    gw_response_zone1.dg(:,:,2) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                    gw_response_zone1.input(2) = ref_gw(i);
                    dg_gw(i,1) = sum(sum(gw_response_zone1.dg(:,:,2)));     % sum all pixels
                else                                                        % Do not compute the matrix for all GW values in between min. and max.
                    dg_gw(i,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
                end
            end
        else                                                                % For geodetic coordinates
            dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180); % compute spherical distance
            dist = dist.*180/pi;                                            % transform to degrees       
            density(dist>radius{1}) = 0;clear dist;                         % remove points/pixels out of zone 1 (zet density to zero)
            if output_file{5} == 0                                          % if matrix/scalar output is required
                dg_gw(i,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol1lat,6371008,1);
            else 
                if ref_gw(i) == min(ref_gw)                                 % for minimum GW
                    gw_response_zone1.dg(:,:,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol1lat,6371008,0);
                    gw_response_zone1.input(1) = ref_gw(i);
                    dg_gw(i,1) = sum(sum(gw_response_zone1.dg(:,:,1)));       % sum all pixels
                elseif ref_gw(i) == max(ref_gw)                             % for maximum GW
                    gw_response_zone1.dg(:,:,2) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol1lat,6371008,0);
                    gw_response_zone1.input(2) = ref_gw(i);
                    dg_gw(i,1) = sum(sum(gw_response_zone1.dg(:,:,2)));       % sum all pixels
                else                                                        % for all other cases (same as output_file{5} == 0)
                    dg_gw(i,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol1lat,6371008,1);
                end
            end
        end
        clear density max_gw min_gw                                         % delete used variables
        
        % ZONE 2
        if ~isempty(dem{2})                                                 % only if zone 2 is selected
            if gw{4} == 2
                min_gw = dem2.height-max(ref_gw)-1;                         % Lower boundary = fixed in time (depth input)
            else
                min_gw = dem2.height - gw{5} - 1;                           % Lower boundary = fixed in time (altitude input) 
                min_gw(min_gw>min(ref_gw)-1) = min(ref_gw) - 1;
            end
            density = dem2.height;density(:,:) = gw{3}*10;                  % Denity in kg/m^3 (specific yield in Vol%*10)
            density(dem2.height>=9e+35) = 0;                                % Set former NaNs to zero
            resolution = resol2;                                            % spatial resolution of the second zone
            if gw{4} == 2                                                   % Copy terrain in distance given by GW depth
                max_gw = dem2.height-ref_gw(i);                             % ref_gw == GW Depth (not altitude)!
            else                                                            % Copy terrain in constant distance or max altitude given by model
                max_gw = dem2.height - gw{5};                               % ref_gw == GW Altitude (not depth)!
                max_gw(max_gw>ref_gw(i)) = ref_gw(i);                       % Set values above the max. value to current GW level
            end
            if ~isempty(polyg_file)                                         % use exclusion polygon if selected
                switch polyg_switcher
                    case 'shp'
                        for p = 1:length(polyg)
                            id = inpolygon(dem2.x,dem2.y,polyg(p,1).X,polyg(p,1).Y); % find grid cells for each polygon in SHP file
                            density(id) = 0;	                            % set grid cells 
                        end
                    case 'txt'
                        if numel(polyg) == 1
                            switch coor_system
                                case 1
                                    id = find(sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2)<=polyg);
                                case 2
                                    id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi<=polyg);
                            end
                        else
                            id = inpolygon(dem2.x,dem2.y,polyg(:,1),polyg(:,2));
                        end
                        density(id) = 0;	                                % Exclude grid cells within loaded polygon
                end
                clear id p
            end

            if coor_system == 1                                             % For cartesian coordinates
                dist = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2); % temp. (plane) distance between pixels and Gravimeter
                density(dist<=radius{1}) = 0;                               % remove points/pixels out of zone 2 (zet density to zero)
                density(dist>radius{2}) = 0;clear dist
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_gw(i,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
                else 
                    if ref_gw(i) == min(ref_gw)                             % Compute matrix only for min. and max. GW!
                        gw_response_zone2.dg(:,:,1) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                        gw_response_zone2.input(1) = ref_gw(i);
                        dg_gw(i,2) = sum(sum(gw_response_zone2.dg(:,:,1))); % sum all pixels
                    elseif ref_gw(i) == max(ref_gw)                         % For max. GW value
                        gw_response_zone2.dg(:,:,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                        gw_response_zone2.input(2) = ref_gw(i);
                        dg_gw(i,2) = sum(sum(gw_response_zone2.dg(:,:,2))); % sum all pixels
                    else                                                    % Do not compute the matrix for all GW values in between min. and max.
                        dg_gw(i,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
                    end
                end
            else                                                            % For geodetic coordinates
                dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180); % compute spherical distance
                dist = dist.*180/pi;                                        % transform to degrees       
                density(dist<=radius{1}) = 0;                               % remove points/pixels out of zone 2 (zet density to zero)
                density(dist>radius{2}) = 0;clear dist
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_gw(i,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol2lat,6371008,1);
                else 
                    if ref_gw(i) == min(ref_gw)
                        gw_response_zone2.dg(:,:,1) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol2lat,6371008,0);
                        gw_response_zone2.input(1) = ref_gw(i);
                        dg_gw(i,2) = sum(sum(gw_response_zone2.dg(:,:,1))); % sum all pixels
                    elseif ref_gw(i) == max(ref_gw)
                        gw_response_zone2.dg(:,:,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol2lat,6371008,0);
                        gw_response_zone2.input(2) = ref_gw(i);
                        dg_gw(i,2) = sum(sum(gw_response_zone2.dg(:,:,2))); % sum all pixels
                    else
                        dg_gw(i,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol2lat,6371008,1);
                    end
                end
            end
            clear density max_gw min_gw
        end
        
        % ZONE 3
        if ~isempty(dem{3})
            if gw{4} == 2
                min_gw(:,:) = dem3.height-max(ref_gw)-1;                    % Lower boundary = fixed in time (depth input)
            else
                min_gw = dem3.height - gw{5} - 1;                           % Lower boundary = fixed in time (altitude input) 
                min_gw(min_gw>min(ref_gw)-1) = min(ref_gw) - 1;
            end
            density = dem3.height;density(:,:) = gw{3}*10;                  % Denity in kg/m^3 (specific yield in Vol%*10)
            density(dem3.height>=9e+35) = 0;                                % Set former NaNs to zero
            resolution = resol3;                                            % spatial resolution of the third zone
            if gw{4} == 2                                                   % Copy terrain in distance given by GW depth
                max_gw = dem3.height-ref_gw(i);                             % ref_gw == GW Depth (not altitude)!
            else                                                            % Copy terrain in constant distance or max altitude given by model
                max_gw = dem3.height - gw{5};                               % ref_gw == GW Altitude (not depth)!
                max_gw(max_gw>ref_gw(i)) = ref_gw(i);
            end
            if ~isempty(polyg_file)                                         % use exclusion polygon if selected
                switch polyg_switcher
                    case 'shp'
                        for p = 1:length(polyg)
                            id = inpolygon(dem3.x,dem3.y,polyg(p,1).X,polyg(p,1).Y);% find grid cells for each polygon in SHP file
                            density(id) = 0;                                % set grid cells 
                        end
                    case 'txt'
                        if numel(polyg) == 1
                            switch coor_system
                                case 1
                                    id = find(sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2)<=polyg);
                                case 2
                                    id = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi<=polyg);
                            end
                        else
                            id = inpolygon(dem3.x,dem3.y,polyg(:,1),polyg(:,2));
                        end
                        density(id) = 0;                                   % Exclude grid cells within loaded polygon
                end
                clear id p
            end

            if coor_system == 1                                             % For cartesian coordinates
                dist = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2); % temp. (plane) distance between pixels and Gravimeter
                density(dist<=radius{2}) = 0;                               % remove points/pixels out of zone 3 (set density to zero)
                density(dist>radius{3}) = 0;clear dist
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_gw(i,3) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
                else 
                    if ref_gw(i) == min(ref_gw)                                               % Compute matrix only for min. and max. GW!
                        gw_response_zone3.dg(:,:,1) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                        gw_response_zone3.input(1) = ref_gw(i);
                        dg_gw(i,3) = sum(sum(gw_response_zone3.dg(:,:,1))); % sum all pixels
                    elseif ref_gw(i) == max(ref_gw)                         % For max. GW value
                        gw_response_zone3.dg(:,:,2) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,max_gw,density,max_gw-min_gw,resolution,0); % matrix output
                        gw_response_zone3.input(2) = ref_gw(i);
                        dg_gw(i,3) = sum(sum(gw_response_zone3.dg(:,:,2))); % sum all pixels
                    else                                                    % Do not compute the matrix for all GW values in between min. and max.
                        dg_gw(i,3) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,max_gw,density,max_gw-min_gw,resolution,1); % sum of all pixels
                    end
                end
            else                                                            % For geodetic coordinates
                dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180); % compute spherical distance
                dist = dist.*180/pi;                                        % transform to degrees       
                density(dist<=radius{2}) = 0;                               % remove points/pixels out of zone 3 (set density to zero)
                density(dist>radius{3}) = 0;clear dist
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_gw(i,3) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol3lat,6371008,1);
                else 
                    if ref_gw(i) == min(ref_gw)
                        gw_response_zone3.dg(:,:,1) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol3lat,6371008,0);
                        gw_response_zone3.input(1) = ref_gw(i);
                        dg_gw(i,3) = sum(sum(gw_response_zone3.dg(:,:,1))); % sum all pixels
                    elseif ref_gw(i) == max(ref_gw)
                        gw_response_zone3.dg(:,:,2) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol3lat,6371008,0);
                        gw_response_zone3.input(2) = ref_gw(i);
                        dg_gw(i,3) = sum(sum(gw_response_zone3.dg(:,:,2))); % sum all pixels
                    else
                        dg_gw(i,3) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,min_gw,max_gw-min_gw,density,resolution(1),resol3lat,6371008,1);
                    end
                end
            end
            clear density max_gw min_gw 
        end
        set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: groundwater (%2.0f%%)...',round(i/length(ref_gw)*100)));drawnow; % send message to status bar
    end
    
    % Output table (Groundwater)
    if length(ref_time) > 20
        out.gw_zone1_dg(:,1) = interp1(ref_gw,dg_gw(:,1),out.gw_input)*1e+9; % interpolate gravity effect for zone 1
        out.gw_zone2_dg(:,1) = interp1(ref_gw,dg_gw(:,2),out.gw_input)*1e+9;
        out.gw_zone3_dg(:,1) = interp1(ref_gw,dg_gw(:,3),out.gw_input)*1e+9;
    else
        out.gw_zone1_dg = dg_gw(:,1)*1e+9;                                  % equals the computed gravity effect transpormed to nm/s^2
        out.gw_zone2_dg = dg_gw(:,2)*1e+9;
        out.gw_zone3_dg = dg_gw(:,3)*1e+9;
    end
    % Prepare Snow variable (filled with zeros)
    out.snow_zone1_dg(1:length(ref_time),1) = 0;
    out.snow_zone2_dg(1:length(ref_time),1) = 0;
    out.snow_zone3_dg(1:length(ref_time),1) = 0;
    % Prepare SM variable for writing
    if sm_file_exist == 0                                                   % fill with zeros only if not already computed
        out.sm_zone1_dg(1:length(ref_time),1) = 0;
        out.sm_zone2_dg(1:length(ref_time),1) = 0;
        out.sm_zone3_dg(1:length(ref_time),1) = 0;
    end
        
    clear ref_gw dg_gw polyg i resolution
    
    % Save matrices + delete
    if output_file{5} == 1
        out_name = output_file{1};
        if ~exist([out_name(1:end-3),'mat'],'file');                        % check if the file already exists = if the SM has been computed. Warning this also means that in will write to an existin file regardless of file origin!!
            switch exist('gw_response_zone1','var')+exist('gw_response_zone2','var')+exist('gw_response_zone3','var') % check how many zones have been computed
                case 1
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','InputUse'); % One zone (new file = overwrite)
                    clear gw_response_zone1
                case 2
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','gw_response_zone2','InputUse'); % two zones (new file)
                    clear gw_response_zone1 gw_response_zone2
                case 3
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','gw_response_zone2','gw_response_zone3','InputUse'); % three zones (new file)
                    clear gw_response_zone1 gw_response_zone2 gw_response3
            end
        else
            switch exist('gw_response_zone1','var')+exist('gw_response_zone2','var')+exist('gw_response_zone3','var')
                case 1
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','-append'); % One zone (append variables)
                    clear gw_response_zone1
                case 2
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','gw_response_zone2','-append'); % two zones (append variables)
                    clear gw_response_zone1 gw_response_zone2
                case 3
                    save([out_name(1:end-3),'mat'],'gw_response_zone1','gw_response_zone2','gw_response_zone3','-append'); % 3 zones (append variables)
                    clear gw_response_zone1 gw_response_zone2 gw_response3
            end
        end
    end
end


%% Snow
if snow{1}==1 && ~isempty(snow{6})
    set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: snow...'));drawnow; % send message to status bar
    if snow{8}<=0
        data_snow = mLocal_load1d(snow{6},1,snow{7});                       % Load Snow data = only snow depth (one dataset for all zones)
    else
        data_snow = mLocal_load1d(snow{6},1,[snow{7},snow{8}]);             % Load Snow data = snow depth + snow water equivalent (SWE), both in mm
    end
    if ~exist('ref_time','var')
        ref_time = data_snow(:,1);                                          % Create ref.time if no SM or GW has been selected
        out.snow_input = data_snow(:,2:end);                                % use loaded Snow data as input (no interpolation)
    elseif length(data_snow(:,1)) == length(ref_time)                       % Interpolate only if the other file has different time sampling
       if sum(data_snow(:,1)-ref_time)==0                                   % => same time
           out.snow_input = data_snow(:,2:end);                             % no interpolation as both time series have the same time vector
       else
           for j = 2:size(data_snow,2)                                      % For both, height and SWE
                out.snow_input(1:length(ref_time),j-1) = interp1(data_snow(:,1),data_snow(:,j),ref_time); % time vectors are not equal => interpolate
           end
           clear j;
       end
    else
       for j = 2:size(data_snow,2)                                          % For both, height and SWE
            out.snow_input(1:length(ref_time),j-1) = interp1(data_snow(:,1),data_snow(:,j),ref_time); % time vectors are not equal => interpolate
       end
       clear j;
    end
    
    dg_snow(1:length(ref_time),1:3) = 0;                                    % Prepare dg variable (function of time)
    if output_file{5} == 1                                                  % Create matrix for gravity response (only one, i.e. for max snow height)
        snow_response_zone1.x = dem1.x;snow_response_zone1.y = dem1.y;
        snow_response_zone1.dg(1:size(dem1.x,1),1:size(dem1.x,2)) = 0;      % Only 2D matrix = for max. snow height (not max. SWE!)
        snow_response_zone1.input_units = 'mm';snow_response_zone1.dg_units = 'm/s^2';
        if coor_system == 1
            snow_response_zone1.dist = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2);
            snow_response_zone1.dist_units = 'm';
        else
            snow_response_zone1.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi; % compute spherical distance
            snow_response_zone1.ySphere = dem1.ySphere;
            snow_response_zone1.dist_units = 'deg';
        end
        if ~isempty(dem{2})                                                 % Second zone
            snow_response_zone2.x = dem2.x;snow_response_zone2.y = dem2.y;
            snow_response_zone2.dg(1:size(dem2.x,1),1:size(dem2.x,2)) = 0;
            snow_response_zone2.input_units = 'mm';snow_response_zone2.dg_units = 'm/s^2';
            if coor_system == 1
                snow_response_zone2.dist = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2);
                snow_response_zone2.dist_units = 'm';
            else
                snow_response_zone2.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi; % compute spherical distance
                snow_response_zone2.ySphere = dem2.ySphere;
                snow_response_zone2.dist_units = 'deg';
            end
        end
        if ~isempty(dem{3})                                                 % Third zone
            snow_response_zone3.x = dem3.x;snow_response_zone3.y = dem3.y;
            snow_response_zone3.dg(1:size(dem3.x,1),1:size(dem3.x,2)) = 0;
            snow_response_zone3.input_units = 'mm';snow_response_zone3.dg_units = 'm/s^2';
            if coor_system == 1
                snow_response_zone3.dist = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2);
                snow_response_zone3.dist_units = 'm';
            else
                snow_response_zone3.dist = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi; % compute spherical distance
                snow_response_zone3.ySphere = dem3.ySphere;
                snow_response_zone3.dist_units = 'deg';
            end
        end
    end
    
    polyg_file = snow{2};                                                   % Selected exclusion polygon
    if ~isempty(polyg_file)                                                 % Load Exclusion polygon (same for all zones)
        polyg_switcher = polyg_file(end-2:end);                             % file extension
        switch polyg_switcher
            case 'txt'
                polyg = load(polyg_file);
        end
    end
    
    if ~isempty(polyg_file)                                         % use exclusion polygon if selected
        switch polyg_switcher
            case 'shp'
                id1 = dem1.height.*0;
                for p = 1:length(polyg)
                    id1 = id1 + double(inpolygon(dem1.x,dem1.y,polyg(p,1).X,polyg(p,1).Y));% find grid cells for each polygon in SHP file
                end
                id1 = logical(id1);
            case 'txt'
                if numel(polyg) == 1
                    switch coor_system
                        case 1
                            id1 = find(sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2)<=polyg);
                        case 2
                            id1 = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180)*180/pi<=polyg);
                    end
                else
                    id1 = inpolygon(dem1.x,dem1.y,polyg(:,1),polyg(:,2));
                end
        end
        clear p
    end
    
    if coor_system == 1 
        dist1 = sqrt((dem1.x-InputUse(1)).^2+(dem1.y-InputUse(2)).^2); % temp. (plane) distance between pixels and Gravimeter
    else
        dist1 = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem1.x.*pi/180,dem1.y.*pi/180); % compute spherical distance
        dist1 = dist1.*180/pi;                                        % transform to degrees  
    end
    
    if ~isempty(dem{2})
        switch polyg_switcher
            case 'shp'
                id2 = dem2.height.*0;
                for p = 1:length(polyg)
                    id2 = id2 + double(inpolygon(dem2.x,dem2.y,polyg(p,1).X,polyg(p,1).Y));% find grid cells for each polygon in SHP file
                end
                id2 = logical(id2);
            case 'txt'
                if numel(polyg) == 1
                switch coor_system
                    case 1
                        id2 = find(sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2)<=polyg);
                    case 2
                        id2 = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180)*180/pi<=polyg);
                end
                else
                    id2 = inpolygon(dem2.x,dem2.y,polyg(:,1),polyg(:,2));
                end
        end
        if coor_system == 1 
            dist2 = sqrt((dem2.x-InputUse(1)).^2+(dem2.y-InputUse(2)).^2); % temp. (plane) distance between pixels and Gravimeter
        else
            dist2 = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem2.x.*pi/180,dem2.y.*pi/180); % compute spherical distance
            dist2 = dist2.*180/pi;                                        % transform to degrees  
        end
    end
    if ~isempty(dem{3})
        if ~isempty(polyg_file)                                     % use exclusion polygon if selected
            switch polyg_switcher
                case 'shp'
                    id3 = dem3.height.*0;
                    for p = 1:length(polyg)
                        id3 = id3 + double(inpolygon(dem3.x,dem3.y,polyg(p,1).X,polyg(p,1).Y));% find grid cells for each polygon in SHP file
                    end
                    id3 = logical(id3);
                case 'txt'
                    if numel(polyg) == 1
                    switch coor_system
                        case 1
                            id3 = find(sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2)<=polyg);
                        case 2
                            id3 = find(LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180)*180/pi<=polyg);
                    end
                    else
                        id3 = inpolygon(dem3.x,dem3.y,polyg(:,1),polyg(:,2));
                    end
            end
            clear p
        end
        if coor_system == 1
            dist3 = sqrt((dem3.x-InputUse(1)).^2+(dem3.y-InputUse(2)).^2); % temp. (plane) distance between pixels and Gravimeter
        else
            dist3 = LonLat2psi(Input(1)*pi/180,Input(2)*pi/180,dem3.x.*pi/180,dem3.y.*pi/180); % compute spherical distance
            dist3 = dist3.*180/pi;                                    % transform to degrees   
        end
    end
    
            
    for i = 1:length(ref_time)
        if out.snow_input(i,1) > 0                                         % Compute only if snow height ~= zero
            % ZONE 1
            resolution = resol1;
            density = dem1.height.*0;                                          % Prepare variable (get final dimensions)
            if snow{5} > 0                                                  % If snow gradient is not equal to zero
                temp_height = out.snow_input(i,1)/1000 - (Input(3)-dem1.height)/100*snow{5}/1000; % new snow height matrix (m). 
                temp_height(temp_height<0) = 0;                             % remove possible negative values (the resulting Inf will be removed later)
            else
                temp_height = out.snow_input(i,1)/1000;                     % Constant snow height (m)
            end
            if snow{8} <= 0                                                 % Compute snow density if SWE is not loaded
                density(:,:) = snow{3} + snow{4}*temp_height;
            else                                                            % Otherwise use SWE for the compuation of snow density
                density(:,:) = out.snow_input(i,2)./temp_height;            % kg/m^3 (mm/m => kg/m^3)
            end
            density(dem1.height>=9e+35) = 0;                                % set former NaNs to zero
            temp_height(isinf(temp_height)) = 0;                            % remove posible Inf values
            density(temp_height<=0) = 0;                                    % remove posible Inf values
            density(id1) = 0;                                                   % set grid cells
            
            if coor_system == 1                                             % For cartesian coordinates
                density(dist1>radius{1}) = 0;                                   % remove points/pixels out of zone 1 (set density to zero)
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_snow(i,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                else 
                    if out.snow_input(i,1) == max(out.snow_input(:,1))      % Compute matrix only for max. snow height
                        snow_response_zone1.dg = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height+temp_height,density,temp_height,resolution,0); % matrix output
                        snow_response_zone1.input = out.snow_input(i,:);
                        dg_snow(i,1) = sum(sum(snow_response_zone1.dg));    % sum all pixels
                    else
                        dg_snow(i,1) = mLocal_prismEffect(InputUse,dem1.x,dem1.y,dem1.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                    end
                end
            else                                                            % For geodetic coordinates
                density(dist1>radius{1}) = 0;                               % remove points/pixels out of zone 1 (set density to zero)
                if output_file{5} == 0                                      % if matrix/scalar output is required
                    dg_snow(i,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height,temp_height,density,resolution(1),resol1lat,6371008,1);
                else 
                    if out.snow_input(i,1) == max(out.snow_input(:,1))
                        snow_response_zone1.dg = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height,temp_height,density,resolution(1),resol1lat,6371008,0);
                        snow_response_zone1.input = out.snow_input(i,:);
                        dg_snow(i,1) = sum(sum(snow_response_zone1.dg));    % sum all pixels
                    else
                        dg_snow(i,1) = mLocal_tesseroid(InputUse,dem1.x,dem1.ySphere,dem1.height,temp_height,density,resolution(1),resol1lat,6371008,1);
                    end
                end
            end
           clear density resolution temp_height                             % remove used variables 
           
           % ZONE 2
           if ~isempty(dem{2})
                resolution = resol2;
                density = dem2.height.*0;                                      % Prepare variable (get final dimensions)
                if snow{5} > 0                                              % If snow gradient is not equal to zero
                    temp_height = out.snow_input(i,1)/1000 - (Input(3)-dem2.height)/100*snow{5}/1000; % new snow height matrix (m). 
                    temp_height(temp_height<0) = 0;                         % remove possible negative values (the resulting Inf will be removed later)
                else
                    temp_height = out.snow_input(i,1)/1000;                 % Constant snow height (m)
                end
                if snow{8} <= 0                                             % Compute snow density if SWE is not loaded
                    density(:,:) = snow{3} + snow{4}*temp_height;
                else                                                        % Otherwise use SWE for the compuation of snow density
                    density(:,:) = out.snow_input(i,2)./temp_height;        % kg/m^3 
                end
                density(dem2.height>=9e+35) = 0;                            % set former NaNs to zero
                temp_height(isinf(temp_height)) = 0;                        % remove posible Inf values
                density(temp_height<=0) = 0;                                % remove posible Inf values

                density(id2) = 0;                            % set grid cells 
                if coor_system == 1                                         % For cartesian coordinates
                    density(dist2<=radius{1}) = 0;                           % remove points/pixels out of zone 2 (zet density to zero)
                    density(dist2>radius{2}) = 0;
                    if output_file{5} == 0                                  % if matrix/scalar output is required
                        dg_snow(i,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                    else 
                        if out.snow_input(i,1) == max(out.snow_input(:,1))  % Compute matrix only for max. snow height
                            snow_response_zone2.dg = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height+temp_height,density,temp_height,resolution,0); % matrix output
                            snow_response_zone2.input = out.snow_input(i,:);
                            dg_snow(i,2) = sum(sum(snow_response_zone2.dg)); % sum all pixels
                        else
                            dg_snow(i,2) = mLocal_prismEffect(InputUse,dem2.x,dem2.y,dem2.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                        end
                    end
                else                                                        % For geodetic coordinates    
                    density(dist2<=radius{1}) = 0;                           % remove points/pixels out of zone 2 (zet density to zero)
                    density(dist2>radius{2}) = 0;
                    if output_file{5} == 0                                  % if matrix/scalar output is required
                        dg_snow(i,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height,temp_height,density,resolution(1),resol2lat,6371008,1);
                    else 
                        if out.snow_input(i,1) == max(out.snow_input(:,1))
                            snow_response_zone2.dg = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height,temp_height,density,resolution(1),resol2lat,6371008,0);
                            snow_response_zone2.input = out.snow_input(i,:);
                            dg_snow(i,2) = sum(sum(snow_response_zone2.dg));% sum all pixels
                        else
                            dg_snow(i,2) = mLocal_tesseroid(InputUse,dem2.x,dem2.ySphere,dem2.height,temp_height,density,resolution(1),resol2lat,6371008,1);
                        end
                    end
                end
                clear density resolution temp_height                        % remove used variables 
           end
           
           % ZONE 3
           if ~isempty(dem{3})
                resolution = resol3;
                density = dem3.height.*0;                                      % Prepare variable (get final dimensions)
                if snow{5} > 0                                              % If snow gradient is not equal to zero
                    temp_height = out.snow_input(i,1)/1000 - (Input(3)-dem3.height)/100*snow{5}/1000; % new snow height matrix (m). 
                    temp_height(temp_height<0) = 0;                         % remove possible negative values (the resulting Inf will be removed later)
                else
                    temp_height = out.snow_input(i,1)/1000;                 % Constant snow height (m)
                end
                if snow{8} <= 0                                             % Compute snow density if SWE is not loaded
                    density(:,:) = snow{3} + snow{4}*temp_height;
                else                                                        % Otherwise use SWE for the compuation of snow density
                    density(:,:) = out.snow_input(i,2)./temp_height;        % kg/m^3 
                end
                density(dem3.height>=9e+35) = 0;                            % set former NaNs to zero
                temp_height(isinf(temp_height)) = 0;                        % remove posible Inf values
                density(temp_height<=0) = 0;                                % remove posible Inf values
                density(id3) = 0;                                           % set grid cells in polygon to zero
                
                if coor_system == 1                                         % For cartesian coordinates
                    density(dist3<=radius{2}) = 0;                           % remove points/pixels out of zone 3 (set density to zero)
                    density(dist3>radius{3}) = 0;
                    if output_file{5} == 0                                  % if matrix/scalar output is required
                        dg_snow(i,3) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                    else 
                        if out.snow_input(i,1) == max(out.snow_input(:,1))  % Compute matrix only for max. snow height
                            snow_response_zone3.dg = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height+temp_height,density,temp_height,resolution,0); % matrix output
                            snow_response_zone3.input = out.snow_input(i,:);
                            dg_snow(i,3) = sum(sum(snow_response_zone3.dg)); % sum all pixels
                        else
                            dg_snow(i,3) = mLocal_prismEffect(InputUse,dem3.x,dem3.y,dem3.height+temp_height,density,temp_height,resolution,1); % sum of all pixels
                        end
                    end
                else                                                        % For geodetic coordinates    
                    density(dist3<=radius{2}) = 0;                               % remove points/pixels out of zone 3 (set density to zero)
                    density(dist3>radius{3}) = 0;
                    if output_file{5} == 0                                  % if matrix/scalar output is required
                        dg_snow(i,3) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height,temp_height,density,resolution(1),resol3lat,6371008,1);
                    else 
                        if out.snow_input(i,1) == max(out.snow_input(:,1))
                            snow_response_zone3.dg = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height,temp_height,density,resolution(1),resol3lat,6371008,0);
                            snow_response_zone3.input = out.snow_input(i,:);
                            dg_snow(i,3) = sum(sum(snow_response_zone3.dg));% sum all pixels
                        else
                            dg_snow(i,3) = mLocal_tesseroid(InputUse,dem3.x,dem3.ySphere,dem3.height,temp_height,density,resolution(1),resol3lat,6371008,1);
                        end
                    end
                end
                clear density resolution temp_height                        % remove used variables 
           end
        else
            dg_snow(i,1:3) = 0;                                             % zero if the input snow height = 0
        end
        
        % Output table (Snow)
        out.snow_zone1_dg(i,1) = dg_snow(i,1)*1e+9;                         % to nm/s^2
        out.snow_zone2_dg(i,1) = dg_snow(i,2)*1e+9;
        out.snow_zone3_dg(i,1) = dg_snow(i,3)*1e+9;
        if gw{1}==0                                                         % Prepare variable for writing
            out.gw_zone1_dg(i,1) = 0;
            out.gw_zone2_dg(i,1) = 0;
            out.gw_zone3_dg(i,1) = 0;
        end
        % Prepare SM variable for writing
        if sm_file_exist == 0                                               % fill with zeros only if not already computed
            out.sm_zone1_dg(i,1) = 0;
            out.sm_zone2_dg(i,1) = 0;
            out.sm_zone3_dg(i,1) = 0;
        end
        set(findobj('Tag','text_status_mLocal'),'String',sprintf('Gravity response: snow (%2.0f%%)...',round(i/length(ref_time)*100)));drawnow; % send message to status bar
    end
    clear i id1 id2 id3 dist1 dist2 dist3;
    
    % Save matrices + delete
    if output_file{5} == 1
        out_name = output_file{1};
        if ~exist([out_name(1:end-3),'mat'],'file');
            switch exist('snow_response_zone1','var')+exist('snow_response_zone2','var')+exist('snow_response_zone3','var')
                case 1
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','InputUse');
                    clear snow_response_zone1
                case 2
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','snow_response_zone2','InputUse');
                    clear snow_response_zone1 snow_response_zone2
                case 3
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','snow_response_zone2','snow_response_zone3','InputUse');
                    clear snow_response_zone1 snow_response_zone2 snow_response3
            end
        else
            switch exist('snow_response_zone1','var')+exist('snow_response_zone2','var')+exist('snow_response_zone3','var')
                case 1
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','-append');
                    clear snow_response_zone1
                case 2
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','snow_response_zone2','-append');
                    clear snow_response_zone1 snow_response_zone2
                case 3
                    save([out_name(1:end-3),'mat'],'snow_response_zone1','snow_response_zone2','snow_response_zone3','-append');
                    clear snow_response_zone1 snow_response_zone2 snow_response3
            end
        end
    end
end
   
%% WRITE XLS FILE
% Sum all gravity effects
out.snow_total_dg = out.snow_zone1_dg+out.snow_zone2_dg+out.snow_zone3_dg;  % equal to zero if not computed
out.gw_total_dg = out.gw_zone1_dg+out.gw_zone2_dg+out.gw_zone3_dg;
out.sm_total_dg = sum(out.sm_zone1_dg,2)+sum(out.sm_zone2_dg,2)+sum(out.sm_zone3_dg,2);
out.total = sum(out.snow_total_dg,2)+sum(out.gw_total_dg,2)+sum(out.sm_total_dg,2);              % total gravity effect of all zones and compartments

ctime = clock;                                                              % current time
output_name = output_file{1};                                               % output file name
try 
    if output_file{3} == 1
        set(findobj('Tag','text_status_mLocal'),'String','Gravity response: writing (xls) output...');drawnow;
        table ={'mLocal Gravity Response'};
        if exist('warning_message','var')
            table(1,3) = {sprintf('%s',warning_message)};
        end
        if coor_system == 1
            table(2,1) = {'x coordinate (m)'};table(2,2) = num2cell(InputUse(1));
            table(3,1) = {'y coordinate (m)'};table(3,2) = num2cell(InputUse(2));
        else
            table(2,1) = {'Longitude (deg)'};table(2,2) = num2cell(InputUse(1));
            table(3,1) = {'Latitude_ellipsoid (deg)'};table(3,2) = num2cell(Input(2)); % Write the ellipsoidal not the spherical latitude
            table(3,3) = {'Latitude_sphere (deg)'};table(3,4) = num2cell(InputUse(2)); % Write the spherical latitude
        end
        table(4,1) = {'Height given'};table(4,2) = num2cell(Input(3));table(4,3) = {'sensor'};table(4,4) = num2cell(Input(4));
        table(5,1) = {'Height used (+sensor)'};table(5,2) = num2cell(InputUse(3));
        table(6,1) = {'Computation zones:'};
        table(7,1) = {'Zone 1 DEM:'};table(7,2) = {dem{1}};
        table(8,1) = {'Zone 2 DEM:'};table(8,2) = {dem{2}};
        table(9,1) = {'Zone 3 DEM:'};table(9,2) = {dem{3}};
        if coor_system == 1
            table(10,1) = {'Zone 1 Resolution (m):'};table(10,2) = num2cell(mean(resol1));table(10,3) = {'Radius (m):'};table(10,4) = {radius{1}};
            table(11,1) = {'Zone 2 Resolution (m):'};table(11,2) = num2cell(mean(resol2));table(11,3) = {'Radius (m):'};table(11,4) = {radius{2}};
            table(12,1) = {'Zone 3 Resolution (m):'};table(12,2) = num2cell(mean(resol3));table(12,3) = {'Radius (m):'};table(12,4) = {radius{3}};
        else
            table(10,1) = {'Zone 1 Resolution (deg):'};table(10,2) = num2cell(mean(resol1));table(10,3) = {'Radius (deg):'};table(10,4) = {radius{1}};
            table(11,1) = {'Zone 2 Resolution (deg):'};table(11,2) = num2cell(mean(resol2));table(11,3) = {'Radius (deg):'};table(11,4) = {radius{2}};
            table(12,1) = {'Zone 3 Resolution (deg):'};table(12,2) = num2cell(mean(resol3));table(12,3) = {'Radius (deg):'};table(12,4) = {radius{3}};
        end
        table(13,1) = {'Soil moisture (sm) computed for:'};table(13,2)= num2cell(sm_exist);table(13,3) = {'layers'};
        r = 14;
        for s = 1:size(sm,1)
            if sm{s,1} == 1 && ~isempty(sm{s,5})
                table(r,1) = {'SM layer'}; table(r,2) = num2cell(s); table(r,3) = {'from'}; table(r,4) = {sm{s,2}}; table(r,5) = {'to'}; table(r,6) = {sm{s,3}}; table(r,7) = {'Input file'}; table(r,8) =  {sm{s,5}}; table(r,9) = {'column'}; table(r,10) = {sm{s,6}}; table(r,11) = {'exclusion polygon:'}; table(r,12) =  {sm{s,4}};
                r=r+1;
            elseif sm{s,1} == 1 && isempty(sm{s,5})
                table(r,1) = {'SM layer'}; table(r,2) = num2cell(s); table(r,3) = {'from'}; table(r,4) = {sm{s,2}}; table(r,5) = {'to'}; table(r,6) = {sm{s,3}}; table(r,7) = {'computed for constant soil moisture, i.e. 1 Vol%'}; table(r,11) = {'exclusion polygon:'}; table(r,12) =  {sm{s,4}};
                r=r+1;
            end
        end
        if gw{1}==1 && ~isempty(gw{6})
           table(r,1) = {'Groundwater input:'}; table(r,2) = {gw{6}}; table(r,3) = {'column'}; table(r,4) = {gw{7}}; table(r,5) = {'specific yield'};table(r,6) = {gw{3}};
           if gw{4} == 1
               table(r,7) = {'groundwater copies the terrain in (m)'}; table(r,8) = {gw{5}}; table(r,9) = {'exclusion polygon'}; table(r,10) = {gw{2}}; table(r,11) = {'lower boundary = min. observed GW altitude - 1 m'};
               r = r + 1;
           else
               table(r,7) = {'groundwater copies the terrain in variable depth'}; table(r,9) = {'exclusion polygon'}; table(r,10) = {gw{2}}; table(r,11) = {'lower boundary = max. observed GW depth - 1 m'};
               r = r + 1;
           end
        else
           table(r,1) = {'No groundwater effect computed\n'}; 
           r = r + 1;
        end

        if snow{1}==1 && ~isempty(snow{6})
           table(r,1) = {'Snow input:'};table(r,2) = {snow{6}}; table(r,3) = {snow{7}};
           if snow{8} >= 1
               table(r,4) = {'snow density computed using loaded SWE'};table(r,5) = {'column:'}; table(r,6) = {snow{8}};
           else
               table(r,4) = {'snow density computed as follows:'}; table(r,5) = {snow{3}}; table(r,6) = {'+'}; table(r,7) = {snow{4}}; table(r,8) = {'*depth_mm/1000.'};
           end
           if snow{5} >0
               table(r,9) = {'Snow gradient used ='};table(r,10) =  {snow{5}};table(r,11) = {'mm per 100 m'};table(r,12) = {'exclusion polygon:'};table(r,13) = {snow{2}};
               r = r + 1;
           else
               table(r,9) = {'No snow gradient used'};table(r,12) = {'exclusion polygon:'};table(r,13) = {snow{2}};
               r = r + 1;
           end
        else
           table(r,1) = {'No snow effect computed'}; 
           r = r + 1;
        end
        table(r,1) = {'Computed'};table(r,2:7) = num2cell(ctime);r = r + 1;
        table(r,1) = {'Gravity effects (in nm/s^2):'};r = r + 1;
        table(r,1) = {'Matlab_Time'};table(r,2) = {'year'};table(r,3) = {'month'};table(r,4) = {'day'};table(r,5) = {'hour'};table(r,6) = {'minute'};table(r,7) = {'second'};table(r,8) = {'total_effect'};table(r,9) = {'total_sm'};table(r,10) = {'total_gw'};table(r,11) = {'total_snow'};
        c = 12;
        for z = 1:3
            for i = 1:size(sm,1)
                if sm{i,1} == 1
                    table(r,c) = {sprintf('%9s%1d_L%02d','sm_zone',z,i)};c = c + 1;
                end
            end
            if gw{1}==1 && ~isempty(gw{6})
                table(r,c) = {sprintf('%13s%1d','gw_zone',z)};c = c + 1;
            end
            if snow{1}==1 && ~isempty(snow{6})
                table(r,c) = {sprintf('%13s%1d','snow_zone',z)};c = c + 1;
            end
        end
        r = r + 1;
        if exist('ref_time','var')
            [year,month,day,hour,minute,second] = datevec(ref_time);
            for i = 1:length(ref_time)
                table(r,1:7) = num2cell([ref_time(i),year(i),month(i),day(i),hour(i),minute(i),round(second(i))]);
                table(r,8:11) = num2cell([out.total(i),out.sm_total_dg(i),out.gw_total_dg(i),out.snow_total_dg(i)]);
                c = 12;
                for z = 1:3                                                 % for all three zoned
                    for j = 1:size(sm,1)                                    % for all SM layers
                        if sm{j,1} == 1                                     % use only selected layers
                            switch z
                                case 1
                                    table(r,c) = num2cell(out.sm_zone1_dg(i,j));c = c + 1;
                                case 2
                                    table(r,c) = num2cell(out.sm_zone2_dg(i,j));c = c + 1;
                                case 3
                                    table(r,c) = num2cell(out.sm_zone3_dg(i,j));c = c + 1;
                            end
                        end
                    end
                    if gw{1}==1 && ~isempty(gw{6})
                        switch z
                            case 1
                                table(r,c) = num2cell(out.gw_zone1_dg(i));c = c + 1;
                            case 2
                                table(r,c) = num2cell(out.gw_zone2_dg(i));c = c + 1;
                            case 3
                                table(r,c) = num2cell(out.gw_zone3_dg(i));c = c + 1;
                        end
                    end
                    if snow{1}==1 && ~isempty(snow{6})
                        switch z
                            case 1
                                table(r,c) = num2cell(out.snow_zone1_dg(i));c = c + 1;
                            case 2
                                table(r,c) = num2cell(out.snow_zone2_dg(i));c = c + 1;
                            case 3
                                table(r,c) = num2cell(out.snow_zone3_dg(i));c = c + 1;
                        end
                    end
                end 
                r = r + 1;
            end
        else                                                                    % Write output for constant soil moisture layers
    %         fprintf(fid,'%13.6f  %8d  %6d',NaN,NaN,NaN);%datenum(ctime),ctime(1),ctime(2),ctime(3),ctime(4),ctime(5),round(ctime(6)));          % No time or date
            table(r,8:11) = num2cell([sum(out.sm_zone1_dg+out.sm_zone2_dg+out.sm_zone3_dg),sum(out.sm_zone1_dg+out.sm_zone2_dg+out.sm_zone3_dg),0,0]);
            c = 12;
            for i = 1:length(out.sm_zone1_dg)
                if sm{i,1} == 1
                    table(r,c) = num2cell(out.sm_zone1_dg(i));c = c + 1;
                end
            end
            for i = 1:length(out.sm_zone2_dg)
                if sm{i,1} == 1
                    table(r,c) = num2cell(out.sm_zone2_dg(i));c = c + 1;
                end
            end
            for i = 1:length(out.sm_zone3_dg)
                if sm{i,1} == 1
                    table(r,c) = num2cell(out.sm_zone3_dg(i));c = c + 1;
                end
            end
        end
        xlswrite([output_name(1:end-3),'xls'],table);
    end
catch
    fprintf('Could not write the xls file, see txt for all results\n');
    output_file{2} = 1;                                                     % turn of TXT if error occurs during XLS writing
end

%% WRITE TXT FILE
if output_file{2} == 1
    set(findobj('Tag','text_status_mLocal'),'String','Gravity response: writing (txt) output...');drawnow;
    fid = fopen([output_name(1:end-3),'txt'],'w');
    fprintf(fid,'%% mLocal Gravity Response');
    if exist('warning_message','var')
        fprintf(fid,' %s\n',warning_message);
    else
        fprintf(fid,'\n');
    end
    if coor_system == 1
        fprintf(fid,'%% x coordinate (m) = %12.3f\n',InputUse(1));
        fprintf(fid,'%% y coordinate (m) = %12.3f\n',InputUse(2));
    else
        fprintf(fid,'%% Longitude (deg) = %14.10f\n',Input(1));
        fprintf(fid,'%% Latitude_ellipsoid (deg) = %14.10f,\t Latitude_sphere (deg) = %14.10f\n',Input(2),InputUse(2));
    end
    fprintf(fid,'%% Height given = %9.3f  (sensor = %6.3f)\n',Input(3),Input(4));
    fprintf(fid,'%% Height used (+sensor) = %9.3f\n',InputUse(3));
    fprintf(fid,'%% Computation zones:\n');
    fprintf(fid,'%% Zone 1 DEM: %s\n',dem{1});
    fprintf(fid,'%% Zone 2 DEM: %s\n',dem{2});
    fprintf(fid,'%% Zone 3 DEM: %s\n',dem{3});
    if coor_system == 1
        fprintf(fid,'%% Zone 1 Resolution (m):%12.4f      Radius (m):%12.2f\n',mean(resol1),radius{1});
        fprintf(fid,'%% Zone 2 Resolution (m):%12.4f      Radius (m):%12.2f\n',mean(resol2),radius{2});
        fprintf(fid,'%% Zone 3 Resolution (m):%12.4f      Radius (m):%12.2f\n',mean(resol3),radius{3});
    else
        fprintf(fid,'%% Zone 1 Resolution (deg):%12.9f      Radius (deg):%12.9f\n',mean(resol1),radius{1});
        fprintf(fid,'%% Zone 2 Resolution (deg):%12.9f      Radius (deg):%12.9f\n',mean(resol2),radius{2});
        fprintf(fid,'%% Zone 3 Resolution (deg):%12.9f      Radius (deg):%12.9f\n',mean(resol3),radius{3});
    end
    fprintf(fid,'%% Soil moisture (sm) computed for %2d layer(s):\n',sm_exist);
    for s = 1:size(sm,1)
        if sm{s,1} == 1 && ~isempty(sm{s,5})
            fprintf(fid,'%% Soil moisture layer %1d (from %5.2f to %5.2f m): %s (column %2d) (exclusion polygon: %s)\n',s,sm{s,2},sm{s,3},sm{s,5},sm{s,6},sm{s,4});
        elseif sm{s,1} == 1 && isempty(sm{s,5})
            fprintf(fid,'%% Soil moisture layer %1d (from %5.2f to %5.2f m): computed for constant soil moisture, i.e. 1 Vol% (exclusion polygon: %s)\n',s,sm{s,2},sm{s,3},sm{s,4});
        end
    end
    if gw{1}==1 && ~isempty(gw{6})
       fprintf(fid,'%% Groundwater (gw) input: %s (column %2d),  specific yield = %3.1f%%,',gw{6},gw{7},gw{3});
       if gw{4} == 1
           fprintf(fid,' groundwater copies the terrain in %4.1f m (exclusion polygon: %s, lower boundary = min. observed GW altitude - 1 m)\n',gw{5},gw{2});
       else
           fprintf(fid,' groundwater copies the terrain in variable depth (exclusion polygon: %s, lower boundary = max. observed GW depth - 1 m)\n',gw{2});
       end
    else
       fprintf(fid,'%% No groundwater effect computed\n'); 
    end
    
    if snow{1}==1 && ~isempty(snow{6})
       fprintf(fid,'%% Snow input: %s (depth column %2d),',snow{6},snow{7});
       if snow{8} >= 1
           fprintf(fid,' snow density computed using loaded SWE (column %2d).',snow{8});
       else
           fprintf(fid,' snow density computed as follows: %5.2f + %5.2f*depth_mm/1000.',snow{3},snow{4});
       end
       if snow{5} >0
           fprintf(fid,' Snow gradient used = %5.2f mm per 100 m (exclusion polygon: %s)\n',snow{5}, snow{2});
       else
           fprintf(fid,' No snow gradient used (exclusion polygon: %s)\n', snow{2});
       end
    else
       fprintf(fid,'%% No snow effect computed\n'); 
    end
    fprintf(fid,'%% Computed %04d/%02d/%02d %02d:%02d:%02d\n',ctime(1),ctime(2),ctime(3),ctime(4),ctime(5),round(ctime(6)));
    fprintf(fid,'%% Gravity effects (in nm/s^2):\n');
    fprintf(fid,'%% %11s%10s%8s%14s%14s%14s%14s','Matlab_Time','Date','Time','total_effect','total_sm','total_gw','total_snow');
    for z = 1:3
        for i = 1:size(sm,1)
            if sm{i,1} == 1
                fprintf(fid,'%9s%1d_L%02d','sm_zone',z,i);
            end
        end
        if gw{1}==1 && ~isempty(gw{6})
            fprintf(fid,'%13s%1d','gw_zone',z);
        end
        if snow{1}==1 && ~isempty(snow{6})
            fprintf(fid,'%13s%1d','snow_zone',z);
        end
    end
    fprintf(fid,'\n');
    if exist('ref_time','var')
        [year,month,day,hour,minute,second] = datevec(ref_time);
        for i = 1:length(ref_time)
            fprintf(fid,'%13.6f  %4d%02d%02d  %02d%02d%02d',ref_time(i),year(i),month(i),day(i),hour(i),minute(i),round(second(i)));
            fprintf(fid,'%14.3f%14.3f%14.3f%14.3f',out.total(i),...
                    out.sm_total_dg(i),out.gw_total_dg(i),out.snow_total_dg(i));
            for z = 1:3
                for j = 1:size(sm,1)
                    if sm{j,1} == 1
                        switch z
                            case 1
                                fprintf(fid,'%14.3f',out.sm_zone1_dg(i,j));
                            case 2
                                fprintf(fid,'%14.3f',out.sm_zone2_dg(i,j));
                            case 3
                                fprintf(fid,'%14.3f',out.sm_zone3_dg(i,j));
                        end
                    end
                end
                if gw{1}==1 && ~isempty(gw{6})
                    switch z
                        case 1
                            fprintf(fid,'%14.3f',out.gw_zone1_dg(i,1));
                        case 2
                            fprintf(fid,'%14.3f',out.gw_zone2_dg(i,1));
                        case 3
                            fprintf(fid,'%14.3f',out.gw_zone3_dg(i,1));
                    end
                end
                if snow{1}==1 && ~isempty(snow{6})
                    switch z
                        case 1
                            fprintf(fid,'%14.3f',out.snow_zone1_dg(i,1));
                        case 2
                            fprintf(fid,'%14.3f',out.snow_zone2_dg(i,1));
                        case 3
                            fprintf(fid,'%14.3f',out.snow_zone3_dg(i,1));
                    end
                end
            end 
            fprintf(fid,'\n');
        end
        clear i j z;
    else                                                                    % Write output for constant soil moisture layers
        fprintf(fid,'%13.6f  %8d  %6d',NaN,NaN,NaN);%datenum(ctime),ctime(1),ctime(2),ctime(3),ctime(4),ctime(5),round(ctime(6)));          % No time or date
        fprintf(fid,'%14.3f%14.3f%14.3f%14.3f',...
            sum(out.sm_zone1_dg+out.sm_zone2_dg+out.sm_zone3_dg),sum(out.sm_zone1_dg+out.sm_zone2_dg+out.sm_zone3_dg),0,0);
        for i = 1:length(out.sm_zone1_dg)
            if sm{i,1} == 1
                fprintf(fid,'%14.3f',out.sm_zone1_dg(i));
            end
        end
        for i = 1:length(out.sm_zone2_dg)
            if sm{i,1} == 1
                fprintf(fid,'%14.3f',out.sm_zone2_dg(i));
            end
        end
        for i = 1:length(out.sm_zone3_dg)
            if sm{i,1} == 1
                fprintf(fid,'%14.3f',out.sm_zone3_dg(i));
            end
        end
        fprintf(fid,'\n');
    end
    fclose('all');
end

%% WRITE TSF FILE
if output_file{4} == 1
    if exist('ref_time','var')                                              % Write only for time variable effect
        set(findobj('Tag','text_status_mLocal'),'String','Gravity response: writing (tsf) output...');drawnow;
        fid = fopen([output_name(1:end-3),'tsf'],'w');
        fprintf(fid,'[TSF-file] v01.0\n\n');
        fprintf(fid,'[UNDETVAL] 1234567.89\n\n');
        out.total(isnan(out.total)) = 1234567.89;
        out.sm_zone1_dg(isnan(out.sm_zone1_dg)) = 1234567.89;
        out.sm_zone2_dg(isnan(out.sm_zone2_dg)) = 1234567.89;
        out.sm_zone3_dg(isnan(out.sm_zone3_dg)) = 1234567.89;
        out.sm_total_dg(isnan(out.sm_total_dg)) = 1234567.89;
        out.gw_zone1_dg(isnan(out.gw_zone1_dg)) = 1234567.89;
        out.gw_zone2_dg(isnan(out.gw_zone2_dg)) = 1234567.89;
        out.gw_zone3_dg(isnan(out.gw_zone3_dg)) = 1234567.89;
        out.gw_total_dg(isnan(out.gw_total_dg)) = 1234567.89;
        out.snow_zone1_dg(isnan(out.snow_zone1_dg)) = 1234567.89;
        out.snow_zone2_dg(isnan(out.snow_zone2_dg)) = 1234567.89;
        out.snow_zone3_dg(isnan(out.snow_zone3_dg)) = 1234567.89;
        out.snow_total_dg(isnan(out.snow_total_dg)) = 1234567.89;
        fprintf(fid,'[TIMEFORMAT] DATETIME\n\n');
        fprintf(fid,'[INCREMENT] %8.3f\n\n',(ref_time(2)-ref_time(1))*24*60*60);
        fprintf(fid,'[CHANNELS]\n');
        fprintf(fid,'  Location:mLocal:total_effect\n'); 
        fprintf(fid,'  Location:mLocal:total_sm\n'); 
        fprintf(fid,'  Location:mLocal:total_gw\n'); 
        fprintf(fid,'  Location:mLocal:total_snow\n'); 
        for z = 1:3
            for i = 1:size(sm,1)
                if sm{i,1} == 1
                    fprintf(fid,'  Location:mLocal:%s%1d_L%02d\n','sm_zone',z,i);
                end
            end
            if gw{1}==1 && ~isempty(gw{6})
                fprintf(fid,'  Location:mLocal:%s%1d\n','gw_zone',z);
            end
            if snow{1}==1 && ~isempty(snow{6})
                fprintf(fid,'  Location:mLocal:%s%1d\n','snow_zone',z);
            end
        end
        fprintf(fid,'\n');
        fprintf(fid,'[UNITS]\n  nm/s^2\n  nm/s^2\n  nm/s^2\n  nm/s^2\n');
        for z = 1:3
            for i = 1:size(sm,1)
                if sm{i,1} == 1
                    fprintf(fid,'  nm/s^2\n');
                end
            end
            if gw{1}==1 && ~isempty(gw{6})
                fprintf(fid,'  nm/s^2\n');
            end
            if snow{1}==1 && ~isempty(snow{6})
                fprintf(fid,'  nm/s^2\n');
            end
        end
        fprintf(fid,'\n[COMMENT]\n\n');
        fprintf(fid,'[COUNTINFO] %8.0f\n\n',length(ref_time));
        fprintf(fid,'[DATA]\n');
        
        [year,month,day,hour,minute,second] = datevec(ref_time);
        for i = 1:length(ref_time)
            fprintf(fid,'%4d %02d %02d  %02d %02d %02d',year(i),month(i),day(i),hour(i),minute(i),round(second(i)));
            fprintf(fid,'%14.3f%14.3f%14.3f%14.3f',out.total(i),...
                    out.sm_total_dg(i),out.gw_total_dg(i),out.snow_total_dg(i));
            for z = 1:3
                for j = 1:size(sm,1)
                    if sm{j,1} == 1
                        switch z
                            case 1
                                fprintf(fid,'%14.3f',out.sm_zone1_dg(i,j));
                            case 2
                                fprintf(fid,'%14.3f',out.sm_zone2_dg(i,j));
                            case 3
                                fprintf(fid,'%14.3f',out.sm_zone3_dg(i,j));
                        end
                    end
                end
                if gw{1}==1 && ~isempty(gw{6})
                    switch z
                        case 1
                            fprintf(fid,'%14.3f',out.gw_zone1_dg(i,1));
                        case 2
                            fprintf(fid,'%14.3f',out.gw_zone2_dg(i,1));
                        case 3
                            fprintf(fid,'%14.3f',out.gw_zone3_dg(i,1));
                    end
                end
                if snow{1}==1 && ~isempty(snow{6})
                    switch z
                        case 1
                            fprintf(fid,'%14.3f',out.snow_zone1_dg(i,1));
                        case 2
                            fprintf(fid,'%14.3f',out.snow_zone2_dg(i,1));
                        case 3
                            fprintf(fid,'%14.3f',out.snow_zone3_dg(i,1));
                    end
                end
            end 
            fprintf(fid,'\n');
        end
        fclose('all');
    end
end

end

