// compsim.C Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 09:35:51 tac>
//
// 区画モデルによる血中濃度のシミュレーション。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法:g++ -O -s -W compsim compsim.C -I/usr/local/include
//             -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(mat.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "compsim  Version 7.42        Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:compsim  [<オプション> …]",
  "機能:区画モデルによる血中濃度のシミュレーションを行う。",
  "<オプション>",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -f0=<吸収率>:吸収率を指定する。(0＜吸収率≦100)",
  "  -h[=<ページャ名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -k12=<分布速度定数>:分布速度定数を指定する。(0＜分布速度定数)",
  "  -k21=<分布速度定数>:分布速度定数を指定する。(0＜分布速度定数)",
  "  -ka=<吸収速度定数>:吸収速度定数を指定する。(0＜吸収速度定数)",
  "  -ke=<排出速度定数>:排出速度定数を指定する。(0＜排出速度定数)",
  "  -kr=<放出速度定数>:放出速度定数を指定する。(0＜放出速度定数)",
  "  -model[={iv1|iv1if|iv2|iv22|iv2T|iv2e|iv2if|oi1r0|oi1r1|po1",
  "    |po1T|po1e|po1lag|po1r0|po1r1|po2|po2lag}]:モデルを指定します。",
  "  -q=<投与量1>[,<投与量2>,…]:投与量を指定します。(0＜投与量)",
  "  -t=<投与時間1>[,<投与時間2>,…]:投与時間を指定する。(0≦投与時間)",
  "  -tif=<持続投与時間>:持続投与時間を指定する。(0＜持続投与時間)",
  "  -tlag=<ラグタイム>:ラグタイムを指定する。(0≦ラグタイム)",
  "  -vd=<分布容積>:分布容積を指定する。(0＜分布容積)"
};

const string Title = "=== 区画モデルによる血中濃度シミュレーション ==="; // 表題
constexpr double F0def = 100.0; // 吸収率省略時値
constexpr double K12def = 0.4;  // k12省略時値
constexpr double K21def = 0.6;  // k21省略時値
constexpr double KAdef = 1.0;   // ka省略時値
constexpr double KEdef = 0.5;   // ke省略時値
constexpr double KRdef = 0.1;   // kr省略時値
constexpr double TIFdef = 1.0;  // tif省略時値
constexpr double TLAGdef = 0.0; // tlag省略時値
constexpr double VDdef = 1.0;   // vd省略時値

