メインコンテンツ

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

MATLAB Function ブロックを使用したらせん状の銀河形成シミュレーション

この例では、MATLAB Function ブロックを使用して渦巻腕銀河の形成をシミュレーションする方法について説明します。この例では、Toomre and Toomre によって説明された相互作用をモデル化します。この説明では、円盤状の銀河がどのようにして渦巻腕を発達させることができるかについて議論されています [1]。この例では、遠く離れた 2 つの円盤状の銀河が、衝突しそうなほど互いに接近しています。銀河が互いに接近すると、相互重力がはたらいて銀河の渦巻腕が形成されます。

モデルの表示

SpiralGalaxyExample モデルを開いて設計を表示します。

陰影が付けられた領域内のブロックは銀河の構造と方向を初期化し、MATLAB Function ブロックはデータを処理してプロットします。Apply Newtonian Gravitation ブロックは相互作用をモデル化し、Plot Galaxies ブロックは銀河のアニメーションをプロットします。

Plot Galaxies ブロックのプロット ルーチンを除き、このモデル内のすべての MATLAB Function ブロックで、Simulink® Coder™ と Embedded Coder® を使用したコード生成がサポートされています。

銀河の初期化

シミュレーションを開始すると、モデルはサブシステムの Constant ブロックを使用して 2 つの銀河の初期プロパティを指定します。

  • Radius ブロックは銀河の半径をパーセク単位で指定します。

  • Center of Mass ブロックは、銀河の質量を太陽質量単位で指定します。

  • Position ブロックは銀河の位置をパーセク単位で指定します。

  • Velocity ブロックは、銀河の速度をメートル/秒単位で指定します。

モデルは Constant ブロックを使用して各銀河の初期プロパティを指定します。たとえば、Radius ブロックを開いて銀河の初期半径を表示します。

初期プロパティにより、シミュレーション中に銀河が互いにすれ違います。

銀河の作成

この例では、For Each Subsystem ブロック Construct Galaxies を使用して、初期プロパティから 2 つの銀河を作成します。Construct Galaxies は、各銀河の初期プロパティに対して MATLAB Function ブロック Construct Galaxy を実行します。

実行するプロパティを決定するために、For Each ブロックは入力を分割入力として使用します。For Each ブロックを開いて、[入力分割] タブの入力を表示します。

Construct Galaxy ブロックを開いて、銀河を作成する MATLAB® コードを表示します。典型的な銀河では、質量の大半が超大質量ブラック ホール、星団、またはその両方の組み合わせとして銀河の中心部に集中しています。この例では、ブロックは、半径が r で大半の質量が半径 r/3 の内円に集中している円盤状のものとして各銀河をモデル化します。このブロックは、超大質量の中心核をもつ入力銀河をモデル化し、その中心核を 349 個の星で囲みます。合計で 350 個のオブジェクトが存在します。それぞれの星の質量は、4 ~ 24 太陽質量です。ブロックは、中心から r/3 から r の距離に星をランダムに配置します。これらの星は最初、銀河の中心の周りの円軌道上を進んでいます。各オブジェクトには質量、位置、および速度があります。

function bodies = constructGalaxy(rp,cm,pos,vel)
persistent bs;
numberOfBodies = 350;
if isempty(bs)
    rng('default');
    bs = constructGalaxy0(rp,cm,pos,vel,numberOfBodies);
end
bodies = bs;
function bodies = constructGalaxy0(rp,cm,pos,vel,n)
SolarMass = 1.9891e+30; % In kg
G = 6.672E-11; % Nm^2/kg^2 (Gravitational constant)
SpeedOfLight = 299792458; % in m/s
YearInSeconds = 365*24*60*60;
LightYear = SpeedOfLight*YearInSeconds;
Parsec = 3.26*LightYear;
radiusOuter = rp*Parsec;
radiusInner = (rp/3)*Parsec;
% Each star has X,Y,Z,VX,VY,VZ
% X,Y,Z position in Cartesian coordinates
% VX,VY,VZ velocity in Cartesian coordinates
cm = cm*SolarMass;
bodies = zeros(n,8);
bodies(1,1) = cm;
bodies(1,2) = pos(1)*Parsec;
bodies(1,3) = pos(2)*Parsec;
bodies(1,4) = pos(3)*Parsec;
bodies(1,5) = vel(1);
bodies(1,6) = vel(2);
bodies(1,7) = vel(3);
bodies(1,8) = 'r';
if n > 1
    for i = 2:n
        m0 = rand*20+4;
        m = m0*SolarMass;
        r = rand*(radiusOuter - radiusInner) + radiusInner;
        arg = rand*(2*pi);
        x = r*cos(arg);
        y = r*sin(arg);
        z = 0;
        dx = cos(arg+pi/2);
        dy = sin(arg+pi/2);
        dz = 0;
        % Compute free fall velocity
        v = sqrt(G*cm/r);
        bodies(i,1) = m;
        bodies(i,2) = x+pos(1)*Parsec;
        bodies(i,3) = y+pos(2)*Parsec;
        bodies(i,4) = z+pos(3)*Parsec;
        bodies(i,5) = dx*v+vel(1);
        bodies(i,6) = dy*v+vel(2);
        bodies(i,7) = dz*v+vel(3);
        bodies(i,8) = 'r';
    end
end

銀河のダイナミクスの計算

Construct Galaxies サブシステムは両方の銀河に関する情報をバスとして出力します。MATLAB Function ブロック Apply Newtonian Gravitation は、ニュートン力学をバスに適用して、物体間の物理的相互作用を確立します。コードには次の 3 つのサブセクションがあります。

  • Partition bodies into heavy and light

  • Apply Newtonian gravity

  • Merge heavy and light bodies

