雑感等

音楽,数学,語学,その他に関するメモを記す.

扇型に内接する楕円

扇型に内接する楕円を求めた

――かったが、実質的には二等辺三角形に内接する楕円の式が求まった。

扇型は、半径を1、中心角を \theta とする。

楕円は、一方(扇型の径方向)の半径を a 、もう一方の半径を b とする。

各図形は図の座標に配置し、 \theta  a を既知の値とすると、 b は次式で与えられる。

 \displaystyle
b = \frac{ \frac{a^2}{a-1} -a+1}{\tan {\left( \frac{ \pi - \theta}{2} \right) } \sqrt{\frac{1-2a}{(a-1)^2}} }

この b を使い、上図の楕円の式は次式で与えられる。

 \displaystyle
\frac{x^2}{b^2} + \frac{(y+1-a)^2}{a^2} = 1

このページで式を動かせる:二等辺三角形に内接する楕円 – GeoGebra

丸かっこで囲まれた部分以外にマッチする正規表現

丸かっこで囲まれた部分以外にマッチする正規表現

[^\(\)]+(?=\([^\(\)]+\))|[^\(\)]+$|^[^\(\)]+

テキストエディタMeryで動作確認済み。

説明

この正規表現は、orを意味するパイプ記号"|"と、3つの部分(下記ABC)で構成されている。

A | B | C

 ↓

[^\(\)]+(?=\([^\(\)]+\)) | [^\(\)]+$ | ^[^\(\)]+

つまり、「AまたはBまたはC」のパターンを表す。

まず、比較的単純なBから見る:[^\(\)]+$

これは、さらに3つの要素に分けられる。

D E F

 ↓

[^\(\)] + $

最初の角括弧[]は、「指定した文字を含まない」という意味の表現:[^~~~~]だ。

"[^"と"]"に挟まれた各文字以外の文字にマッチする。ここでいう各文字は \( と \)だ。エスケープされているが要は丸括弧"("または")"を意味する。

つまり、[^\(\)]は「丸括弧でない何らかの1文字」を意味する。

次、+は「前の文字が1文字以上連続する」を意味する。

次、$は、「行末」を意味する。

ここまで合わせて考えると[^\(\)]+$は「丸括弧以外の文字が1文字以上連続し、その直後に行末が来る」パターンを表す。

簡単に言うと、「丸括弧から行末までの文字 ※丸括弧は含まず」にマッチする。