enum modeltype                  // 区画モデルの種類
{
  modelnon,                     // 指定なし
  modeliv1,                     // 静注1区画モデル
  modeliv1if,                   // 静注1区画持続注入モデル
  modeliv2,                     // 静注2区画モデル
  modeliv22,                    // 静注2区画モデル組織区画関数
  modeliv2T,                    // 静注2区画モデル全体関数
  modeliv2e,                    // 静注2区画モデル排出関数
  modeliv2if,                   // 静注2区画持続注入モデル
  modeloi1r0,                   // 軟膏1区画0次徐放モデル
  modeloi1r1,                   // 軟膏1区画1次徐放モデル
  modelpo1,                     // 内服1区画モデル
  modelpo1T,                    // 内服1区画モデル体全体関数
  modelpo1e,                    // 内服1区画モデル排出関数
  modelpo1lag,                  // 内服1区画ラグタイムモデル
  modelpo1r0,                   // 内服1区画0次徐放モデル
  modelpo1r1,                   // 内服1区画1次徐放モデル
  modelpo2,                     // 内服2区画モデル
  modelpo2lag,                  // 内服2区画ラグタイムモデル
  modelukm1,                    // 透析1区画モデル
  modelukm2,                    // 透析2区画モデル
  modeletc                      // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_f0;                  // 吸収率(%)
  double v_k12;                 // 分布速度定数
  double v_k21;                 // 分布速度定数
  double v_ka;                  // 吸収速度定数
  double v_ke;                  // 排出速度定数
  double v_kr;                  // 放出速度定数
  modeltype v_model;            // 区画モデルの種類
  lstr qlst;                    // 投与量リスト
  vecdbl qvec;                  // 投与量ベクトル
  lstr tlst;                    // 投与時間リスト
  vecdbl tvec;                  // 投与時間ベクトル
  int v_nt;                     // 投与回数
  double v_tif;                 // 持続投与時間
  double v_tlag;                // ラグタイム
  double v_vd;                  // 分布容積
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double f0() const             // 吸収率(%)を返す
  { return v_f0; }
  double k12() const            // 分布速度定数を返す
  { return v_k12; }
  double k21() const            // 分布速度定数を返す
  { return v_k21; }
  double ka() const             // 吸収速度定数を返す
  { return v_ka; }
  double k0(size_t idx) const   // 0次放出速度を返す
  { return (idx < qvec.size()) ? qvec[idx] / tif() : 0.0; }
  double ke() const             // 排出速度定数を返す
  { return v_ke; }
  double kr() const             // 放出速度定数を返す
  { return v_kr; }
  modeltype model() const       // 区画モデルの種類を返す
  { return v_model; }
  const string modelname() const; // 区画モデルの名称を返す
  double q(size_t idx) const      // 投与量を返す
  { return (idx < qvec.size()) ? qvec[idx] : 0.0; }
  double t(size_t idx) const    // 投与時間を返す
  { return (idx < tvec.size()) ? tvec[idx] : 0.0; }
  int nt() const                // 投与回数を返す
  { return v_nt; }
  double tif() const            // 持続投与時間を返す
  { return v_tif; }
  double tlag() const           // ラグタイムを返す
  { return v_tlag; }
  double vd() const             // 分布容積を返す
  { return v_vd; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                  // オプションのチェックとセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_f0(F0def), v_k12(K12def), v_k21(K21def), v_ka(KAdef), v_ke(KEdef), v_kr(KRdef), v_model(modeletc), qlst(), qvec(), tlst(), tvec(), v_nt(0), v_tif(TIFdef), v_tlag(TLAGdef), v_vd(VDdef)
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "f")            // オプション指定ファイルから読み込み
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "f0")      // 吸収率の指定
  {
    v_f0 = opt[1].empty() ? F0def : strtodf(opt[1]);
    if (v_f0 <= 0.0 || 100.0 < v_f0) err = badopt(pname(), optstr, USAGE, 0.0, 100.0);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "k12")     // 分布速度定数の指定
  {
    v_k12 = opt[1].empty() ? K12def : strtodf(opt[1]);
    if (v_k12 <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "k21")     // 分布速度定数の指定
  {
    v_k21 = opt[1].empty() ? K21def : strtodf(opt[1]);
    if (v_k21 <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ka")      // 吸収速度定数の指定
  {
    v_ka = opt[1].empty() ? KAdef : strtodf(opt[1]);
    if (v_ka <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ke")      // 排出速度定数の指定
  {
    v_ke = opt[1].empty() ? KEdef : strtodf(opt[1]);
    if (v_ke <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "kr")      // 放出速度定数の指定
  {
    v_kr = opt[1].empty() ? KRdef : strtodf(opt[1]);
    if (v_kr <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "model")   // 区画モデルの指定
  {
    if (opt[1].empty() || opt[1] == "iv1") v_model = modeliv1;
    else if (opt[1] == "iv1if") v_model = modeliv1if;
    else if (opt[1] == "iv2") v_model = modeliv2;
    else if (opt[1] == "iv22") v_model = modeliv22;
    else if (opt[1] == "iv2T") v_model = modeliv2T;
    else if (opt[1] == "iv2e") v_model = modeliv2e;
    else if (opt[1] == "iv2if") v_model = modeliv2if;
    else if (opt[1] == "oi1r0") v_model = modeloi1r0;
    else if (opt[1] == "oi1r1") v_model = modeloi1r1;
    else if (opt[1] == "po1") v_model = modelpo1;
    else if (opt[1] == "po1T") v_model = modelpo1T;
    else if (opt[1] == "po1e") v_model = modelpo1e;
    else if (opt[1] == "po1lag") v_model = modelpo1lag;
    else if (opt[1] == "po1r0") v_model = modelpo1r0;
    else if (opt[1] == "po1r1") v_model = modelpo1r1;
    else if (opt[1] == "po2") v_model = modelpo2;
    else if (opt[1] == "po2lag") v_model = modelpo2lag;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "q")       // 投与量の指定
  {
    if (strsplit(qlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "t")       // 投与時間の指定
  {
    if (strsplit(tlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "tif")     // 持続投与時間の指定
  {
    v_tif = opt[1].empty() ? TIFdef : strtodf(opt[1]);
    if (v_tif <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "tlag")    // ラグタイムの指定
  {
    v_tlag = opt[1].empty() ? TLAGdef : strtodf(opt[1]);
    if (v_tlag < 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "vd")      // 分布容積の指定
  {
    v_vd = opt[1].empty() ? VDdef : strtodf(opt[1]);
    if (v_vd <= 0.0) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check()             // オプションのチェックとセット
{
  v_nt = tlst.size();
  int i = 0;
  if (v_nt < 1)                 // 投与時間が指定されていない時
  { v_nt = 1;  tvec.resize(1, 0.0); }
  else
  {
    tvec.resize(v_nt);
    double tpre = 0.0;
    for (lstr_citr pos = tlst.begin(); pos != tlst.end(); ++pos, ++i)
    {
      double t = strtodf(*pos);
      if (t < tpre || t < 0.0 || t == BLANKf) return disperr(pname(), "投与時間の指定が不適切！");
      tvec[i] = t;  tpre = t;
    }
  }
  if (qlst.size() > size_t(v_nt)) return disperr(pname(), "投与量指定が多すぎる！");
  qvec.resize(v_nt, 1.0);       // 省略時投与量を1にする
  i = 0;
  for (lstr_citr pos = qlst.begin(); pos != qlst.end(); ++pos, ++i)
  {
    qvec[i] = strtodf(*pos);
    if (qvec[i] <= 0.0 || qvec[i] == BLANKf) return disperr(pname(), "投与量が不適切！"); 
  }
  for (i = i; i < v_nt; ++i) qvec[i] = qvec[0]; // 残りの投与量を自動設定する
  if (v_model == modeloi1r0 || v_model == modeloi1r1 || v_model == modelpo1 || v_model == modelpo1T || v_model == modelpo1e || v_model == modelpo1lag || v_model == modelpo1r0 || v_model == modelpo1r1 || v_model == modelpo2 || v_model == modelpo2lag)
  {
    if (v_ka < v_ke)  return disperr(pname(), "kaとkeの指定が不適切！");
    if (v_model == modeloi1r1 || v_model == modelpo1r1)
    {
      if (v_ke < v_kr)  return disperr(pname(), "keとkrの指定が不適切！");
    }
  }
  if (v_model == modeliv2 || v_model == modeliv2if || v_model == modelpo2 || v_model == modelpo2lag)
  {
    if (sgn(v_k12 - v_k21) == 0.0)  return disperr(pname(), "k11とk21の指定が不適切！");
  }
  return 0;                     // エラーなし
}

const string option::modelname() const // 区画モデルの名称を返す
{
  switch (model())
  {
  case modeliv1if : return "静注1区画持続注入モデル";
  case modeliv2 : return "静注2区画モデル";
  case modeliv22 : return "静注2区画モデル組織区画関数";
  case modeliv2T : return "静注2区画モデル全体関数";
  case modeliv2e : return "静注2区画モデル排出関数";
  case modeliv2if : return "静注2区画持続注入モデル";
  case modeloi1r0 : return "軟膏1区画0次徐放モデル(放出速度一定)";
  case modeloi1r1 : return "軟膏1区画1次徐放モデル";
  case modelpo1 : return "内服1区画モデル";
  case modelpo1T : return "内服1区画モデル体全体関数";
  case modelpo1e : return "内服1区画モデル排出関数";
  case modelpo1lag : return "内服1区画ラグタイムモデル";
  case modelpo1r0 : return "内服1区画0次徐放モデル(放出速度一定)";
  case modelpo1r1 : return "内服1区画1次徐放モデル";
  case modelpo2 : return "内服2区画モデル";
  case modelpo2lag : return "内服2区画ラグタイムモデル";
  case modeliv1 :
  default : return "静注1区画モデル";
  }
}

class compsimclass              // メインクラス
{
private:
  void plot_iv1(const option& opt) const; // 静注1区画モデルによる血中濃度曲線出力
  void plot_iv1if(const option& opt) const;
  // 静注1区画持続注入モデルによる血中濃度曲線出力
  void plot_iv2(const option& opt) const; // 静注2区画モデルによる血中濃度曲線出力
  void plot_iv2if(const option& opt) const;
  // 静注2区画持続注入モデルによる血中濃度曲線出力
  void plot_oi1r0(const option& opt) const;
  // 軟膏1区画0次徐放モデルによる血中濃度曲線出力
  void plot_oi1r1(const option& opt) const;
  // 軟膏1区画1次徐放モデルによる血中濃度曲線出力
  void plot_po1(const option& opt) const; // 内服1区画モデルによる血中濃度曲線出力
  void plot_po1lag(const option& opt) const;
  // 内服1区画ラグタイムモデルによる血中濃度曲線出力
  void plot_po1r0(const option& opt) const;
  // 内服1区画0次徐放モデルによる血中濃度曲線出力
  void plot_po1r1(const option& opt) const;
  // 内服1区画1次徐放モデルによる血中濃度曲線出力
  void plot_po2(const option& opt) const;
  // 内服2区画モデルによる血中濃度曲線出力
  void plot_po2lag(const option& opt) const;
  // 内服2区画ラグタイムモデルによる血中濃度曲線出力
  void plot_comp(const option& opt) const; // 血中濃度曲線出力
  double fsigmac(const option& opt, double t) const; // 累積血中濃度関数
  double fcomp(const option& opt, double t, double q) const; // 血中濃度関数
  double iv2para(double& a, double& k1, double& b, double& k2, double& d, double q, double vd, double ke, double k12, double k21) const;
  // 静注2区画モデルのパラメーターを返す
  double iv2ifpara(double& a, double& k1, double& b, double& k2, double q, double vd, double ke, double k12, double k21, double tif) const;
  // 静注2区画持続注入モデルのパラメーターを返す
public:
  double po2para(double& a, double& k1, double& b, double& k2, double q, double f0, double vd, double ka, double ke, double k12, double k21) const;
  // 内服2区画モデルのパラメーターを返す
  compsimclass()                // コンストラクタ
  {}
  ~compsimclass()               // デストラクタ
  {}
  void put(option& opt);     // 区画モデルによる血中濃度曲線を出力する
};

void compsimclass::plot_iv1(const option& opt) const
// 静注1区画モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=(q/Vd)exp(-ke・t)\n";
  cout << "%分布容積Vd=" << opt.vd() << "  排泄速度定数ke=" << opt.ke() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double ee = exp(-opt.ke() * tint);
    double cmin = c0 * ee / (1.0 - ee);
    double cmax = c0 / (1.0 - ee);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_iv1if(const option& opt) const
// 静注1区画持続注入モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=(q/Vd)[exp{-ke(t-tif)}-exp(-ke・t)]\n";
  cout << "%分布容積Vd=" << opt.vd() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%持続注入時間tif=" << opt.tif() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < opt.nt(); ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "  注入速度k0=" << opt.k0(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q / opt.vd();
    double a = c0 / opt.tif() / opt.ke();
    double cmean = c0 / opt.ke() / tint;
    double ee = exp(-opt.ke() * tint);
    double tt = tint - opt.tif();
    double cmin = a * (exp(-opt.ke() * tt) - ee) / (1.0 - ee);
    double cmax = a * (1.0 - exp(-opt.ke() * opt.tif())) / (1.0 - ee);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_iv2(const option& opt) const
// 静注2区画モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A・exp(-α・t)+B・exp(-β・t)\n";
  cout << "%血液区画分布容積Vc=" << opt.vd() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%分布速度定数k12=" << opt.k12() << "  分布速度定数k21=" << opt.k21() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double a = 0.0, k1 = 0.0, b = 0.0, k2 = 0.0, d = 0.0;
    double c0 = iv2para(a, k1, b, k2, d, q, opt.vd(), opt.ke(), opt.k12(), opt.k21());
    double cmean = c0 / opt.ke() / tint;
    double e1 = exp(-k1 * tint), e2 = exp(-k2 * tint);
    double cmin = a * e1 / (1.0 - e1) + b * e2 / (1.0 - e2);
    double cmax = a / (1.0 - e1) + b / (1.0 - e2);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / k2 << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / k2 << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_iv2if(const option& opt) const
// 静注2区画持続注入モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A[exp{-α(t-tif)}-exp(-α・t)]+B[exp{-β(t-tif)}-exp(-β・t)]\n";
  cout << "%血液区画分布容積Vc=" << opt.vd() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%分布速度定数k12=" << opt.k12() << "  分布速度定数k21=" << opt.k21() << "\n";
  cout << "%持続注入時間tif=" << opt.tif() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "  注入速度k0=" << opt.k0(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double a = 0.0, k1 = 0.0, b = 0.0, k2 = 0.0;
    double c0 = iv2ifpara(a, k1, b, k2, q, opt.vd(), opt.ke(), opt.k12(), opt.k21(), opt.tif());
    double cmean = c0 / opt.ke() / tint;
    double e1 = exp(-k1 * tint), e2 = exp(-k2 * tint);
    double tt = tint - opt.tif();
    double cmin = a * (exp(-k1 * tt) - e1) / (1.0 - e1) + b * (exp(-k2 * tt) - e2) / (1.0 - e2);
    double cmax = a * (1.0 - exp(-k1 * opt.tif())) / (1.0 - e1) + b * (1.0 - exp(-k2 * opt.tif())) / (1.0 - e2);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / k2 << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / k2 << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_oi1r0(const option& opt) const
// 軟膏1区画0次徐放モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A[ka{exp(-ke・T)-exp(-ke・t)}-ke{exp(-ka・T)-exp(-ka・t)}]\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%放出時間tr=" << opt.tif() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "  放出速度k0=" << opt.k0(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double a = c0 / opt.tif() / opt.ke() / (opt.ka() - opt.ke());
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double tt = tint - opt.tif();
    double cmin = a * (opt.ka() * (exp(-opt.ke() * tt) - ee) / (1.0 - ee) - opt.ke() * (exp(-opt.ka() * tt) - ea) / (1.0 - ea));
    double cmax = a * (opt.ka() * (1.0 - exp(-opt.ke() * opt.tif())) / (1.0 - ee) - opt.ke() * (1.0 - exp(-opt.ka() * opt.tif())) / (1.0 - ea));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_oi1r1(const option& opt) const
// 軟膏1区画1次徐放モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A{(ka-kr)exp(-kr・tif-ke・T)-(ke-kr)exp(-kr・tif-ka・T)+(ke-kr)exp(-ka・t)-(ka-kr)exp(-ke・t)}\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%放出速度定数kr=" << opt.kr() << "  放出時間tr=" << opt.tif() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double k1 = opt.ka() - opt.kr(), k2 = opt.ke() - opt.kr();
    double k3 = opt.ka() - opt.ke();
    double a = c0 * opt.kr() * opt.ka() / k1 / k2 / k3;
    double tt = tint - opt.tif();
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double eet = exp(-opt.ke() * tt), eat = exp(-opt.ka() * tt);
    double err = exp(-opt.kr() * opt.tif());
    double cmin = a * (k1 * err * eet / (1.0 - ee) - k2 * err * eat / (1.0 - ea) + k2 * ea / (1.0 - ea) - k1 * ee / (1.0 - ee));
    double cmax = a * (k1 * err / (1.0 - ee) - k2 * err / (1.0 - ea) + k2 * exp(-opt.ka() * opt.tif()) / (1.0 - ea) - k1 * exp(-opt.ke() * opt.tif()) / (1.0 - ee));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po1(const option& opt) const
// 内服1区画モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A{exp(-ke・t)-exp(-ka・t)}\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double a = c0 * opt.ka() / (opt.ka() - opt.ke());
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double cmin = a * (ee / (1.0 - ee) - ea / (1.0 - ea));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po1lag(const option& opt) const
// 内服1区画ラグタイムモデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A[exp{-ke(t-tlag)}-exp{-ka(t-tlag)}]\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%ラグタイムtlag=" << opt.tlag() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double a = c0 * opt.ka() / (opt.ka() - opt.ke());
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double cmin = a * (ee / (1.0 - ee) - ea / (1.0 - ea));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() + opt.tlag() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() + opt.tlag() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po1r0(const option& opt) const
// 内服1区画0次徐放モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A[ka{exp(-ke・T)-exp(-ke・t)}-ke{exp(-ka・T)-exp(-ka・t)}]\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%放出時間tr=" << opt.tif() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "  放出速度k0=" << opt.k0(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double a = c0 / opt.tif() / opt.ke() / (opt.ka() - opt.ke());
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double tt = tint - opt.tif();
    double cmin = a * (opt.ka() * (exp(-opt.ke() * tt) - ee) / (1.0 - ee) - opt.ke() * (exp(-opt.ka() * tt) - ea) / (1.0 - ea));
    double cmax = a * (opt.ka() * (1.0 - exp(-opt.ke() * opt.tif())) / (1.0 - ee) - opt.ke() * (1.0 - exp(-opt.ka() * opt.tif())) / (1.0 - ea));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%定常状態における最高血中濃度=" << cmax << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.ke() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.ke() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po1r1(const option& opt) const
// 内服1区画1次徐放モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A{(ka-kr)exp(-kr・t)+(ke-kr)exp(-ka・t)-(ka-kr)exp(-ke・t)}\n";
  cout << "%分布容積Vd=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%放出速度定数kr=" << opt.kr() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double c0 = q * opt.f0() * 0.01 / opt.vd();
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double k1 = opt.ka() - opt.kr(), k2 = opt.ke() - opt.kr();
    double k3 = opt.ka() - opt.ke();
    double a = c0 * opt.kr() * opt.ka() / k1 / k2 / k3;
    double ee = exp(-opt.ke() * tint), ea = exp(-opt.ka() * tint);
    double er = exp(-opt.kr() * tint);
    double cmin = a * (k3 * er / (1.0 - er) + k2 * ea / (1.0 - ea) - k1 * ee / (1.0 - ee));
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / opt.kr() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / opt.kr() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po2(const option& opt) const
// 内服2区画モデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A・exp(-α・t)+B・exp(-β・t)-(A+B)exp(-ka・t)\n";
  cout << "%血液区画分布容積Vc=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%分布速度定数k12=" << opt.k12() << "  分布速度定数k21=" << opt.k21() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double a = 0.0, k1 = 0.0, b = 0.0, k2 = 0.0;
    double c0 = po2para(a, k1, b, k2, q, opt.f0(), opt.vd(), opt.ka(), opt.ke(), opt.k12(), opt.k21());
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double e1 = exp(-k1 * tint), e2 = exp(-k2 * tint);
    double ea = exp(-opt.ka() * tint);
    double cmin = a * e1 / (1.0 - e1) + b * e2 / (1.0 - e2) - (a + b) * ea / (1.0 - ea);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / k2 << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / k2 << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_po2lag(const option& opt) const
// 内服2区画ラグタイムモデルによる血中濃度曲線出力
{
  cout << "-s1 -dcurve,spline -w80\n";
  cout << "%";  puttitle(cout, Title);
  cout << "%モデル:" << opt.modelname() << "\n";
  cout << "%モデル関数:c(t)=A・exp(-α・t)+B・exp(-β・t)-(A+B)exp(-ka・t)\n";
  cout << "%血液区画分布容積Vc=" << opt.vd() << "  吸収率f=" << opt.f0() << "\n";
  cout << "%吸収速度定数ka=" << opt.ka() << "  排泄速度定数ke=" << opt.ke() << "\n";
  cout << "%分布速度定数k12=" << opt.k12() << "  分布速度定数k21=" << opt.k21() << "\n";
  cout << "%ラグタイムtlag=" << opt.tlag() << "\n";
  int nt = opt.nt();
  double q = 0.0;
  for (int i = 0; i < nt; ++i)
  {
    cout << "%" << i + 1 << "回目投与:時間t=" << opt.t(i) << "  用量q=" << opt.q(i) << "\n";
    q += opt.q(i);
  }
  if (nt > 1)
  {
    double tint = (opt.t(nt - 1) - opt.t(0)) / (double)(nt - 1);
    q /= (double)nt;
    double a = 0.0, k1 = 0.0, b = 0.0, k2 = 0.0;
    double c0 = po2para(a, k1, b, k2, q, opt.f0(), opt.vd(), opt.ka(), opt.ke(), opt.k12(), opt.k21());
    double cmean = c0 / opt.ke() / tint; // AUC(無限大)/tint
    double e1 = exp(-k1 * tint), e2 = exp(-k2 * tint);
    double ea = exp(-opt.ka() * tint);
    double cmin = a * e1 / (1.0 - e1) + b * e2 / (1.0 - e2) - (a + b) * ea / (1.0 - ea);
    cout << "%定常状態における平均血中濃度=" << cmean << "\n";
    cout << "%定常状態における最低血中濃度=" << cmin << "\n";
    cout << "%95%定常状態になるまでの時間=" << -log(0.05) / k2 + opt.tlag() << "\n";
    cout << "%99%定常状態になるまでの時間=" << -log(0.01) / k2 + opt.tlag() << "\n";
  }
  plot_comp(opt);
}

void compsimclass::plot_comp(const option& opt) const // 血中濃度曲線出力
{
  constexpr int RDdef = 100;    // グラフの描画範囲分割数省略時値
  double ts = 0, te = 0.0, trange = 8.0, dt = 0.0, t = 0.0;
  for (int i = 0; i < opt.nt(); ++i)
  {
    ts = opt.t(i);
    te = (i + 1 < opt.nt()) ? opt.t(i + 1) : ts + trange;
    trange = te -  ts;
    dt = trange / (double)RDdef;
    bool modelif = false, modellag = false;
    double tif = ts + opt.tif(), tlag = ts + opt.tlag();
    if (opt.model() == modeliv1if || opt.model() == modeliv2if || opt.model() == modelpo1r0 || opt.model() == modeloi1r0 || opt.model() == modeloi1r1) modelif = true;
    if (opt.model() == modelpo1lag || opt.model() == modelpo2lag) modellag = true;
    t = ts;
    for (int i = 0; i < RDdef; ++i, t += dt)
    {
      if (modelif)              // 持続投与モデル
      {
        if (t > tif)
        { cout << tif << " " << fsigmac(opt, tif) << "\n";  modelif = false; }
      }
      if (modellag)              // ラグタイムモデル
      {
        if (t > tlag)
        { cout << tlag << " " << fsigmac(opt, tlag) << "\n";  modellag = false; }
      }
      cout << t << " " << fsigmac(opt, t) << "\n";
    }
  }
  double feps = opt.q(0) / opt.vd() * 0.1 / (double)(RDdef);
  double fc = fsigmac(opt, te);
  t = te;
  do
  {
    cout << t << " " << fc << "\n";
    t += dt;
    fc = fsigmac(opt, t);
  } while ((opt.model() != modelpo1e && opt.model() != modeliv2e && fc > feps) || (opt.model() == modelpo1e && (opt.q(0) * opt.f0() * 0.01 - fc) > feps) || (opt.model() == modeliv2e && (opt.q(0) * 0.01 - fc) > feps));
  cout << t << " " << fc << endl;
}

double compsimclass::fsigmac(const option& opt, double t) const
// 累積血中濃度関数
{
  double fc = 0.0;
  for (int i = 0; i < opt.nt(); ++i) fc += fcomp(opt, t - opt.t(i), opt.q(i));
  return fc;                    // 累積血中濃度を返す
}

double compsimclass::fcomp(const option& opt, double t, double q) const
// 血中濃度関数
{
  if (t < 0.0) return 0.0;
  double qf = q * opt.f0() * 0.01;
  double c0 = qf / opt.vd();    // 初期濃度
  double a = 0.0, b = 0.0, d = 0.0, k1 = 0.0, k2 = 0.0, k3 = 0.0; // パラメータ
  double tt = 0.0;              // 計算用一時時間
  switch (opt.model())
  {
  case modeliv1if : a = q / opt.tif() / opt.vd() / opt.ke();
    if (t <= opt.tif()) return a * (1.0 - exp(-opt.ke() * t));
    tt = t - opt.tif();
    return a * (exp(-opt.ke() * tt) - exp(-opt.ke() * t));
  case modeliv2 : c0 = iv2para(a, k1, b, k2, d, q, opt.vd(), opt.ke(), opt.k12(), opt.k21());
    return a * exp(-k1 * t) + b * exp(-k2 * t);
  case modeliv22 : c0 = iv2para(a, k1, b, k2, d, q, opt.vd(), opt.ke(), opt.k12(), opt.k21());
    return d * (exp(-k2 * t) - exp(-k1 * t));
  case modeliv2T : c0 = iv2para(a, k1, b, k2, d, q, opt.vd(), opt.ke(), opt.k12(), opt.k21());
    return opt.ke() * (a * exp(-k1 * t) / k1 + b * exp(-k2 * t) / k2);
  case modeliv2e : c0 = iv2para(a, k1, b, k2, d, q, opt.vd(), opt.ke(), opt.k12(), opt.k21());
    return opt.ke() * (a * (1.0 - exp(-k1 * t)) / k1 + b * (1.0 - exp(-k2 * t)) / k2);
  case modeliv2if : c0 = iv2ifpara(a, k1, b, k2, q, opt.vd(), opt.ke(), opt.k12(), opt.k21(), opt.tif());
    if (t <= opt.tif()) return a * (1.0 - exp(-k1 * t)) + b * (1.0 - exp(-k2 * t));
    tt = t - opt.tif();
    return a * (exp(-k1 * tt) - exp(-k1 * t)) + b * (exp(-k2 * tt) - exp(-k2 * t));
  case modeloi1r0 :
  case modelpo1r0 : a = c0 / opt.tif() / (opt.ka() - opt.ke()) / opt.ke();
    if (t <= opt.tif()) return a * (opt.ka() * (1.0 - exp(-opt.ke() * t)) - opt.ke() * (1.0 - exp(-opt.ka() * t)));
    tt = t - opt.tif();
    return a * (opt.ka() * (exp(-opt.ke() * tt) - exp(-opt.ke() * t)) - opt.ke() * (exp(-opt.ka() * tt) - exp(-opt.ka() * t)));
  case modeloi1r1 : k1 = opt.ka() - opt.ke();
    k2 = opt.ke() - opt.kr();  k3 = opt.ka() - opt.kr();
    a = c0 * opt.kr() * opt.ka() / k3 / k2 / k1;
    if (t <= opt.tif()) return a * (k1 * exp(-opt.kr() * t) + k2 * exp(-opt.ka() * t) - k3 * exp(-opt.ke() * t));
    tt = t - opt.tif();
    return a * (k3 * exp(-opt.kr() * opt.tif() - opt.ke() * tt) - k2 * exp(-opt.kr() * opt.tif() - opt.ka() * tt) + k2 * exp(-opt.ka() * t) - k3 * exp(-opt.ke() * t));
  case modelpo1r1 : k1 = opt.ka() - opt.ke();
    k2 = opt.ke() - opt.kr();  k3 = opt.ka() - opt.kr();
    a = c0 * opt.kr() * opt.ka() / k3 / k2 / k1;
    return a * (k1 * exp(-opt.kr() * t) + k2 * exp(-opt.ka() * t) - k3 * exp(-opt.ke() * t));
  case modelpo1 : a = c0 * opt.ka() / (opt.ka() - opt.ke());
    return a * (exp(-opt.ke() * t) - exp(-opt.ka() * t));
  case modelpo1T : a = qf * opt.ka() * opt.ke() / (opt.ka() - opt.ke());
    return a * (exp(-opt.ke() * t) / opt.ke() - exp(-opt.ka() * t) / opt.ka());
  case modelpo1e : a = qf * opt.ka() * opt.ke() / (opt.ka() - opt.ke());
    return a * ((1.0 - exp(-opt.ke() * t)) / opt.ke() - (1.0 - exp(-opt.ka() * t)) / opt.ka());
  case modelpo1lag : if (t <= opt.tlag()) return 0.0;
    tt = t - opt.tlag();
    a = c0 * opt.ka() / (opt.ka() - opt.ke());
    return a * (exp(-opt.ke() * tt) - exp(-opt.ka() * tt));
  case modelpo2 : c0 = po2para(a, k1, b, k2, q, opt.f0(), opt.vd(), opt.ka(), opt.ke(), opt.k12(), opt.k21());
    return a * exp(-k1 * t) + b * exp(-k2 * t) - (a + b) * exp(-opt.ka() * t);
  case modelpo2lag : if (t <= opt.tlag()) return 0.0;
    tt = t - opt.tlag();
    c0 = po2para(a, k1, b, k2, q, opt.f0(), opt.vd(), opt.ka(), opt.ke(), opt.k12(), opt.k21());
    return a * exp(-k1 * tt) + b * exp(-k2 * tt) - (a + b) * exp(-opt.ka() * tt);
  case modeliv1 :
  default : a = q / opt.vd();
    return a * exp(-opt.ke() * t);
  }
}

double compsimclass::iv2para(double& a, double& k1, double& b, double& k2, double& d, double q, double vd, double ke, double k12, double k21) const
// 静注2区画モデルのパラメーターを返す
{
  double c0 = q / vd;
  double al = k12 + k21 + ke;
  double be = k21 * ke;
  k1 = (al + sqrt(al * al - 4.0 * be)) * 0.5;
  k2 = (al - sqrt(al * al - 4.0 * be)) * 0.5;
  a = c0 * (k1 - k21) / (k1 - k2);
  b = c0 * (k21 - k2) / (k1 - k2);
  d = c0 * k12 / (k1 - k2);
  return c0;                    // 初期濃度を返す
}

double compsimclass::iv2ifpara(double& a, double& k1, double& b, double& k2, double q, double vd, double ke, double k12, double k21, double tif) const
// 静注2区画持続注入モデルのパラメーターを返す
{
  double c0 = q / vd;
  double k0 = q / tif;
  double al = k12 + k21 + ke;
  double be = k21 * ke;
  k1 = (al + sqrt(al * al - 4.0 * be)) * 0.5;
  k2 = (al - sqrt(al * al - 4.0 * be)) * 0.5;
  a = k0 / vd * (k1 - k21) / (k1 - k2) / k1;
  b = k0 / vd * (k21 - k2) / (k1 - k2) / k2;
  return c0;                    // 初期濃度を返す
}

double compsimclass::po2para(double& a, double& k1, double& b, double& k2, double q, double f0, double vd, double ka, double ke, double k12, double k21) const
// 内服2区画モデルのパラメーターを返す
{
  double c0 = q * f0 * 0.01 / vd;
  double al = k12 + k21 + ke;
  double be = k21 * ke;
  k1 = (al + sqrt(al * al - 4.0 * be)) * 0.5;
  k2 = (al - sqrt(al * al - 4.0 * be)) * 0.5;
  a = c0 * ka * (k1 - k21) / (ka - k1) / (k1 - k2);
  b = c0 * ka * (k21 - k2) / (ka - k2) / (k1 - k2);
  return c0;                    // 初期濃度を返す
}

void compsimclass::put(option& opt)
// 区画モデルによる血中濃度曲線を出力する
{
  if (opt.check() != 0) return;
  switch (opt.model())
  {
  case modeliv1if : plot_iv1if(opt);
    break;
  case modeliv22 :
  case modeliv2T :
  case modeliv2e :
  case modeliv2 : plot_iv2(opt);
    break;
  case modeliv2if : plot_iv2if(opt);
    break;
  case modeloi1r0 : plot_oi1r0(opt);
    break;
  case modeloi1r1 : plot_oi1r1(opt);
    break;
  case modelpo1T :
  case modelpo1e :
  case modelpo1 : plot_po1(opt);
    break;
  case modelpo1lag : plot_po1lag(opt);
    break;
  case modelpo1r0 : plot_po1r0(opt);
    break;
  case modelpo1r1 : plot_po1r1(opt);
    break;
  case modelpo2 : plot_po2(opt);
    break;
  case modelpo2lag : plot_po2lag(opt);
    break;
  case modeliv1 :
  default : plot_iv1(opt);
  }
}

int main(int argc, char* argv[]) // メインルメーチン
{
  compsimclass compsim;         // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    for (int i = 1; i < argc; ++i)
    {
      if (*(argv[i]) == CHRopt) opt.set(argv[i] + 1);
      else throw runtime_error(string("オプション指定文字が違います→") + string(argv[i]));
    }
    compsim.put(opt);
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(opt.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(opt.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
