unwrap

// 1次元の例
// -----------
f = scf();
f.figure_size = [800 1000];
f.figure_position(2) = 0;
f.figure_name = "unwrap() & ""unfold""" + _(": 1-D examples ");
ax = gda();
ax.y_label.font_size=2;
drawlater()
// 元の1次元の外形
t = linspace(-4,4.2,800);
alpha = t.^2 + t -1;
subplot(5,2,1)
titlepage("unwrap(): unwrap | unfold")
subplot(5,2,2)
plot(t,alpha)
t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
ax = gca();
ax.tight_limits = "on";
yT = max(alpha) - strange(alpha)*0.1;
xstring(0,yT,t2)
e = gce();
e.text_box_mode = "centered";
e.font_size = 2;
// 複数の例をループ
for i=1:4
    subplot(5,2,2*i+1)
    if i==1 then
        // atan(tan())でアンラップ
        ralpha = atan(tan(alpha));            // 回復したalpha [pi]
        ylabel("$atan(tan(\alpha))$")
        [u, K] = unwrap(ralpha, %pi);         // arctan
        t2 = "$\text{unwrap(Y, \%pi)}$";
    elseif i==2
        // modulo() + Y-シフトでアンラップ
        c = (rand(1,1)-0.5)*4;
        ralpha = pmodulo(alpha, 5) + c;
        ylabel("$modulo(\alpha,\ 5)"+msprintf("%+5.2f",c)+"$")
        [u, K] = unwrap(ralpha, 0);
        t2 = "$\text{unwrap(Y, 0)}$";
    elseif i==3
        // asin(sin()) + Y-シフトで展開
        ralpha = 1+asin(sin(alpha));          // 回復したalpha [2.pi]
        ylabel("$1+asin(sin(\alpha))$")
        [u, K] = unwrap(ralpha, "unfold");
        t2 = "$\text{unwrap(Y,""unfold"")}$";
    else
        // acos(cos()) + トレンド で折りたたみ
        ralpha = 1+alpha/10+acos(cos(alpha)); // 回復したalpha [2.pi]
        ylabel("$1+\frac{\alpha}{10}+acos(cos(\alpha))$")
        [u, K] = unwrap(ralpha, "unfold");
        t2 = "$\text{unwrap(Y,""unfold"")}$";
    end
    // 処理する外形をプロット
    plot(t, ralpha)
    // 曲線の不連続点または尖点に星印を付けます:
    if K~=[] then
        plot(t(K), ralpha(K),"*")
    end
    // 処理(展開またはアンラップ)した外形をプロットします:
    subplot(5,2,2*i+2)
    plot(t,u)
    ax = gca();
    ax.tight_limits = "on";
    // Adding a legend:
    yT = max(u) - strange(u)*0.2;
    xstring(0,yT,t2)
    e = gce();
    e.text_box_mode = "centered";
    e.font_size = 2;
end
sda();
drawnow()

// 2次元の例
// -------
ax = gda();
ax.title.font_size = 2;
f = scf();
f.color_map = hotcolormap(100);
f.axes_size = [900 450];
f.figure_position(2) = 0;
f.figure_name = "unwrap()" + _(": 2-D examples");
drawlater()

nx = 300;
ny = 400;
rmax = 8.8;
x = linspace(-rmax/2, rmax/2, nx)-1;
y = linspace(-rmax/2, rmax/2, ny)+1;
[X, Y] = meshgrid(x,y);
for ex=0:1   // 例
    // 元の面
        // 面を生成
    if ex==0 then
        z = X.^2 + Y.^2;
    else
        z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
    end

    // 3次元でプロットします
    subplot(2,4,1+4*ex)
    surf(x, y, z)
    title("Original profile Z")
    e = gce();
    e.thickness = 0; // メッシュを削除
    e.parent.tight_limits = "on";

    // ラップされた面 (平面)
    m = 2.8;
    zw = pmodulo(z, m); // ラップする
    subplot(2,4,2+4*ex)
    grayplot(x, y, zw')
    title(msprintf("Zw = pmodulo(Z, %g)  (flat)",m))
    ax0 = gca();
    ax0.tight_limits = "on";

    // アンラップされた面 (平面):
    // in both directions:
    u = unwrap(zw, 0);
    subplot(2,4,3+4*ex)
    grayplot(x, y, u')
    title(msprintf("unwrap(Zw, %g)  (flat)", 0))
    ax3 = gca();
    ax3.tight_limits = "on";

    if ex==0 then
        direc = "r";
    else
        direc = "c";
    end
    // 一方向にアンラップ:
    u = unwrap(zw, m, direc);
    subplot(2,4,4+4*ex)
    grayplot(x, y, u')
    title(msprintf("unwrap(Zw, %g, ""%s"")  (flat)",m,direc))
    ax1 = gca();
    ax1.tight_limits = "on";
end
sda();
drawnow()