この節では、前節での理論を元に、パソコンで動作するダイビングコンピュータのシミュレータのプログラムを解説します。
このシミュレータに必要なデータ構造として、次の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 まで
これを読み込む関数は、次の通り。プログラム上注意するのは、
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;
}
→理論編はこちら
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潜水の計算例 |