function transitionRate = obtainTransitionRate(position, ...
potentialEnergy,wavefunc,eigenEnergy,temperature,hbar,v_sound,density)
%
steplength = position(2) -position(1);
temperature = temperature*0.025852206869652155/300;
distribution = @(x)(exp(hbar*x/temperature)-1)^-1;
%%%%
diff_eigen = abs(eigenEnergy.'-eigenEnergy);
for ii = 1:size(diff_eigen)
diff_eigen(ii,ii) = nan;
end
diff_max = max(max(diff_eigen));
diff_min = min(min(diff_eigen)); % the smallest one is not zero
s = 1; dos_E = @(x)(x.^s*(diff_max/hbar)^(1-s));
%%%%
tr_temp = zeros(size(wavefunc,2));
% for ii = 1:size(wavefunc,2)
for ii = 2:size(wavefunc,2)
for jj = 1:ii-1
diff_E = abs(eigenEnergy(ii)-eigenEnergy(jj))/hbar;
% coeff1 = diff_E/(2*pi*hbar*v_sound^3*density);
coeff1 = dos_E(diff_E)/(2*pi*hbar*v_sound^3*density);
coeff2 = abs(trapz(position(1:end-1),conj(wavefunc(1:end-1,jj)) ...
.*diff(potentialEnergy)/steplength.*wavefunc(1:end-1,ii)))^2;
tr_temp(ii,jj) = coeff1*coeff2*(distribution(diff_E)+1);
tr_temp(jj,ii) = coeff1*coeff2*distribution(diff_E);
end
end
transitionRate = zeros(size(wavefunc,2));
for ii = 1:size(transitionRate,1)
for jj = 1:size(transitionRate,1)
transitionRate(ii,jj) = tr_temp(jj,ii);
end
end
for ii = 1:size(transitionRate,1)
% transitionRate(ii,ii) = -(sum(tr_temp(ii,:))-tr_temp(ii,ii));
transitionRate(ii,ii) = -(sum(tr_temp(ii,:))-tr_temp(ii,ii));
end
end