Partition bodies into heavy and light セクションでは、オブジェクトが重量体と軽量体という 2 つのグループに分けられます。重量体は銀河の中心核であり、軽量体は星です。

%% Partition bodies into heavy and light
SolarMass = 1.9891e+30; % kg
Limit = 100*SolarMass;
n = size(bodies,1);
props = size(bodies,2);
heavy = zeros(n,props);
light = zeros(n,props);
lightIndex = 1;
heavyIndex = 1;
for i = 1:n
    m = bodies(i,1);
    if m < Limit
        light(lightIndex,:) = bodies(i,:);
        lightIndex = lightIndex + 1;
    else
        heavy(heavyIndex,:) = bodies(i,:);
        heavyIndex = heavyIndex + 1;
    end
end

Apply Newtonian gravity セクションでは、各ステップで物体の速度および位置の計算にニュートン力学が使用されます。銀河の中心核と個々の星とでは質量に大きな差があるため、このモデルでは、重量体と重量体、重量体と軽量体との相互作用のみが計算され、軽量体と軽量体との相互作用は無視されます。このモデルでは 700 個の物体のうち 698 個が軽量体であるため、重量体との相互作用のみを計算することで計算量が大幅に減ります。

%% Apply Newtonian gravity
G = 6.672E-11; % Nm^2/kg^2 (Gravitational constant)
YearInSeconds = 365*24*60*60;
timeStep = 2000000*YearInSeconds;
n = size(heavy,1);
heavy1 = heavy;
light1 = light;
for i = 1:n
    mi = heavy(i,1);
    if mi == 0
        break;
    end
    xi = heavy(i,2);
    yi = heavy(i,3);
    zi = heavy(i,4);
    ar = [0 0 0];
    for j = 1:n
        if i ~= j
            mj = heavy(j,1);
            if mj == 0
                break;
            end
            xj = heavy(j,2);
            yj = heavy(j,3);
            zj = heavy(j,4);
            d = [xj yj zj] - [xi yi zi];
            dr2 = d(1)*d(1)+d(2)*d(2)+d(3)*d(3);
            ar = ar + (d/sqrt(dr2))*((G*mj)/dr2);
        end
    end
    for k = 1:3
        heavy1(i,4+k) = heavy(i,4+k) + ar(k)*timeStep;
    end
end
for i = 1:n
    mi = light(i,1);
    if mi == 0
        break;
    end
    xi = light(i,2);
    yi = light(i,3);
    zi = light(i,4);
    ar = [0 0 0];
    for j = 1:n
        mj = heavy(j,1);
        if mj == 0
            break;
        end
        xj = heavy(j,2);
        yj = heavy(j,3);
        zj = heavy(j,4);
        d = [xj yj zj] - [xi yi zi];
        dr2 = d(1)*d(1)+ d(2)*d(2) + d(3)*d(3);
        ar = ar + (d/sqrt(dr2))*((G*mj)/dr2);
    end
    for k = 1:3
        light1(i,4+k) = light(i,4+k) + ar(k)*timeStep;
    end
end
for i = 1:n
    for k = 1:3
        heavy1(i,k+1) = heavy1(i,k+1) + timeStep*heavy1(i,k+4);
    end
    for k = 1:3
        light1(i,k+1) = light1(i,k+1) + timeStep*light1(i,k+4);
    end
end

最後に、Merge heavy and light bodies セクションで、重量体と軽量体に関するデータが単一の配列にマージされます。

%% Merge heavy and light bodies
n = size(heavy1,1);
nProps = size(heavy1,2);
M = zeros(n,nProps);
for i = 1:n
    if heavy1(i,1) == 0
        break
    end
    M(i,:) = heavy1(i,:);
end
n1 = n - i + 1;
for j = 1:n1
    M(j+i-1,:) = light1(j,:);
    if light1(i,1) == 0
        break
    end
end

銀河の相互作用のプロット

MATLAB Function ブロック Plot Galaxies は、統合された銀河相互作用データを図にプロットし、シミュレーションの各タイム ステップで計算された各星の位置を更新します。

Plot Galaxies ブロックは外部関数である関数 delete を使用します。したがって、ブロックは関数 coder.extrinsic を使用して、関数 delete を使用する前に宣言します。ただし、coder.extrinsic はコード生成ではサポートされていません。

シミュレーションの実行

シミュレーションを実行して、らせん状の銀河形成を表示します。

シミュレーション情報のログ記録

この例では、モデルは galaxyBodies 信号をログに記録し、ログ データを outputGalaxy という名前の Dataset オブジェクトとして MATLAB ワークスペースに保存します。次のコードを使用して Simulink.SimulationData.Signal オブジェクトを取得することで、galaxyBodies 信号に関する情報を取得できます。

simResult = sim("SpiralGalaxyExample.slx");
simResult.outputGalaxy.get('galaxyBodies')

信号ログ設定を変更するには、galaxyBodies 信号を右クリックして、Properties を選択します。

追加の銀河シミュレーションの作成

銀河をさらに追加することで、この例を変更できます。追加する各銀河の初期条件に対応する RadiusCenter of MassPosition、および Velocity の各サブシステムに Constant ブロックを追加します。

参照

[1] Toomre, Alar, and Juri Toomre. “Galactic Bridges and Tails.” The Astrophysical Journal 178 (December 1972): 623. https://doi.org/10.1086/151823.

参考

トピック