# Spiral Galaxy Formation Simulation Using MATLAB Function Blocks

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 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`.

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.