例えば入力が「abc0123(hijk」なら、「hijk」にマッチする。

次にC^[^\(\)]+を見る。

最初の^は、「行頭」を表す。※"[^"の"^"とは意味が異なる。

次は上で見た、[^\(\)]+と全く同じだ。

つまり、^[^\(\)]+は「行頭の直後に、丸括弧以外の文字が1文字以上連続する」パターンを表す。

簡単に言うと、「行頭から丸括弧までの文字 ※丸括弧は含まず」にマッチする。

例えば入力が「abc0123(hijk」なら、「abc0123」にマッチする。

最後にA[^\(\)]+(?=\([^\(\)]+\))を見る。

まず、大枠は

G (?= H )

[^\(\)]+ (?= \([^\(\)]+\) )

G (?= H ) のパターンは、肯定先読みと呼ばれる。

意味は、「G。ただし、Hが直後にあるもの」のパターンを表す。

「マッチしてほしいのはGだけど、Gならなんでもいいわけじゃなくて、Hの直前にあるGだけにマッチしてほしい」という事。

さて、GHを見ていく。

Gは上で見た[^\(\)]+と全く同じだ。 つまり、「丸括弧以外の文字が1文字以上連続する」パターンを表す。

次、H\([^\(\)]+\)は、さらに3つの要素に分解できて、 \( [^\(\)]+ \)となる。

\(は開き括弧、 [^\(\)]+は「丸括弧以外の文字が1文字以上連続する」パターン、 \)は閉じ括弧。

Hはつまり、「丸括弧で囲まれた文字列(ただし囲んでいる丸括弧も含む)」パターンを表す。

つまり[^\(\)]+ (?= \([^\(\)]+\) )は、 「丸括弧で囲まれた文字列の直前にある、丸括弧以外の文字列」のパターンを表す。

例えば入力が「)abcd(efgh)0123(4567)ijkl」なら、「abcd」と「0123」にマッチする。

参照

https://ja.stackoverflow.com/questions/9024/%E3%81%8B%E3%81%A3%E3%81%93%E3%81%A7%E5%9B%B2%E3%81%BE%E3%82%8C%E3%81%9F%E6%96%87%E5%AD%97%E4%BB%A5%E5%A4%96%E3%82%92%E3%83%9E%E3%83%83%E3%83%81%E3%81%95%E3%81%9B%E3%81%9F%E3%81%84

https://qiita.com/shotets/items/98f3828b6e5f08d42498

物理モデルで自励振動させようとした

いろいろ調整してみたら自励振動させられてないっぽい。管体部分でループ(ハウリングみたいな状態)していただけだった。

物理モデリング音源の特許(だいたい権利失効してる)の実施例にある、励振部をpythonで実装した。

一応音は出るが、パラメータの調整が甘くか弱い音しか鳴らない。

とりあえず音が出る状態でメモしておく。

pythonソース

import pyaudio
import numpy as np
from collections import deque
import matplotlib.pyplot as plt
import scipy.signal as signal

fs_sound = 44100  # [1/s] sampling freq.
fs_internal = 44100
c = 340  # [m/s] speed of sound


class DelayLine:
    def __init__(self, n_delay):
        n_delay = np.round(n_delay).astype(np.int64)
        self._dl = deque([0] * n_delay)
        pass

    def update(self, x_in):
        self._dl.append(x_in)
        x_out = self._dl.popleft()
        return x_out

    pass


class Filter(DelayLine):
    def __init__(self, coef_fir):
        n_taps = coef_fir.size
        super().__init__(n_taps)
        self.coef_fir = coef_fir[::-1]
        pass

    def update(self, x_in):
        self._dl.append(x_in)
        self._dl.popleft()
        x_out = self._dl @ self.coef_fir

        return x_out

    pass


def sound_play(samples):
    p = pyaudio.PyAudio()
    stream = p.open(format=pyaudio.paFloat32,
                    channels=1,
                    rate=fs_sound,
                    frames_per_buffer=1024,
                    output=True)
    t_end = 1
    n_samples = t_end * fs_sound
    n_t = np.arange(0, fs_sound)
    f = 440

    if np.max(np.abs(samples)) > 1:
        samples = samples / np.max(np.abs(samples)) * 0.9
        pass

    stream.write(samples.astype(np.float32).tobytes())
    stream.close()

    pass


def filter_analy(a, fs):
    freq, res = signal.freqz(a)
    tate = 3
    yoko = 1
    plt.subplot(tate, yoko, 1)
    plt.plot(0.5 * fs * freq / np.pi, np.abs(res))
    plt.subplot(tate, yoko, 2)
    plt.plot(0.5 * fs * freq / np.pi, np.angle(res))
    plt.subplot(tate, yoko, 3)
    plt.plot(a[::-1], ".-")
    plt.show()
    pass


def main():
    # FILTER SETTING ###########################################
    fs = fs_internal
    fc = 2000
    freqs = np.array([0, fc, fc + 1000, fs / 2])
    a = signal.firwin2(numtaps=500, freq=freqs / (fs / 2), gain=[1, 1, 0, 0], nfreqs=fs)
    a = signal.minimum_phase(a, "homomorphic")

    fl1 = Filter(a)
    fl_press2reed = Filter(a)

    fc = 30
    freqs = np.array([0, fc, fc + 50, fs / 2])
    a = signal.firwin2(numtaps=2001, freq=freqs / (fs / 2), gain=[0, 0, 1, 1], nfreqs=fs)
    a = signal.minimum_phase(a, "homomorphic")
    fl_outhpf = Filter(a)
    ############################################
    dl_len = 0.8  # [m] delay line length
    n_dl = fs_internal * dl_len / c
    print(n_dl)
    # filter_analy(a, fs)

    dl1 = DelayLine(n_dl)
    dl2 = DelayLine(n_dl)
    dl_s = DelayLine(1)

    t_simulation = 2
    n_sim = fs_internal * t_simulation

    xs_input = np.zeros(n_sim)

    in_max = 3
    in_min = 0.15
    n_on1 = 10000
    n_on2 = 20000
    n_off1 = 30000
    n_off2 = 80000
    # xs_input[0] = 1
    for n in range(n_sim):
        if n_on1 <= n and n < n_on2:
            x = in_min + (in_max - in_min) / (n_on2 - n_on1) * (n - n_on1)
            pass
        elif n_on2 <= n and n < n_off1:
            x = in_max
            pass
        elif n_off1 <= n and n < n_off2:
            x = in_min + (in_max - in_min) * (1 - 1 / (n_off2 - n_off1) * (n - n_off1))
            pass
        else:
            x = in_min
            pass
        xs_input[n] = x
        pass

    xs_output = np.zeros(n_sim)
    xs_reed = np.zeros(n_sim)
    xs_volvel = np.zeros(n_sim)

    k_refl_reed = 0.1

    # INIT ############################
    do1 = 0.0
    do2 = 0.0
    do_s = 0.0

    def press2reed(x):
        out = x
        if x < -1:
            out = -1
            pass
        elif x > 1:
            out = 1
            pass
        else:
            out = x
            pass
        # return 1 / (1 + np.exp(-x))
        return out
        pass

    def press2vel(x):
        # return np.abs(x)
        return x
        pass

    def volvel2press(x):
        return x
        pass

    for n in range(n_sim):
        x_n = xs_input[n]
        s_press = do_s - x_n
        s_press_reed = s_press - 0.8
        s_reed_area = press2reed(fl_press2reed.update(s_press_reed))
        # s_reed_area = press2reed(s_press_reed)
        xs_reed[n] = s_reed_area

        s_velocity = press2vel(s_press)
        s_volvel = s_reed_area * s_velocity
        s_reedout = volvel2press(s_volvel)

        xs_volvel[n] = s_volvel

        s_forward_d = (1 + k_refl_reed) * s_reedout + (-k_refl_reed) * do2
        do1 = dl1.update(s_forward_d)
        x_n = fl_outhpf.update(do1)

        s_backward_1 = fl1.update(x_n)
        do2 = dl2.update(s_backward_1)
        s_backward_2 = (1 - k_refl_reed) * do2 + k_refl_reed * s_reedout

        do_s = dl_s.update(s_backward_2)

        xs_output[n] = x_n
        pass

    plt.plot(xs_output, label="out")
    plt.plot(xs_input, "--", label="in")
    plt.plot(xs_reed, label="reed")
    plt.plot(xs_volvel, label="volvel")
    plt.legend()
    plt.show()
    sound_play(xs_output)

    pass


if __name__ == '__main__':
    main()

Rustで配列(Vec)の要素数を変数で指定:vec!マクロ

出典

https://users.rust-lang.org/t/attempt-to-use-a-non-constant-value-in-a-constant/32112#:~:text=Sep%20%2719-,Use%20vec!%5BDefault%3A%3Adefault()%3B%20n%5D%20or%20Vec%3A%3Awith_capacity(n)%20and%20push%20to%20it.,-Solution

#

配列の要素数コンパイル時に決定されているため、入力に応じて長さを伸縮できない。

下記はコンパイルエラーになる。(要素数arr_lenが変数のためコンパイルエラー)

let arr_len: usize = 8;
let arr = [0 as u8; arr_len];

一方、Vecだとコンパイル時に要素数が決定されていなくても良い。

下記のようにvec!マクロで宣言すると楽。(要素数はvec_len)

let vec_array = vec![0 as u8; vec_len];

#

下記rustソースでは、

  • 標準入力から文字列を受け取る。
  • 文字列から数値(vec_len)に変換する。
  • vec_len個の要素を持つVecを生成する(vec_array)。
  • vec_arrayの要素数と全要素を表示する。

標準入力に応じてvec_arrayの要素数が変化するのが確認でき、 実行時に(入力に応じて)要素数を決定できる。

use std::io;

fn main() {
    //文字列格納先
    let mut input_str = String::new();

    //標準入力から文字列を読む
    println!("array len:");
    io::stdin()
        .read_line(&mut input_str)
        .expect("failed to read line");

    //usize型に変換
    let vec_len: usize = input_str.trim().parse().expect("failed to parse");

    //u8型を要素に持つVec配列を生成
    let vec_array = vec![0 as u8; vec_len];

    //生成した配列長と、配列を表示
    println!("array len={}\tarray={:?}", vec_array.len(), vec_array);
}

実行例(10を入力)

array len:
10
array len=10    array=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Process finished with exit code 0

C言語で動的配列と調べるとmallocがすぐに出るが、 rustだと検索してもvec!が直球で出てこなかった気がした。(ドキュメントをちゃんと読めば書いてあると思うが)

複数フォルダをそれぞれzipに圧縮(バッチファイル)

↓このあたりのソースを改変して、圧縮後のzipファイルの保存先を指定できるようにした。

rezv.net

oshiete-suitman.com

batプログラムソースは以下。

【保存先フォルダ(フルパス)】のところに、保存先フォルダのパスを書く。

ZIP_PATHには7zipの実行ファイルのパスを指定する。

@echo off
set ZIP_PATH="C:\Program Files\7-Zip\7z.exe"
for %%f in (%*) do (
call :filename %%f
)
exit /b

:filename
%ZIP_PATH% a -tzip "【保存先フォルダ(フルパス)】\%~n1.zip" "%~1"

%%fから、圧縮するフォルダ名だけ切り出す方法がわからなかったので、

サブルーチン:filenameを呼び出して、%~n1でフォルダ名だけ切り出すようにした。

二次元の波動方程式を陰解法+ガウス・ザイデル法で解く(Rust)

バグを直した&メッシュ数に応じて出力画像の比率を変えるようにした。

■出力画像

■Rustソース

use std::time::{Duration, Instant};

use plotters::prelude::*;

const PLT_RATIO: usize = 20;
const T_MAX: usize = 1000;
const X_MAX: usize = 10;
const Y_MAX: usize = 31;

const M_ELEM: usize = 3 * 4 + 4 * 2 * (X_MAX - 2) + 4 * 2 * (Y_MAX - 2) + 5 * (X_MAX - 2) * (Y_MAX - 2);


const T_WINDOW: usize = 4;

struct SparseMatrix {
    col: [usize; M_ELEM],
    row: [usize; M_ELEM],
    val: [f64; M_ELEM],
}

fn sub_main() {
    println!("M_ELEM = {:?}", M_ELEM);
    // plot ////////////////////////////////////////////////////////////////////////////////////////

    let root = BitMapBackend::gif("testrustimage.gif", ((X_MAX * PLT_RATIO) as u32, (Y_MAX * PLT_RATIO) as u32), 50).unwrap()
        .into_drawing_area();

    let cells = root.split_evenly((Y_MAX, X_MAX));

    let mut cell_nx;
    let mut cell_ny;
    let mut val;
    let mut col;


    // calc ////////////////////////////////////////////////////////////////////////////////////////
    const DT: f64 = 1.0 / 44100.0;
    const DX: f64 = 0.05;
    const C: f64 = 347.0;
    const COURANT: f64 = C * DT / DX;

    println!("DT:{},DX:{},C:{},COURANT:{}", DT, DX, C, COURANT);

    //u[time][y][x]
    let mut u: [[[f64; X_MAX]; Y_MAX]; T_WINDOW] = [[[0.0f64; X_MAX]; Y_MAX]; T_WINDOW];

    let mut M: [[f64; X_MAX * Y_MAX]; X_MAX * Y_MAX] = [[0.0f64; X_MAX * Y_MAX]; X_MAX * Y_MAX];
    let mut x: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX];
    let mut b: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX];

    let r = 0.5;
    let nu2 = COURANT * COURANT;


    //define M
    {
        let diag = 4.0 * r * nu2 + 2.0;
        let oth_elem = -r * nu2;
        let mut curr_ind = 0;
        for ny in 0..Y_MAX {
            for nx in 0..X_MAX {
                curr_ind = ny * X_MAX + nx;

                M[curr_ind][curr_ind] = diag;
                if nx > 0 {
                    M[curr_ind][curr_ind - 1] = oth_elem;
                }
                if nx < X_MAX - 1 {
                    M[curr_ind][curr_ind + 1] = oth_elem;
                }
                if ny > 0 {
                    M[curr_ind][curr_ind - X_MAX] = oth_elem;
                }
                if ny < Y_MAX - 1 {
                    M[curr_ind][curr_ind + X_MAX] = oth_elem;
                }
            }
        }
    }
    let M = M;

    println!("u:{}", u[0][0][X_MAX - 1]);

    let mut u_max = 2.0;
    let mut u_max_t: usize = 0;
    let mut u_max_x: usize = 0;
    let mut u_max_y: usize = 0;
    let mut u_min = -1.0;
    let mut u_min_t: usize = 0;
    let mut u_min_x: usize = 0;
    let mut u_min_y: usize = 0;

    for nx in 1..X_MAX - 1 {
        u[1][2][nx] = 1.0;
    }


    let mut t_curr;
    let mut t_prev;
    let mut t_next;

    let loop_start = Instant::now();

    //main loop
    //gauss seidel: solve Mx=b
    {
        let omega = 1.0;
        let mut x_temp = 0.0;
        let mut m_temp = 0.0;
        let mut gs_diff = 0.0;

        let mut curr_ind = 0;
        for nt in 1..T_MAX - 1 {
            t_curr = nt % T_WINDOW;
            t_prev = (nt - 1) % T_WINDOW;
            t_next = (nt + 1) % T_WINDOW;

            //define b
            for ny in 1..=Y_MAX - 2 {
                for nx in 1..=X_MAX - 2 {
                    curr_ind = ny * X_MAX + nx;
                    b[curr_ind] =
                        r * nu2 * (u[t_prev][ny][nx - 1] + u[t_prev][ny - 1][nx] + u[t_prev][ny][nx + 1] + u[t_prev][ny + 1][nx])
                            + 2.0 * (1.0 - r) * nu2 * (u[t_curr][ny][nx - 1] + u[t_curr][ny - 1][nx] + u[t_curr][ny][nx + 1] + u[t_curr][ny + 1][nx])
                            - (4.0 * r * nu2 + 2.0) * u[t_prev][ny][nx] + (8.0 * (r - 1.0) * nu2 + 4.0) * u[t_curr][ny][nx];
                }
            }
            //gauss seidel
            for lp in 0_usize..200_usize {
                gs_diff = 0.0;
                for ny in 0..X_MAX * Y_MAX {
                    for nx in 0..ny {
                        m_temp -= M[ny][nx] * x[nx];
                    }
                    for nx in (ny + 1)..X_MAX * Y_MAX {
                        m_temp -= M[ny][nx] * x[nx];
                    }
                    x_temp = (b[ny] + m_temp) / M[ny][ny];
                    gs_diff += (x_temp - x[ny]).abs();
                    // println!("{:?}:{:?}:{:?}", nt, lp, x_temp - x[ny]);
                    x[ny] = omega * x_temp + (1.0 - omega) * x[ny];

                    m_temp = 0.0;
                }
                if gs_diff < f64::MIN_POSITIVE {
                    println!("break {}:{}", nt, lp);
                    break;
                }
            }
            println!("{:?}:{:?}", nt, gs_diff);
            //x->u
            for ny in 0..Y_MAX {
                for nx in 0..X_MAX {
                    curr_ind = ny * X_MAX + nx;
                    u[t_next][ny][nx] = x[curr_ind];
                }
            }
            // plot ////////////////////////////////////////////////////////////////////////////
            for (area, num) in (&cells).into_iter().zip(0..) {
                cell_nx = num % X_MAX;
                cell_ny = num / X_MAX;
                val = u[t_curr][cell_ny][cell_nx];
                col = (val - u_min) / (u_max - u_min);
                area.fill(&RGBColor((col * 255.0) as u8, (col * 255.0) as u8, (col * 255.0) as u8)).unwrap();
            }
            root.present().expect("unable to write");
        }
    }
    let duration = loop_start.elapsed();
    println!("elapsed:{:?}", duration);


    println!("max u[{}][{}][{}]={}", u_max_t, u_max_y, u_max_x, u_max);
    println!("min u[{}][{}][{}]={}", u_min_t, u_min_y, u_min_x, u_min);
}


