gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/+markov/@bath/corr.m

    function [g, s] = corr(b, x)
% CORR  Bath spectral correlation tensor.
%
%  [g, s] = corr(bath, dH)  % returns Gamma(omega0 * dH)/omega0, see below
%
%  Returns the bath spectral correlation tensor Gamma evaluated at omega0 * dH:
%
%    Gamma(omega0 * dH)/omega0 == 0.5*g +i*s

% Ville Bergholm 2009-2010


if (nargin ~= 2)
  error('Need a bath object and dH.')
end

tol = 1e-8;
max_w = 0.1; % TODO justify


% assume parameters are set and lookup table computed
%s = interp1(b.dH, b.s_table, x, 'linear', 0);

% binary search for the interval in which x falls
first = 1;
last = length(b.dH);
while (first+1 ~= last)
  pivot = round((first + last)/2);
  if (x < b.dH(pivot))
    last = pivot;
  else
    first = pivot;
  end
end
ee = b.dH([first last]); 
tt = b.gs_table(:, [first last]);
% now x is in [ee(1), ee(2))

gap = ee(2) - ee(1);
d1 = abs(x - ee(1));
d2 = abs(x - ee(2));

% close enough to either endpoint?
if (d1 <= tol)
  temp = b.gs_table(:, first);

elseif (d2 <= tol)
  temp = b.gs_table(:, last);

elseif (gap <= max_w + tol) % short enough gap to interpolate?
  temp = interpolate(ee, tt, x);

else % compute a new point
  if (gap <= 2*max_w)
    p = ee(1) +gap/2; % halfway
    if (x < p)
      idx = 2; % which ee p will replace
    else
      idx = 1;
    end
  elseif (d1 <= max_w)
    p = ee(1)+max_w;
    idx = 2;
  elseif (d2 <= max_w)
    p = ee(2)-max_w;
    idx = 1;
  else
    p = x;
    idx = 1;
  end

  % compute new g, S values at p and insert them into the table
  temp = [0; 0];
  temp(2) = S_func(b, p);

  if (abs(p) <= tol)
    temp(1) = b.g0; % limit at p == 0
  else
    temp(1) = b.g_func(p) .* b.cut_func(p);
  end

  b.dH = [b.dH(1:first), p, b.dH(last:end)];
  b.gs_table = [b.gs_table(:, 1:first), temp, b.gs_table(:, last:end)];

  % now interpolate the required value
  ee(idx) = p;
  tt(:, idx) = temp;
  temp = interpolate(ee, tt, x);
end

g = temp(1);
s = temp(2);

end

function y = interpolate(ee, tt, x)
% interp1 does way too many checks
  y = tt(:,1) + ((x - ee(1))/(ee(2) - ee(1)))*(tt(:,2) - tt(:,1));
end