function y = takeda(x, xold, yy, k1, xmax, xmin, yold, k2, ymax, ymin)
  % 変位xに対する修正武田モデルの履歴復元力を計算
  % x:変位　xold:前ステップの変位　yy:降伏耐力
  % k1:初期剛性　xmax:変位の最大値　xmin:変位の最小値
  % （xmaxとxminの初期値は正負の降伏変位とする）
  % yold:前ステップの復元力　k2:降伏後の二次剛性
  % ymax:復元力の最大値　ymin:復元力の最小値

  EPS = 0.0001;              % この値未満であれば変化していないと判断
  xy = yy / k1;
  if abs(xmax - xy) > EPS    % 正側が非線形領域になっていれば除荷剛性を変える
    krp = k1 * (xmax / xy) ^ (-0.4);
  else
    krp = k1;
  end
  if abs(xmin + xy) > EPS    % 負正側が非線形領域になっていれば除荷剛性を変える
    krm = k1 * (-xmin / xy) ^ (-0.4);
  else
    krm = k1;
  end
  fLp = krp * (x - xold) + yold;
  fLm = krm * (x - xold) + yold;
  x0p = xold - yold / krp;
  x0m = xold - yold / krm;
  fplus = k2 * (x - xy) + yy;
  fminus = k2 * (x + xy) - yy;
  if x > xold
   if abs(xmax - x0m) < EPS
    fp = yy;
   else
    fp = (x - x0m) * ymax / (xmax - x0m);
   end
   if abs(xmax - xold) < EPS
    fd = yy;
   else
    fd = (x - xold) * (ymax - yold) / (xmax - xold) + yold;
   end
   fm = median( [fLm fp fd] );
   y = min(fplus, fm);
  else
   if abs(xmin - x0p) < EPS
    fp = -yy;
   else
    fp = (x - x0p) * ymin / (xmin - x0p);
   end
   if abs(xmin - xold) < EPS
    fd = -yy;
   else
    fd = (x - xold) * (ymin - yold) / (xmin - xold) + yold;
   end
   fm = median( [fLp fp fd] );
   y = max(fminus, fm);
  end

end