fn main() {
    {
        const STACK_SIZE: usize = 1024 * 1024 * 1024 * 1024;
        std::thread::Builder::new()
            .stack_size(STACK_SIZE)
            .spawn(sub_main)
            .unwrap()
            .join()
            .unwrap();
    }
}

二次元の波動方程式を陰解法+ガウス・ザイデル法で解く(Rust) ※バグあり

※バグあり メッシュ数を縦≠横にするとおかしくなる

とりあえずソースと結果を貼っておく。

陰解法と書いたが、実際は陰解法と陽解法を平均した差分方程式を使っている。(クランクニコルソン法っぽくした)

↓の記事で導出した式と、陽解法を組み合わせている。

波動方程式の有限差分法を半自動で導出(sympy) - 雑感等

ソースの条件だとクーラン数=0.1573<1で陽解法でも解ける。

const DT: f64 = 1.0 / 44100.0;

const DX: f64 = 0.05;

このあたりを変えるとクーラン数が変わる。

クーラン数を1以上にしても一応解けたが、上げるとガウス・ザイデル法の反復回数が増える(収束が遅くなる)。

下のソースだと反復回数に上限を設けているから、収束し無くても打ち切る場合がある。

ある程度まで反復回数が増えて、一定を超えると、誤差が増えていくはず。

