// Jacobi.v 差分法により熱平衡状態を求める

// まわりの温度(固定値)を定義
#define T_SIDE 0.
#define T_UD   200.

// 温度差の閾値
#define epsilon 1.

module Jacobi(in:input,out:output)
float in on pixel cache 1;
float out on pixel;
{
    // 全体の温度差のうち最も値の大きなもの
    float res;
#vios mutex res off max
    res = 2.;
    // 画像を初期化します
    parallel{
	in[][] = 0.;
    }

    // 新たに求められた温度と過去の温度の差のうち
    // 最も大きなものが閾値以下になるまで繰り返す
    while(epsilon < res){
	// 最大の温度差を初期化
	res = 0.;

	parallel{
	    float left=0.,right=0.,up=0.,down=0.;
	    // 注目点における温度差
	    float localres = 0.;

	    // 上下左右に隣接する点の温度を取得
	    // 左
	    if(hotx(in) == 0)        left = T_SIDE;
	    else                     left = in[-1][];

	    // 右
	    if(hotx(in) == width(in)-1)  right = T_SIDE;
	    else                         right = in[+1][];

	    // 上
	    if(hoty(in) == 0)        up = T_UD;
	    else                     up = in[][-1];

	    // 下
	    if(hoty(in) == height(in)-1) down = T_UD;
	    else                         down = in[][+1];

	    // 上下左右の温度の平均が注目点の温度となる
	    out[][] = (left + right + up + down)/4;
	    // 注目点の以前の温度との差をとる
	    localres = fabs(in[][] - out[][]);
	    // 差型の点と比べ大きければその値で更新
	    (res < localres) && (res = localres);
	}
	// 出力画像に結果書き込み
	parallel{
	     in[][] = out[][];
	}
	// キャッシュ値を更新
	sync_cache(in);
    }
}