find_opt_L.m

function L_opt = find_opt_L(SF,rf_lonlat)

L_test  = 0.5:0.5:100;

D  = [];
D2 = [];

STDs = [];

% Reduce spatial resolution
sires= [500 500];

% in case z is already small (smaller than 500x500)
if size(SF.z,1)<= 500
    sires(1)= size(SF.z,1);
end

if size(SF.z,2)<= 500
    sires(2) = size(SF.z,2);
end

lat2    = imresize(SF.lat,sires);
lon2    = imresize(SF.lon,sires);
z2      = imresize(SF.z, sires);

% Find values of reduction factor/ water levels in rfgrid
P = [];
xy = [lon2(:) lat2(:)];
for ii = 1: size(rf_lonlat,1)

    [~,P(ii,1)] = min(abs(distance(rf_lonlat(ii,2),rf_lonlat(ii,1),xy(:,2),xy(:,1))));
end

for c = 1: length(L_test)

    L           = L_test(c);

    rfgrid      = red_fac_grid_FOL(lon2,lat2,z2,rf_lonlat,L);

    % differences between the optimal interpolation grid and the values at
    % the scalar positions
    D{c,1}      = rfgrid(P) - rf_lonlat(:,3);

    D2(c,1)     = mean(D{c,1});

    STDs(c,1)   = std(std(rfgrid));

end

% normalization
STDs_n      = (STDs - min(STDs))./(max(STDs) - min(STDs));
D2_n        = (D2 - min(D2))./(max(D2) - min(D2));

% criteria for the optimal value of the interpolation length (L)
[~,L_opt]   = max(STDs_n./D2_n);
L_opt       = L_test(L_opt);