Spiral Galaxy Formation Simulation Using MATLAB Function Blocks

Open Model

This example shows how to simulate a spiral arm galaxy formation by using MATLAB Function blocks. This example models the interactions described by Toomre and Toomre, which discusses how disk-shaped galaxies could develop spiral arms [1]. In this example, two disk-shaped galaxies that are far apart fly by each other and almost collide. After the galaxies closely approach each other, mutual gravitational forces cause the galaxies to form spiral arms.

View the Model

Open the SpiralGalaxyExample model to view the design.

The blocks in the shaded area initialize the galaxy structures and orientations, and the MATLAB Function blocks process and plot the data. The Apply Newtonian Gravitation block models the interactions, and the Plot Galaxies block plots the galaxy animation.

Except for the plotting routines in the Plot Galaxies block, all MATLAB Function blocks in this model support code generation with Simulink® Coder™ and Embedded Coder®.

Initialize Galaxies

When you begin the simulation, the model uses Constant blocks in subsystems to specify the initial properties of the two galaxies.

The Radius block specifies the galaxy radius in parsecs.
The Center of Mass block specifies the galaxy mass in solar mass units.
The Position block specifies the galaxy position in parsecs.
The Velocity block specifies the galaxy velocity in meters per second.

The model uses Constant blocks to specify the initial properties of each galaxy. For example, open the Radius block to view the initial radii of the galaxies.

The initial properties cause the galaxies to pass by each other during simulation.

Construct Galaxies

This example creates two galaxies from the initial properties by using the For Each Subsystem block, Construct Galaxies. Construct Galaxies executes the MATLAB Function block, Construct Galaxy, for the initial properties for each galaxy.

To determine which properties to execute, the For Each block uses the inputs as partition inputs. Open the For Each block to view the inputs in the Input Partition tab.

Open the Construct Galaxy block to view the MATLAB® code that constructs the galaxies. In a typical galaxy, most of the mass is concentrated in the center of the galaxy as a super-massive black hole, a cluster of stars, or combination of both. In this example, the block models each galaxy as a disk with a radius of r where most of the mass is concentrated in an inner circle that has a radius of r/3. The block models the input galaxy with a super-massive nucleus, and surrounds the nucleus with 349 stars, for a total of 350 objects. Each star has a mass that ranges from 4 to 24 solar masses. The block randomly positions the stars a distance of r/3 to r away from the center. The stars initially move in circular orbits around the galaxy center. Each object has mass, position, and velocity.

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

Calculate Galaxy Dynamics

The Construct Galaxies subsystem outputs the information about both galaxies as a bus. The MATLAB Function block, Apply Newtonian Gravitation, applies Newtonian mechanics to the bus to establish the physical interactions between the bodies. The code has three subsections:

Partition bodies into heavy and light
Apply Newtonian gravity
Merge heavy and light bodies

The Partition bodies into heavy and light section separates the objects into two groups: heavy bodies and light bodies. The heavy bodies are the galaxy cores and the light bodies are the stars.

%% 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

The Apply Newtonian gravity section uses Newtonian mechanics to compute the velocities and positions of the bodies at each step. Because the galaxy cores are much heavier than individual stars, the model calculates only the heavy-heavy and heavy-light interactions, and ignores the light-light body interactions. Because 698 of the 700 bodies in the model are light, calculating only the heavy interactions significantly reduces the computational complexity.

%% 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

Finally, the Merge heavy and light bodies section merges the data about heavy and light objects together into a single array.

%% 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

Plot Galaxy Interactions

The MATLAB Function block Plot Galaxies plots the consolidated galaxy interaction data in a figure and updates the calculated position of each star at each simulation time step.

The Plot Galaxies block uses the delete function, which is an extrinsic function. Consequentially, the block uses the coder.extrinsic function to declare the delete function before using it. However, coder.extrinsic is not supported for code generation.

Run the Simulation

Run the simulation to view the spiral galaxy formation.

Log Simulation Information

In this example, the model logs the galaxyBodies signal and then saves the logged data to the MATLAB workspace as a Dataset object with the name outputGalaxy. You can retrieve the information on the galaxyBodies signal by retrieving the Simulink.SimulationData.Signal object by using this code:

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

To modify the signal logging settings, right-click the galaxyBodies signal and select Properties.

Create Additional Galaxy Simulations

You can modify this example by adding more galaxies. Add Constant blocks to the Radius, Center of Mass, Position, and Velocity subsystems that correspond to the initial conditions for each galaxy that you want to add.

References

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