function v = shoulder(ii,iia,iib)
  % Computes a smooth function that is 0 at {ii = iia}
  % and 1 at {ii = iib}, monotonic between the two,
  % and constant outside them.
  
  if iia > iib
    % Flip the argumens:
    iim = (iia+iib)/2;
    ii = iim - (ii - iim);
    tmp = iia; iia = iib; iib = tmp;
  end
  sigma = 3.8; % Sharpness factor of ramps.
  if ii <= iia
    v = 0.0;
  elseif ii >= iib
    v = 1.0;
  else
    d = ii - (iia + iib)/2;
    x = sigma*d/((iib - iia)/2);
    v = (1 + erf(x))/2;
    % Debugging:
    if (ii == iia+1) || (ii == iib-1)
      v1 = v
    end
  end