function [mean_bv,std_bv,bv]=bond_valence_data(ion1,ion2,R,varargin)

if nargin==3
    load('bond_valence_values.mat');
else
    Ion_1=varargin{1};
    Ion_2=varargin{2};
    R0=varargin{3};
    b=varargin{4};
end

ind1=find(ismember(Ion_1,ion1));
ind2=find(ismember(Ion_2,ion2));
ind=find(ismember(ind1,ind2));
ind=ind1(ind);

if numel(ind)==0
    ind1=find(ismember(Ion_2,ion1));
    ind2=find(ismember(Ion_1,ion2));
    ind=find(ismember(ind1,ind2));
    ind=ind1(ind);
end

if numel(ind)==0
    disp('Could not find any matching pair...')
end

for i=1:numel(ind)
    bv(i)=exp((R0(ind(i))-R)/b(i));
end

mean_bv=mean(bv);
std_bv=std(bv);