Under-Construction 工作>ダイビングコンピュータを作る>シミュレータの作成

ダイビングコンピュータを作る

最終更新 (2002/05/19)

この節では、前節での理論を元に、パソコンで動作するダイビングコンピュータのシミュレータのプログラムを解説します。

  1. 基本部品
  2. 計算部分
  3. 使い方
  4. ダウンロード

基本部品

データ構造

このシミュレータに必要なデータ構造として、次の2つを定義します。

コンパートメント構造体
この構造体で、コンパートメント一つを表します。

メンバとしては、そのコンパートメントごとの定数を表すもの、そのコンパートメントの、その時点での状態を保持するもの、の2通りあります。
前者としてk値(ハーフタイムτに体して k=ln(1/2)/τ)、M値(M0 と ΔMの2つ)、後者として、現時点での窒素分圧、その組織のシーリング、フロアー、NDLがあります。シーリング、フロアー、NDL はすべての組織の中の最弱のものが、体全体の値として、最終的に採用されます。。

人体全体としては、コンパートメントの集合体なので、この構造体の配列として表すことができます。また、各組織の定数は、最初に読み込まれた時から変化しませんが、コンパートメントの窒素分圧は、各時刻での窒素の溶解と排出によって、再計算去れます。
コンパートメントの数は、配列として、最大100個まで確保しておきますが、設定ファイル(config.dat)を読み込んだ時に、コンパートメント数は数えられ、NUM_COMPARTMENT にセットされます。

struct compartment{
  double k;       /* 各組織の k値 = ln(1/2)/tau */
  double M_a;      /* 各組織の M値 */
  double delta_M;  /* 各組織の M値 */

  double n2;      /* 各組織の現在の窒素分圧 */
  double seal;
  double floor;
  double ndl;
};
int NUM_COMPARTMENT = 9; /* コンパートメントの数 */
#define MAX_NUM_COMPARTMENT 100
struct compartment body[MAX_NUM_COMPARTMENT];
     
深度構造体
この構造体で、その時刻での、周囲環境を表します。

メンバとしては、深度、呼気中の窒素の割合、個忌中のヘリウムの割合です。ヘリウムについては、今は考えていませんが、将来のオプションです。もちろん、他にもいろいろと混合気体をデータとしてもたしてもかまいません。

この深度構造体は、実際のダイビングコンピュータとして利用する時は、配列は必要なくて、現時点でのデータが一つだけあって、それを随時更新していく形でOKです。しかし、シミュレータ上では、時刻を引数とした配列になります。深度データファイル(depth.dat)で初期化されます。今回はメモリの省略のため、1分ごとの配列として用意しました。計算は、1秒おきに行うことを前提としてありますが。

struct depth_pattern{
  double d;       /* 水圧 */
  double n2;      /* 呼気中の窒素の割合 0.0〜1.0       */
  double he;      /* 呼気中のヘリウムの割合 0.0〜1.0 */
};
struct depth_pattern depth[MAX_TIMESPAN];
     

設定ファイルの読み込み

このシミュレータに必要なデータは、1. 各コンパートメントのハーフタイムと M値、2. 深度と混合ガスの時間変化を表したもの、の2つです。便利なように、ファイルから読み込むようにしましょう。

コンパートメントの設定ファイルのフォーマットとしては、次のものとします。コンパートメントの最大数は、プログラム中の #define MAX_NUM_COMPARTMENT の数までで、任意の数を取れるようにします。

# 1カラム目が # の行はコメント
# データ行は、任意数のスペースで区切られる 3列のデータ
# 先頭から、ハーフタイム M_0値 ΔM値
5   31.7 1.8
10  26.8 1.6
20  21.9 1.5
40  17.0 1.4
80  16.4 1.3
120 15.8 1.2
160 15.5 1.16
200 15.5 1.1
240 15.2 1.1
# コンパートメントの数は、プログラム中の MAX_NUM_COMPARTMENT まで

これを読み込む関数は、次の通り。プログラム上注意するのは、

  1. 良い子は gets ではなく fgets を使う。当然 scanf で直接読み取るのも不可。
  2. sscanfの format の文法に注意、空白は「任意個数の空白を表します。あと %lf と %f も注意といえばそうかな。
void read_config(char *filename, struct compartment *body)
{
  FILE *file;
  int i;
  char buf[LINE_BUF];
  double tau, M_a, delta_M;

  
  if(file = fopen(filename, "rt")){
    i = 0;
    while(fgets(buf, LINE_BUF, file)){
      if (buf[0] == '#'){ /* 先頭が # はコメント行 */
	;
      } else {
	sscanf(buf, " %lf %lf %lf", &tau, &M_a, &delta_M);
	body[i].k  = tau2k(tau)/60.0;   /* τは分で記述される */
	body[i].M_a  = M_a;           /* M はmsw  */
	body[i].delta_M  = delta_M;
	i++;
	if (i > MAX_NUM_COMPARTMENT)  return;
      }
    }
    NUM_COMPARTMENT = i;
  } else {
    fprintf(stderr, "Can't load config file %s\n", filename);
    exit(-1);
  }
}

深度のほうも、似たような感じです。ただし、データファイルは 1分ごとに用意します。すこし作成がメンドウですね。単位系に注意です。深度は絶対圧で、それ以外は比率で書きます。次の例は 30m に潜った例です。水面の圧力 10[msw] をくわえていることに注意。

# 深度データファイル
# 深度(絶対圧) 呼気中の窒素の比率 呼気中のヘリウムの比率
10.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
# 10分
40.0 0.8 0.0
40.0 0.8 0.0
40.0 0.8 0.0
....略

このデータファイルを読み取るプログラムは read_config とすごい似てるので省略します。

計算部分

