import UIKit
import SceneKit
import PlaygroundSupport
// Constants
let c: Double = 299792.458 // speed of light in km/s
let H0: Double = 70 // Hubble constant in km/s/Mpc
let G: Double = 6.6743e-11 // gravitational constant in m^3/kg/s^2
let rho_crit: Double = 1.878e-26 // critical density of the universe in kg/m^3
// Parameters
let N: Int = 1000 // number of particles
let t0: Double = 0 // initial time in Gyr
let tf: Double = 13.8 // final time in Gyr
let dt: Double = 0.01 // time step in Gyr
let a0: Double = 1 / (1 + t0 * H0 / 2997.92) // initial scale factor
// Generate initial conditions
let x = (0..<N).map { _ in
return SCNVector3(
CGFloat(Double.random(in: -10...10)),
CGFloat(Double.random(in: -10...10)),
CGFloat(Double.random(in: -10...10))
)
}
let v = (0..<N).map { _ in
return SCNVector3(
CGFloat(Double.random(in: -100...100)),
CGFloat(Double.random(in: -100...100)),
CGFloat(Double.random(in: -100...100))
)
}
let m = (0..<N).map { _ in
return CGFloat(Double.random(in: 1e11...2e11))
}
// Define Friedmann equations
func HubbleParameter(_ a: Double) -> Double {
return H0 * sqrt((rho_crit / 3) * (1 / pow(a, 3)))
}
func Acceleration(_ x: SCNVector3, _ v: SCNVector3) -> SCNVector3 {
let r = x.length()
let a = SCNVector3(
-G * x.x / pow(r, 3),
-G * x.y / pow(r, 3),
-G * x.z / pow(r, 3)
) * m.reduce(0, +)
a -= SCNVector3(
Float(HubbleParameter(1) * Double(x.x)),
Float(HubbleParameter(1) * Double(x.y)),
Float(HubbleParameter(1) * Double(x.z))
)
return a
}
// Initialize arrays for storing data
let t = stride(from: t0, to: tf, by: dt).map { $0 }
var a = Array(repeating: 0.0, count: t.count)
var v = Array(repeating: SCNVector3Zero, count: t.count)
var x = Array(repeating: SCNVector3Zero, count: t.count)
// Set initial conditions
a[0] = a0
v[0] = v[0]
x[0] = x[0]
// Integrate equations of motion using Verlet algorithm
for i in 1..<t.count {
x[i] = x[i-1] + v[i-1] * Float(dt) + Acceleration(x[i-1], v[i-1]) * pow(Float(dt), 2) / 2
a[i] = a0 * (1 + HubbleParameter(a0) * (t[i] - t0)) // scale factor as a function of time
v[i] = v[i-1]
1
0
358