■出力結果

■Rustプログラムのソース

use std::time::{Duration, Instant};

use plotters::prelude::*;

const T_MAX: usize = 100;
const X_MAX: usize = 50;
const Y_MAX: usize = 50;

const M_ELEM: usize = 3 * 4 + 4 * 2 * (X_MAX - 2) + 4 * 2 * (Y_MAX - 2) + 5 * (X_MAX - 2) * (Y_MAX - 2);


const T_WINDOW: usize = 4;

struct SparseMatrix {
    col: [usize; M_ELEM],
    row: [usize; M_ELEM],
    val: [f64; M_ELEM],
}

fn sub_main() {
    println!("M_ELEM = {:?}", M_ELEM);
    // plot ////////////////////////////////////////////////////////////////////////////////////////
    let root = BitMapBackend::gif("testrustimage.gif", (200, 200), 50).unwrap()
        .into_drawing_area();

    let cells = root.split_evenly((X_MAX, Y_MAX));

    let mut cell_nx;
    let mut cell_ny;
    let mut val;
    let mut col;


    // calc ////////////////////////////////////////////////////////////////////////////////////////
    const DT: f64 = 1.0 / 44100.0;
    const DX: f64 = 0.05;
    const C: f64 = 347.0;
    const COURANT: f64 = C * DT / DX;

    println!("DT:{},DX:{},C:{},COURANT:{}", DT, DX, C, COURANT);

    //u[time][y][x]
    let mut u: [[[f64; X_MAX]; Y_MAX]; T_WINDOW] = [[[0.0f64; X_MAX]; Y_MAX]; T_WINDOW];

    let mut M: [[f64; X_MAX * Y_MAX]; X_MAX * Y_MAX] = [[0.0f64; X_MAX * Y_MAX]; X_MAX * Y_MAX];
    let mut x: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX];
    let mut b: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX];

    let r = 0.5;
    let nu2 = COURANT * COURANT;


    //define M
    {
        let diag = 4.0 * r * nu2 + 2.0;
        let oth_elem = -r * nu2;
        let mut curr_ind = 0;
        for ny in 0..Y_MAX {
            for nx in 0..X_MAX {
                curr_ind = ny * X_MAX + nx;
                
                M[curr_ind][curr_ind] = diag;
                if nx > 0 {
                    M[curr_ind][curr_ind - 1] = oth_elem;
                }
                if nx < X_MAX - 1 {
                    M[curr_ind][curr_ind + 1] = oth_elem;
                }
                if ny > 0 {
                    M[curr_ind][curr_ind - X_MAX] = oth_elem;
                }
                if ny < Y_MAX - 1 {
                    M[curr_ind][curr_ind + X_MAX] = oth_elem;
                }
            }
        }
    }
    let M = M;

    println!("u:{}", u[0][0][X_MAX - 1]);

    let mut u_max = 2.0;
    let mut u_max_t: usize = 0;
    let mut u_max_x: usize = 0;
    let mut u_max_y: usize = 0;
    let mut u_min = -1.0;
    let mut u_min_t: usize = 0;
    let mut u_min_x: usize = 0;
    let mut u_min_y: usize = 0;

    u[1][10][20] = 1.0;


    let mut t_curr;
    let mut t_prev;
    let mut t_next;

    let loop_start = Instant::now();

    //main loop
    //gauss seidel: solve Mx=b
    {
        let omega = 1.0;
        let mut x_temp = 0.0;
        let mut m_temp = 0.0;
        let mut gs_diff = 0.0;

        let mut curr_ind = 0;
        for nt in 1..T_MAX - 1 {
            t_curr = nt % T_WINDOW;
            t_prev = (nt - 1) % T_WINDOW;
            t_next = (nt + 1) % T_WINDOW;

            //define b
            for ny in 1..=Y_MAX - 2 {
                for nx in 1..=X_MAX - 2 {
                    curr_ind = ny * X_MAX + nx;
                    b[curr_ind] =
                        r * nu2 * (u[t_prev][ny][nx - 1] + u[t_prev][ny - 1][nx] + u[t_prev][ny][nx + 1] + u[t_prev][ny + 1][nx])
                            + 2.0 * (1.0 - r) * nu2 * (u[t_curr][ny][nx - 1] + u[t_curr][ny - 1][nx] + u[t_curr][ny][nx + 1] + u[t_curr][ny + 1][nx])
                            - (4.0 * r * nu2 + 2.0) * u[t_prev][ny][nx] + (8.0 * (r - 1.0) * nu2 + 4.0) * u[t_curr][ny][nx];
                }
            }
            //gauss seidel
            for lp in 0_usize..200_usize {
                gs_diff = 0.0;
                for ny in 0..X_MAX * Y_MAX {
                    for nx in 0..ny {
                        m_temp -= M[ny][nx] * x[nx];
                    }
                    for nx in (ny + 1)..X_MAX * Y_MAX {
                        m_temp -= M[ny][nx] * x[nx];
                    }
                    x_temp = (b[ny] + m_temp) / M[ny][ny];
                    gs_diff += (x_temp - x[ny]).abs();
                    // println!("{:?}:{:?}:{:?}", nt, lp, x_temp - x[ny]);
                    x[ny] = omega * x_temp + (1.0 - omega) * x[ny];

                    m_temp = 0.0;
                }
                if gs_diff < f64::MIN_POSITIVE {
                    println!("break {}:{}", nt, lp);
                    break;
                }
            }
            println!("{:?}:{:?}", nt, gs_diff);
            //x->u
            for ny in 0..Y_MAX {
                for nx in 0..X_MAX {
                    curr_ind = ny * X_MAX + nx;
                    u[t_next][ny][nx] = x[curr_ind];
                }
            }
            // plot ////////////////////////////////////////////////////////////////////////////
            for (area, num) in (&cells).into_iter().zip(0..) {
                cell_nx = num % X_MAX;
                cell_ny = num / X_MAX;
                val = u[t_curr][cell_nx][cell_ny];
                col = (val - u_min) / (u_max - u_min);
                area.fill(&RGBColor((col * 255.0) as u8, (col * 255.0) as u8, (col * 255.0) as u8)).unwrap();
            }
            root.present().expect("unable to write");
        }
    }
    let duration = loop_start.elapsed();
    println!("elapsed:{:?}", duration);


    println!("max u[{}][{}][{}]={}", u_max_t, u_max_y, u_max_x, u_max);
    println!("min u[{}][{}][{}]={}", u_min_t, u_min_y, u_min_x, u_min);
}


fn main() {
    {
        const STACK_SIZE: usize = 1024 * 1024 * 1024 * 1024;
        std::thread::Builder::new()
            .stack_size(STACK_SIZE)
            .spawn(sub_main)
            .unwrap()
            .join()
            .unwrap();
    }
}