窒素の溶解と排出

→理論編はこちら

体内窒素量は、各コンパートメント毎に、現在の体内窒素分圧、ハーフタイムτ、呼気窒素分圧より計算できます。(よくわからないですが)呼気全体圧から水の飽和蒸気圧を引くらしいです。body[].k は、ハーフタイムから、あらかじめ、一秒ごとの増分になるように変換しておきます。単位系は要注意ですが、全体を通して絶対圧での水深[msw]であつかいます。

いずれの計算関数も、計算対象は、一つのコンパートメントです。第一引数に、コンパートメント構造体のポインタを与え、第二引数には深度構造体のポインタを与えます。計算結果は、コンパートメント構造体を書き換えることによって返します

計算自体は理論編そのままの簡単なものです。

/* 各コンパートメントの現在の状態と、現在の呼気窒素分圧から
   次の時点の各コンパートメントの状態を算出します */
void calc(struct compartment *comp,      /* コンパートメント */
	  struct depth_pattern *depth)   /* 深度             */
{
  comp->n2 += ((depth->d - PH2O) * depth->n2 - comp->n2) * comp->k;
}

NDL計算

→理論編はこちら

void calc_ndl(struct compartment *comp,       /* コンパートメント */
	      struct depth_pattern *depth)    /* 深度             */
{
  double pamb;
  double po;
  double pi;
  double mo;
  double det;

  /* http://www.abysmal.com/pages/articles/calculating_no_stop_time.html */
  pamb = depth->d;
  pi   = (pamb - PH2O)*depth->n2;
  po   = comp->n2 - PH2O;
  mo   = comp->M_a;

  det = (pi - mo) / (pi - po);

  if (det < 0.0){                         /* ln の中が負→ ndl が複素数 */
    if ((mo > pi) && (pi > po)){          /* a. 何分潜っても無限圧 */
      comp->ndl = 59.0;
    } else if ((mo < pi) && (pi < po)){   /* b. */
      comp->ndl = -9.0;
    }
  } else if (det > 1.0){                  /* ln の中が1以上→ ndl が負 */
    if ((mo < pi) && (pi > po)){          /* c. 減圧潜水に突入 */
      comp->ndl = (-1.0/(comp->k * 60.0)) * log(det);
    } else if ((mo > pi) && (pi < po)){   /* d. 減圧終了、浮上可能 */
      comp->ndl = 55.0;
    }
  } else {                                /* ln の中が1以下→ ndl が正 */
    if ((pi > mo) && (mo > po)){   /* A. 通常の潜水(窒素蓄積) */
      comp->ndl = (-1.0/(comp->k * 60.0)) * log(det);
    } else if ((pi < mo) && (mo < po)){   /* B. 通常の潜水(窒素排出) */
      comp->ndl = (-1.0/(comp->k * 60.0)) * log(det);
    } else {                              /* ここには来ない */
      printf("mo = %lf, pi = %lf, po = %lf det = %lf\n", mo, pi, po, det);
      assert(0);
    }
  }
}

シーリング

→理論編はこちら

/* 各コンパートメントの窒素分圧からM値を引いた浮上可能水深のうち
   最も深いものを返します */
void calc_seal(struct compartment *comp,      /* コンパートメント */
		    struct depth_pattern *depth) /* 深度             */
{
  comp->seal = (comp->n2 - comp->M_a) / comp->delta_M + 10.0;
}

フロアー

→理論編はこちら

/* 各コンパートメントの体内窒素分圧と呼気窒素分圧から、
   窒素放出が始まる水深(絶対圧換算)を返します */
void calc_floor(struct compartment *comp,      /* コンパートメント */
		  struct depth_pattern *depth) /* 深度             */
{
  comp->floor = comp->n2 / depth->n2;
}  

最弱値の選定

計算関数は、各コンパートメントごとに計算してあるので、その中の最弱値を、体全体の値として選定して、見易く表示する必要があります。

    /* 最大値の決定 */
    seal_max  = -10000.0;
    floor_max = -10000.0;
    ndl_min   =  10000.0;
    for(j = 0; j < NUM_COMPARTMENT; j++){
      if (seal_max  < body[j].seal) { seal_max  = body[j].seal;  }
      if (floor_max < body[j].floor){ floor_max = body[j].floor; }
      if (ndl_min   > body[j].ndl)  { ndl_min   = body[j].ndl;   }
    }

使い方

Makefile と algo.c を用意して、

$ make

の make 一発でプログラムは作成されます(``$'' はシェルのプロンプト)。実行方法は、第一引数に、M値/ハーフタイム設定ファイルを、第二引数に深度データファイルを設定して実行すると、標準出力に計算結果を出力します。それを、result というファイルに取ると、plot*.gp ファイル群によって gnuplot でグラフ化できます。

$ ./algo config.dat depth.dat > result
$ ./plotA.gp

plotA.gp は、体全体の様子を表示します。plot?.gp(?は1〜9)は、?番目の組織について表示します。いずれも、ターミナルでなにか入力するとグラフ表示は終了します。

ダウンロード

algo.c プログラムのソースファイルです。
Makefile
plot.gp GNUPLOT 用のプロットコマンドファイル。実行可能。
plot1.gp
plot2.gp
plot3.gp
plot4.gp
plot5.gp
plot6.gp
plot7.gp
plot8.gp
plot9.gp
plotA.gp
config.dat M値、ハーフタイムの設定ファイル例
depth_100m.dat 100m潜水の深度データファイル(ハーフタイムチェック用)
depth_30m.dat 減圧のチェック用
depth_30m_e.dat ナイトロックスを利用した30m潜水の計算例
depth_30m_edeco.dat 減圧ミックスを利用した30m潜水の計算例

注意

募集



近藤靖浩