"""Spiral Galaxies
Galactic rotation curves without dark matter.
**A new paradigm for the universe.**
https://msp.warwick.ac.uk/~cpr/paradigm
ISBN: 9781973129868
I keep dipping into this book.
Each time with new understanding. Each time with a new perspective.
It is a wonderful work, with compelling arguments.
Chapter 2, Sciama's principle finishes with:
Sciama's initiative, to base a dynamical theory on Mach's principle as
formulated in Sciama's principle, has never been followed up and this
approach to dynamics remains dormant. One of the aims of this book is to
reawaken this approach.
One of the aims of this project is to help understanding of such theories.
In particular, help my own understanding with computer simulations and visualisations of astrophysical data.
This project will explore the mathematics of the book, aiming to make it accessible to a wider audience.
The mathematics of the book is accessible to anyone with a good
understanding of high school calculus. This project aims to extend
that group, to include myself.
The project will require `astropy`, `matpotlib`, `scipy` and of course `numpy`.
Create a Sciama Principle version of the galactic gravitational potential.
Simulations of gravitational wave from a new galactic arrival.
2021 update
===========
The Geometry of the Universe
Colin Rourke
https://www.worldscientific.com/worldscibooks/10.1142/12195
ISBN: 978-981-123-386-9 (hardcover)
A new book, with a bold claim. A stunning story, a revival of some
old ideas and new brilliance too.
Spiral
======
Spiral galaxies? Why spirals?
The problem noted by observational scientist Vera Rubin was that stars
in the outer reaches of spiral galaxies are moving far too fast.
This observation is what lead to the dark matter assumption: there
must be some matter we can't see that is dragging the stars along.
Colin Rourke, says nonesense, there must be a giant rotating mass at
the centre, dragging things along.
More specifically, as we have learnt from the Pulsar Timing Array and
the LIGO/Virgo/Kagra (LVK) detections of gravitational waves is that
space time is itself oscillating.
It's a medium through which waves propagate. Much of the time these
oscillations can be ignored, but sometimes they matter.
It's like boats on the water. A cruise ship can ignore the waves from
a passing dingy, but the wake from a cruise ship is very noticable in
the dingy.
Now General Relativity captures the curvature that is generated by the
high frequency oscillations of matter such as protons and electrons.
Second order effects such as the rotations of bodies can largely be
ignored.
Indeed, under General Relativity, the Kerr metric provides the unique
solution, assuming space is not a vacuum. That metric has the induced
rotation on the surrounding space time proportional to 1/r**3, at
distance r from the rotating body. Provided r is large, it can be
ignored.
The catch is that Kerr just gives the locally induced rotation of
space time. This rotation of space-time will then propagate at the
speed of light, according to the Sciama Principle.
Now in the case of a mature galaxy, a central black hole with mass 10s
or 100s of billions of solar masses is required.
There is growing evidence that galaxies do indeed harbour such
supermassive black holes.
The detection of nanohertz gravitational waves is going to be taken as
evidence of a period in the early universe when galaxies were merging,
hence seeing the gravitational wave background that we are seeing.
My estimate is that the Pulsar Timing Array will see larger and larger
amplitude waves, leading to belief in larger and larger mergers,
earlier and earlier in the universe.
Further, it will be assumed that this is the period in which dark
matter was created.
And the theory is >almost< correct, if we take dark matter as the
wobbles of space time, not some mysterious particle waiting to be
detected.
And if we also recognise that the waves are coming from all the
matter in our visible universe and have been pretty much the same for
a very long time. The Perfect Copernican Principle.
They are not coming from an un-seen phase of the Big Bang theory, but
rather from all the matter that we can see.
Rourke points out an important point about the magnitude of the
gravitational wave background. There is a thing called the
Rees-Sciama effect. Imagine a photon crossing the universe. It
experiences the wobbles of space time.
If the trough of a wave is deeper than the peak, the photon will lose
energy, if it is the reverse, it will gain energy.
So the temperature of the Cosmic Microwave Background should be in
balance with the background waves.
In turn, the magnitude of these waves impacts how far a photon can
travel before being thoroughly diverted from its original course.
This is important, since it in turn has an impact on the temperature
of the Cosmic Microwave Background, the intensity of gamma-ray
bursts and the strength of gravitational waves detected by LVK.
These waves reduce the intensity of gamma-ray bursts since there is a
limit to how long a wave can travel through space before being
thoroughly diverted from its original direction.
Einstein's general relativity withstands double pulsar's Scrutiny
=================================================================
https://physics.aps.org/articles/v14/173
Strong-Field Gravity Tests with the Double Pulsar
M. Kramer et al.
Phys. Rev. X 11, 041050 (2021)
Published December 13, 2021
Binary pulsars provide strong tests of general relativity. The idea
is that there are two, rapidly rotating, neutron stars in a close
orbit around each other.
The Hulse-Taylor binary system has long been the best observed system
and the strongest test.
A new system *PSR J0737-3039A/B* has now taken that crown as the best
observed system. Lets call it Jumb0 737, from the J designate, or
Jumbo for short.
The rotations of the Jumbo are incredibly stable, at the level of an
atomic clock, allowing very precise measurements of the time delays of
the pulses of energy as the system makes each rotation.
In addition, Jumbo is nearer than the Hulse-Taylor system and there is
much less uncertainty in our estimate of its distance from earth.
It all gets complicated, since although the system is in our own
galaxy, the intervening spacetime modulates the signal we receive.
The new paper explains the precise measurements that can now be made
and how variations on general relativity, can, in some cases, be
tested.
In particular, the paper now mentions that rotation of the faster
rotating body now has to be taken into account in the equations of
state of the entire system. The other body is rotating around one
hundred times slower so the effect can largely be ignored.
The aproach that is taken is to assume the Lense-Thiring effect.
This is discussed in some detail in chapter 3 of the book, *the
biggest blunder...*. Particular attention should be paid to the
discussion around page 60.
There are some quite technical arguments. The Lense Thiring effect is
intuitively compelling. General relativity tells us that if
*space-time is a vacuum*, then the Kerr metric is the only solution
that fits Einstein's equations.
This is problematic for the book's arguments as the Kerr metric falls
off as 1/(r**3) whereas the book argues it should fall off as 1/r.
Rourke makes the assumption that linear motion has no inertial effect
and notes that you can change angular momentum by adding a linear
motion, whereas angular velocity cannot be changed in this way.
Is space time a vacuum?
-----------------------
The impasse is resolved if we assume that space is not in fact a
vacuum.
For starters, there is an awful lot of microwave radiation buzzing
around: the cosmic microwave background.
The nature of this background depends critically on our assumptions
about the geometry of the universe: big bang, or static universe?
What would a static universe look like given the effects of special
relativity, simplifying a universe to the set of galaxies above a
certain size?
This model will instead show how by assuming this 1/r relation for the
effect of a masses rotation on its surroundings, produces galactic
rotation curves very much in line with observations.
What about Jumbo?
=================
So the cool news is that we now get estimates of the angular velocity
of each body and it's moment of inertia.
The latter adds some uncertainty to model fit as there is uncertainty
of the exact distribution of the matter within each neutron star,
which is important to know, as the model being used assumed it is
angular momentum that drives the rotational frame dragging.
As noted above, Rourke argues that it is angular velocity, rather than
angular momentum that matters in the calculation, in short the matter
distribution within the black hole is not required for his model.
Further, that the Lense-Thiring effect drops off as 1/r, not the
1/r**3 that is presumably being used by the new paper. The question
is at what point does the Sciama effect become large enough to offset
the Kerr metric, which induces a reduction of the angular momentum of
the system, leading to a belief in black hole inspirals as a source of
gravitational waves.
This may mean that a fit using Rourke's model reduces the
uncertainty in all the parameters that the fitting process estimates,
since the uncertainty in the matter distribution no longer comes into
play, just the uncertainty in angular velocity.
Another project would be to fit the model to the latest Hulse-Taylor
data and see what changes. Calculations show that the objects are
currently too far apart for the Sciama effects to matter. It should
also be noted that the objects have strong magnetic fields.
Calculations also show the potential merger time of some systems is of
the order of a handful of millions, or tens of millions of years. As
this is short on cosmological timescales it raises questions of why
there are so many systems that have not already merged?
I believe the explanation for binary pulsars is that the neutron stars
are the result of supernova explosions and hence the system has
periodically kicked itself apart.
I am also curious how Rourke's model affects the long term evolution
of these systems. My hunch (actually I think I read something along
these lines in the book) is that the feedback from the rotation keeps
these binary systems stable and that it is highly improbable that they
will in-spiral and coalesce.
Which raises the interesting question of what is the source of the
waves that our gravitational wave detectors are seeing?
Which reminds me, I need to work on the grb module.
What about the spiral module?
=============================
The idea is to create a visualisation of the formation of spiral
galaxies with a *omega m / r* model.
It would be good to also be able to model binary systems while we are at it.
2024 July Update
================
This module has now hit 2000 lines and has evolved into a joint
exploration of the Sciama Principle and de Sitter Space.
There are simulations of the gravitational wave background from
distant galaxies and a collection of pieces of code that I have
written to help my own understanding of de Sitter Space.
The full de Sitter Space is a space-time with the Perfect Copernican
Principle: there are no special places or times, it has looked like it
does now for some considerable time.
de Sitter Space has constant negative curvature, and it is this
curvature that means that geodesics, the paths of typical galaxies
through our visible universe follow hyperbolae.
Pairs of geodesics separate exponentially in both forwards and
backwards time. It is this forwards separation that explains redshift.
The separation in backwards time is generally overlooked, but this
balances things out with blue shifted galaxies, new arrivals in our
visible universe.
"""
import traceback
import sys
import argparse
import random
import math
import time
# live dangerously, i almost nver import *, but for math I will make an exception
from math import *
import statistics
from collections import deque, Counter
from functools import partial
from pathlib import Path
from scipy import fft, integrate
import astropy.units as u
import astropy.constants as c
import astropy.coordinates as coord
from astropy import cosmology
from blume import magic
from blume import farm as fm
from blume import taybell
from matplotlib import pyplot as plt
from matplotlib import colors
import numpy as np
from scipy import integrate, signal, optimize
from gotu import sd
[docs]
def mass_length():
""" Equivalence of mass and length based on Schwartzchild radius
This is for use with astropy constants.
for example:
>>> u.solMass.compose(mass_length())
"""
ee = [u.kg, u.m, lambda x: (2 * x * c.G / (c.c * c.c)).value, lambda x:(c.c * c.c * x / (2 * c.G)).value]
def kgtom(x):
return (2 * x * c.G / (c.c * c.c)).value
def mtokg(x):
return (c.c * c.c * x / (2 * c.G)).value
return u.Equivalency([
[u.kg, u.m, kgtom, mtokg]], 'mass_length')
[docs]
def planck_radiance_law_wavelength(wavelength, T=None):
""" Energy emitted in wavelength at temperature T
wavelength: a float, wavelength in meters
T: temperature, a quantity with units K
returns astropy Quantity object with units
"""
if T is None:
T = Cosmo().Tcmb0
hh = c.h
cc = c.c
ee = math.exp
kT = c.k_B * T
lam = wavelength * u.meter
#print(lam)
#print(h*c/lam*kT)
# in terms of frequency v
blam = (2*hh*cc**2/(lam**5)) / (ee(hh*cc/(lam*kT)) - 1)
return blam.value
[docs]
def hubble_tension(cmb=cosmology.Planck18.H0, near=None):
""" Explore hubbe tension """
if near is None:
# I should calculate this based on the full model
near = 73 * u.km / (u.Mpc * u.s)
return near - cmb
[docs]
class Cosmo:
"""Mimic an astropy Cosmology object
It's a Sciama-DeSitter Universe.
Static, but doesn't look static.
Old and big in a state of equilibrium.
The O*0 attributes give the current estimates of the share of the
critical mass (also at z=0) in various buckets.
Ogamma for photons (primarial CMB)
Om for non-realtivistic matter
Ob for baryonic
Odm for dark matter.
Ode for dark energy.
Dark matter is not part of this cosmology.
Instead assume this turns up as black holes in
the centre of galaxies.
In the default astopy cosmology, many functions take a redshift
parameter or z value as a parameter. In FLRW, big bang,
cosmologies, an exact Hubble-law is assumed, so a corresponds
to a time and a distance as well as just redshift.
In a Sciama-DeSitter universe there is only an asymptotic relation
between distance and redshift, for each emitter an observer sees.
There is also an asymptotic relation in backwards time, with
emitters arriving highly blue shifted in our visible universe.
Where z is in effect a time, we return the value for z=0.
"""
def __init__(self, cosmo=None):
cosmo = cosmo or cosmology.default_cosmology.get()
self.cosmo = cosmo
# Our Cosmos is essentially static so all these
# Omega values are the same regardless of z
# Code below is magic to set up these functions
for attr in dir(cosmo):
if attr.endswith('0'):
# don't over-ride anything se set below
if hasattr(self, attr):
continue
def f(z, attr=attr):
return getattr(self, attr)
setattr(self, attr, getattr(cosmo, attr))
setattr(self, attr[:-1], f)
# scale factor for z is useful in de Sitter
# due to curvature things appear nearer and
# smaller than they actually are.
# In FLRW, scale_factor gives the relative size at a given z.
# as 1/(z+1)
# In a de Sitter model this is the factor you have to apply
# to convert an actual distance to an apparent distance.
# due to curvature, everything we see appears to be
# within the hubble_distance.
self.scale_factor = cosmo.scale_factor
# turn dark matter and dark energy into matter
# I now believe dark matter is better captured as "the energy in the oscillations of
# space time itself".
# dark energy is not really a thing - what we can show is what the value of dakr
# energy would be if you mistakenly assume a de Sitter Universe is a Big Bang universe.
self.Ob0 += self.Ode0 + self.Odm0
self.Obh0 = self.Ode0 + self.Odm0 # mass of central black holes
self.Ode0 = 0.
self.Odm0 = 0.
@property
def hubble_distance(self):
htime = (1/self.H0).to(u.year)
distance = htime.value * u.lightyear
return distance
[docs]
def check_critical_density(self):
""" Critical density formulae
The critical density only depends on the Hubble distance.
It is the mass per unit volume such that the total mass has
Schwartzchild radius equal to the Hubble distance.
As such, it is determined by the Hubble constant.
"""
pi = math.pi
hd = self.hubble_distance << u.m
cd = self.critical_density0 << u.kg / u.m**3
# mass of the universe [1]
M = cd * (4/3) * pi * hd * hd * hd
print(hd)
print(cd)
print(M)
# Schwartzchild radius
print((c.G * 2 * M / (c.c * c.c)) << u.m)
# hd = 2 * M * c.G / (c.c * c.c)
# M = hd * c.c * c.c / (2 * c.G)
sc = (c.G * 2 * (cd * (4/3) * pi * hd * hd * hd) / (c.c * c.c))
# substitute from [1]
cd2 = (3/4) * c.c * c.c / (2 * c.G * hd * hd * pi)
print(cd, cd2, M)
[docs]
def is_flat(self):
""" de Sitter Space has constant negative curvature
Hence, it is not flat.
"""
return False
[docs]
class SkyMap(magic.Ball):
"""Yet another table viewer?
parsing csv and figuring out what it means.
or just give me lists of fields i can pick form for
any purpose?
Or yet another universe viewer.
And so we descend into the world of coordinate systems.
Something astronomers understand very well.
See :mod:`gotu.wits` for more on this subject.
Maybe SkyMap allows systems to evolve over time, according to the
paradigm.
"""
def __init__(self, gals=None, n=1566, fudge=None):
"""
gals: list of Spiral() objects to initialise with
n: if gals is None, then generate a sample of size n
fudge: value for self.fudge
"""
super().__init__()
self.offset = 0.
# set defaults from Cosmo object
cosmo = Cosmo()
self.cosmo = cosmo
# Conversion factor from stellar mass to black hole mass
# Assume the dark matter and dark energy turn up as central
# galactic black holes.
# more conservative, just use the dark matter.
# self.m_bh = (cosmo.Ode0 + cosmo.Odm0) / cosmo.Ob0
self.m_bh = cosmo.Obh0
# the CMB is 45 times brighter than you would expect based
# on size of visible universe. Hitchhikers should note
# that setting fudge to 42 produces a viable model.
self.fudge = fudge or 45.
self.n = n or 1566
# Conversion from Holmberg radius mass to full stellar mass
self.minmass = 6. # min mass to include, log10(solar_mass)
hd = cosmo.hubble_distance << u.m
volume_of_universe = (4/3) * math.pi * (hd**3)
self.volume_of_universe = volume_of_universe
critical_density = cosmo.critical_density0 << u.kg / u.m**3
mass_of_universe = critical_density * volume_of_universe
self.mass_of_universe = mass_of_universe
self.period = 30 # years
self.fftlen = 32
self.cmb_min = 1e-4 # u.m
self.cmb_max = 0.004 # u.m
self.delta_t = 0.001 # time inc in natural units
self.t = 0
# local simulation only
self.max_distance = 1e8 * u.lightyear
self.filter = True
self.minz = 1.
self.min_phi = 0.001
self.max_phi = 2.0
self.max_t0 = self.delta_t * 1000
self.maxtheta = None
self.mintheta = 0
self.tborigin = True
self.toff = 0.00001 # small offset to avoid overflow when calculating
# uoft for t = tstar
self.scale_for_curvature = False
self.create_sample(gals)
self.good = dict()
self.tablecounts = magic.TableCounts(
xname='z', yname='distance', maxx=10, maxy=10)
self.tasks = magic.Tasks()
self.tasks.add(self.local_mode_sim)
self.tasks.add(self.cmbsim)
self.tasks.add(self.cmb_gwb)
self.tasks.add(self.gravity_waves)
self.tasks.add(self.uoft)
#await self.local_mode_sim()
#await self.cmbsim()
#await self.cmb_gwb()
[docs]
def create_sample(self, gals=None, n=None, t0=0):
"""Create a new sample of galaxies
This also initialises some parameters, based on the sample.
The goal is to set the stellar masses to match the baryonic
matter as given by the cosmology Ob0. This is then scaled so
that the sum of the central masses equals the mass of the
universe.
gals: if given, a list of galaxies to include.
"""
self.balls = gals or []
dt = self.delta_t
#dt = 0
self.balls += list(sample_galaxies(
n or self.n, fudge=self.max_phi**2, t0=t0, delta_t=dt,
min_phi=self.min_phi,
cosmo=self.cosmo))
# recalculate mean sample_mass
mean_sample_mass = self.sample_mass()
stellar_mass = self.cosmo.Ob0 * self.mass_of_universe << u.solMass
# Given number of galaxies, we can get the factor needed to scale
# the sample mass up to the mass of baryonic matter.
self.N = 1e11
self.h_factor = (stellar_mass / self.N) / mean_sample_mass
# initialise Mcent -- done this way because we need
# to calculate h_factor based on mean sample mass
self.set_mcent()
[docs]
def set_mcent(self):
""" Set up central mass of all spheres based on stellar mass """
msun = schwartzchild(c.M_sun)
for gal in self.balls:
stellar = (msun * gal.Mstellar) << u.lightyear
black_hole = stellar.value * self.h_factor * self.m_bh
gal.Mcent = black_hole.value
[docs]
def sample_mass(self):
""" Return the sum of the masses of the sample
In solar masses.
"""
mass = 0
total = 0
for gal in self.balls:
m26 = gal.Mstellar
if math.log10(m26) < self.minmass:
continue
total += 1
mass += m26
return mass / total << u.solMass
[docs]
async def local_mode_sim(self):
""" Simulation of the gravity waves from the local cluster """
wavelengths = []
amps = []
for gal in self.balls:
if not gal.Mstellar: continue
distance = gal.distance
if distance > self.max_distance: continue
stellar_mass = self.h_factor * gal.Mstellar
black_hole = gal.lightyear_to_kg()
sc = (2 * black_hole * c.G / (c.c*c.c)) << u.lightyear
wavelengths.append(sc.value)
amps.append((sc / distance).value)
points, wave = waves(
amps, wavelengths,
period=self.period)
ax = await self.get()
ax.plot(points, wave)
ax.set_title('Local Gravitational Wave Background')
ax.show()
[docs]
def waves(self):
""" What waves would be generated by the catalog? """
sample = []
total_mass = 0.
for gal in self.balls:
stellar = gal.Mstellar
name = gal.name
if math.log10(stellar) < self.minmass:
continue
sample.append(gal)
bh = gal.lightyear_to_kg() / c.M_sun
total_mass += bh.value + stellar
print('total mass', math.log10(total_mass), total_mass * c.M_sun)
mean_mass = math.log10(total_mass / len(sample))
print('mean mass', mean_mass)
print('sample size', len(sample))
print(f'total mass given {self.N} galaxies',
self.N * 10**mean_mass * c.M_sun)
mou = self.mass_of_universe
print("number of galaxies based on mbocd",
f'{mou / (c.M_sun * total_mass/len(sample)):.2e}')
# weight for each galaxy
weight = self.N / len(sample)
R = self.cosmo.hubble_distance << u.m
for gal in sample:
sc = gal.schwartzchild() << u.m
gal.amplitude = 1.5 * (weight * sc)/R
return sample
def decra2rad(self, dec, ra):
return math.radians(dec), math.radians(ra)
ra = (ra - 12) * math.pi / 12.
while ra > math.pi:
ra -= 2 * math.pi
return dec * math.pi / 180., ra
def spinra(self, ra):
ra += self.offset
while ra > math.pi:
ra -= 2 * math.pi
while ra < math.pi * -1:
ra += 2 * math.pi
return ra
[docs]
def tofu(self, u, balls=None, name=None):
""" Work in natural units, c=G=Hubble Distance=1
Time in Hubble times.
"""
balls = balls or self.balls
a = np.cosh([ball.phi for ball in balls])
d = np.cos([ball.theta for ball in balls])
# work with U = e**u and T=e**t
A = D = (a - d)/2
B = C = (a + d)/2
# first time visible is tstar given by
name = name or 'tofu'
tstar = np.log(np.sqrt(A/B))
ax = self.get_nowait()
ax.plot(tstar, 'o', label=f'{name[0]}*')
# the last time u that the source can be seen
# notice e^umax = 1/(e**tstar = e**(-tstar)
# so umax = -1 * tstar.
# t for umax is infinite and u for tstar is -infinity
umax = np.log(np.sqrt(B/A))
ax.plot(umax, 'o', label=f'{name[3]}max')
ax.legend()
ax.show()
U = e**u
T = ((U + np.sqrt(A*B + ((1 - B*B - A*A) * U*U) + A * B* U*U*U*U))
/ (B - A * U*U))
return np.log(T)
def uoft(self, t=0, balls=None):
# use distance as t0, first time visible in our galaxy.
hd = self.cosmo.hubble_distance
balls = balls or self.balls
t0 = [ball.t0 for ball in balls]
tb = [ball.tb() for ball in balls]
phi = [ball.phi for ball in balls]
theta = [ball.theta for ball in balls]
a = np.cosh(phi)
d = np.cos(theta)
A = D = (a - d)/2
B = C = (a + d)/2
# first time visible is tstar given by
tstar = np.log(np.sqrt(A/B))
# the last time u that the source can be seen
# notice e^umax = 1/(e**tstar = e**(-tstar)
# so umax = -1 * tstar.
# t for umax is infinite and u for tstar is -infinity
umax = np.log(np.sqrt(B/A))
# we want u for tt = tstar + t - t0
tt = tb + tstar + t - t0
uu = [ball.uoft(t) for t, ball in zip(tt, balls)]
zz = []
xx = []
for uuu, ttt, ball in zip(uu, tt, balls):
z, x = ball.zandx(ttt, uuu)
zz.append(z)
xx.append(x)
tb = [ball.tb() for ball in balls]
results = dict(distance=xx, z=zz, u=uu, t=tt, t0=t0, tb=tb,
phi=phi, theta=theta)
return results
[docs]
def find_distant(self, results):
""" Find objects we no longer care about """
max_distance = sqrt(self.fudge)
removals = set()
xx = results['distance']
zz = results['z']
dd = set()
for distance, ball in zip(xx, self.balls):
if distance > max_distance:
removals.add(ball)
print('Removals', len(removals), len(dd), len(dd & removals), max_distance)
return removals
def remove(self, balls):
for ball in balls:
self.balls.remove(ball)
async def ticker(self):
delta_t = self.delta_t / 100.
self.t = 0.
self.balls = []
while True:
self.t += delta_t
done = await self.tick()
if done: return
async def counter(self):
self.create_sample()
maxtheta = self.maxtheta
mintheta = self.mintheta
tt = 0.
xname, yname = 'z', 'r'
for ball in self.balls:
t1 = time.time()
if maxtheta or mintheta:
maxcos, mincos = cos(maxtheta), cos(mintheta)
# need uniform numbers in range [maxcos, mincos]
if maxcos < mincos:
maxcos, mincos = mincos, maxcos
rand = (random.random() * (maxcos-mincos)) + mincos
ball.theta = math.acos(rand)
t = ball.tstar() + self.toff
if self.tborigin:
t += ball.tb()
while True:
z, x = ball.zandx(t)
# adjust distance for scale factor, adjust z too...
# ie repace x by x/1+x, ditto z.
if self.scale_for_curvature:
xname, yname = 'z/(1+z)', 'r/(1+r)'
x *= self.cosmo.scale_factor(x)
z *= self.cosmo.scale_factor(z)
if z < 0:
z = z/(1+z)
x = x/(1+x)
if z > self.tablecounts.maxx:
break
if z > 0 and x > self.tablecounts.maxy:
break
#weight = 1/(x*x)
weight = 1
self.tablecounts.update([z], [x], weight)
t += self.delta_t
#print(f'{t:8.2f} {ball.theta:6.3f} {ball.phi:6.2f} {z:6.2f} {x:6.2f}')
t2 = time.time()
tt += t2-t1
if tt > self.sleep:
await self.tablecounts.show(xname=xname, yname=yname)
tt = 0.
# one last show
await self.tablecounts.show(xname=xname, yname=yname)
async def walker(self, inc=0.1, t0=0., maxtheta=None):
self.create_sample()
for ball in self.balls:
ball.t0 = 0.
if maxtheta:
maxacos = cos(maxtheta)
xxx = 1./(1-maxacos)
ball.theta = math.acos(1.0-(random.random()/xxx))
self.good = {}
self.t = t0 or 0.
while True:
self.t += inc
done = await self.tick()
if done: return
await magic.sleep(self.sleep)
async def tick(self):
#self.create_sample()
#t = self.delta_t
t = self.t
result = self.uoft(t)
# filter sample to get those nearer than the sqrt(self.fudge)
max_distance = sqrt(self.fudge)
# get array of bools indicating values we want to keep
if self.filter:
keep = [(redshift < max_distance) or
(distance <= max_distance and t0 < self.max_t0)
for redshift, distance, t0
in zip(result['z'], result['distance'], result['t0'])]
if np.any(keep):
if self.good:
print(f'Keeping {len(keep)} t={t} good={len(self.good["theta"])}')
else:
print(f'Keeping {len(keep)} t={t}')
for key, value in result.items():
keepers = np.array(value)[keep]
if key in self.good:
self.good[key] = np.concatenate((self.good[key], keepers))
else:
self.good[key] = keepers
else:
return True
# delete any balls not keeping
for keeper, ball in zip(keep, self.balls):
if not keeper:
self.balls.remove(ball)
result = self.good
await self.showzdimage()
return
zz = np.array(result['z'])
t0 = np.array(result['t0'])
distance = np.array(result['distance'])
# do some plotting
ax = await self.get()
colours = magic.random_colour()
ax.scatter(zz, distance, c=result['theta'],
cmap=colours)
zline = np.linspace(min(distance), max(distance), 100)
line = statistics.linear_regression(zz, distance, proportional=True)
ax.plot(zline, zline)
ax.plot(zline, (zline * line.slope) + line.intercept)
ax.grid(True)
ax.set_title(
f't={t:.3f} z v distance, t0={min(t0):.4} to {max(t0):.4}' +
f' max(phi)={self.max_phi:.2} slope:{line.slope:.2}')
ax.show()
mask = (zz>0) & (distance < 2)
redszz = zz[mask]
redsdist = distance[mask]
if len(redszz) > 2:
axz = await self.get()
axz.scatter(redszz, redsdist, c=result['theta'][mask],
cmap=colours)
zline = np.linspace(min(redsdist), max(redsdist), 100)
line = statistics.linear_regression(redszz, redsdist, proportional=True)
axz.plot(zline, zline)
axz.plot(zline, (zline * line.slope) + line.intercept)
axz.set_title(f"Best fit slope: {line.slope:.4}")
axz.grid(True)
axz.show()
ax = await self.get()
ax.set_title('Theta')
ax.hist(np.array(result['theta']), 20)
ax.show()
ax = await self.get()
ax.set_title('Phi')
ax.hist(np.array(result['phi']), 20)
ax.show()
await self.showzdimage()
return
async def showzdimage(self):
result = self.good
zz = np.array(result['z'])
t0 = np.array(result['t0'])
distance = np.array(result['distance'])
size = 100, 100
self.grid = np.zeros(size)
grid = self.grid
minz, maxz = 0, sqrt(self.fudge)
mind, maxd = minz, maxz
zd = distance.clip(mind, maxd)
zz = zz.clip(minz, maxz)
dsize = (maxd-mind)/size[0]
zsize = (maxz-minz)/size[1]
for d, z in zip(zd, zz):
#print(z,d)
zbucket = int((z-minz) // zsize)
dbucket = int((d-mind) // dsize)
grid[dbucket, zbucket] += 1
grid = grid[1:-1, 1:-1]
norms = grid / (sum(grid)+1)
#print(sum(grid))
ax = await self.get()
extent = (minz, maxz, mind, maxd)
ax.imshow(norms,
origin='lower',
aspect='auto',
extent=extent,
cmap=magic.random_colour())
ax.show()
# see what a grid sample looks like
width, height = grid.shape
csize = width * height
choices = list(range(csize))
sample = random.choices(choices, weights=grid.flatten(), k=1566)
ax = await self.get()
xx = [x%width for x in sample]
yy = [int(x/width) for x in sample]
ax.scatter(xx, yy)
#ax.hist(sample)
ax.show()
ax = await self.get()
ax.imshow(grid[1:-1, 1:-1], origin='lower', aspect='auto')
ax.show()
async def explore(self, thetas=None, phis=None, thetaphis=None, samples=None):
ax1 = await self.get()
ax2 = await self.get()
ax3 = await self.get()
ax4 = await self.get()
# ignore thetas, phis
if samples:
balls = random.sample(self.balls, samples)
thetaphis = sorted([ (ball.theta, ball.phi) for ball in balls])
if thetas is not None:
thetaphis = [(theta, phi) for theta in thetas for phi in phis]
for theta, phi in thetaphis:
await self.plotthetaphi(theta, phi, self.max_t0,
ax1, ax2, ax3, ax4)
await magic.sleep(0)
async def plotthetaphi(self, theta, phi, tmax,
ax1=None, ax2=None, ax3=None, ax4=None):
a = cosh(phi)
d = cos(theta)
A = D = (a - d)/2
B = C = (a + d)/2
epsilon = 1e-20
# first time visible is tstar given by
tstar = np.log(np.sqrt(A/B))
tb = log((sqrt(a+1) + sqrt(1-d))/sqrt(a-d))
if self.tborigin:
t = np.linspace(tb+tstar, tb+tstar+tmax, 5000)
ttb = t-(tstar+tb)
title = 't-(t*+tb) v '
else:
t = np.linspace(epsilon+tstar, tstar+tmax, 5000)
ttb = t - tstar
title = 't-t* v '
try:
uu = [uoft(tt, theta, phi) for tt in t]
zx = [zandx(tt, u, theta, phi) for u, tt in zip(uu,t)]
zz = np.array([zzx[0] for zzx in zx])
xx = np.array([zzx[1] for zzx in zx])
except Exception as e:
traceback.print_exception(e)
print('in uoft', t[0], theta, phi)
return
mask = xx < self.max_distance.value
ax = ax1 or await self.get()
thetaphi = f'phi={phi:.3} theta={theta/pi:.2} tb={tb:.3}'
ax.set_title('t-t* v u')
ax.plot(t-tstar, uu, label=thetaphi)
ax.legend()
ax.show()
ax = ax2 or await self.get()
ax.set_title(title + 'z')
ax.plot(ttb, zz, label=thetaphi)
ax.legend()
ax.show()
ax = ax3 or await self.get()
ax.set_title(title + 'd')
ax.plot(ttb, xx, label=thetaphi)
ax.legend()
ax.show()
ax = ax4 or await self.get()
ax.set_title(f'z v d')
ax.plot(zz[mask], xx[mask], 'o', label=thetaphi)
ax.mlegend()
ax.show()
async def log10hist(self, values, n=10, title=None):
ax = await self.get()
ax.hist([math.log10(x) for x in values])
ticks = ax.get_xticks()
ax.set_xticklabels([10**x for x in ticks])
ax.set_title(title)
ax.show()
async def run(self):
self.set_mcent()
#ax = fig.add_axes((0,0,1,1), projection='mollweide')
ax = await self.get()
ax.projection('mollweide')
locs = [self.decra2rad(ball.dec, ball.ra)
for ball in self.balls]
ball_colours = [x.distance.value for x in self.balls]
ax.set_facecolor('black')
ax.scatter([self.spinra(xx[1]) for xx in locs],
[xx[0] for xx in locs],
c=ball_colours,
s=[math.log10(x.Mstellar) or 1
for x in self.balls])
ax.axis('off')
ax.show()
ax = await self.get()
sample = self.waves()
ax.set_title('Schwartzchild v amplitude')
ax.scatter([x.schwartzchild().value for x in sample],
[x.amplitude.value for x in sample])
await self.log10hist(
[x.amplitude.value * self.fudge for x in sample],
title='Amplitude(m)')
await self.log10hist(
[x.schwartzchild().value for x in sample],
title='Schwartzschild (lightyear)')
ax.show()
ax = await self.get()
ax.hist([math.log10(x.lightyear_to_kg() / c.M_sun) for x in sample])
ax.set_title('Black hole mass (log10(suns))')
ax.show()
ax = await self.get()
ax.scatter([x.Mcent for x in sample],
[x.amplitude.value for x in sample])
ax.set_title('Black hole mass v amplitude')
ax.show()
await super().run()
async def gravity_waves(self):
# try and simulate the wave
period = self.period
sample = self.waves()
ax = await self.get()
points, wave = waves(
[x.amplitude.value * self.fudge for x in sample],
[x.schwartzchild().value for x in sample],
period=self.period)
ax.plot(points, wave)
ax.set_title('Gravitational Wave Background')
ax.show()
f, t, Sxx = signal.spectrogram(wave, nperseg=self.fftlen)
print(len(wave), len(f), len(t), Sxx.shape)
ax = await self.get()
ax.plot(t)
ax.set_title('t for spect')
ax.show()
ax = await self.get()
ax.plot(f)
ax.set_title('f for spect')
ax.show()
ax = await self.get()
ax.pcolormesh(t, f, Sxx, shading='gouraud')
ax.ylabel('Frequency [Hz]')
ax.xlabel('Time [sec]')
ax.show()
ax = await self.get()
ax.imshow(Sxx, cmap=magic.random_colour())
ax.set_title('imshow of spectrogram')
ax.show()
[docs]
async def cmbsim(self):
""" Calculate some CMB related numbers
This is really just a test of
planck_radiance_law_wavelength, checking that the
energy gives the same Ogamm0 as is implied by
the temperature.
"""
# Mass per unit volume, based on critical density and Ogamma
# Assumes CMB makes up the vast majority of photons.
cosmo = self.cosmo
cmb_mpuv = cosmo.Ogamma0 * cosmo.critical_density0 << u.kg/u.m**3
print("CMB equivalent mass per unit volume", cmb_mpuv)
# What is the energy
mwave = 1e-3 * u.m
mass_mwave = (c.h * c.c / mwave).to(u.kg, equivalencies=u.mass_energy())
print(f"energy of {mwave} photon {mass_mwave}")
print("photon density", cmb_mpuv / mass_mwave)
# energy according to planck radiance
energy, error = integrate.quad(
partial(planck_radiance_law_wavelength, T=cosmo.Tcmb0),
self.cmb_min, self.cmb_max)
# this should relate to cmb per unit volume
# It is the energy per unit solid angle,
# per square meter per second.
# so we want to scale by: 4 * math.pi / c.c
cmb_energy_density = (4 * math.pi * energy / c.c).value
print("Energy per cubic meter", cmb_energy_density)
# mass equivalent
mass_equiv = (cmb_energy_density / (c.c * c.c)).value
print("Mass equivalent per cubic meter:", mass_equiv)
cmbgamma = (mass_equiv /
self.cosmo.critical_density0.to(u.kg/u.m**3)).value
# compare to critical density
print("CMB fraction of critical, Ogamma", cmbgamma)
print("Ogamm0", self.cosmo.Ogamma0)
async def cmb_gwb(self):
cosmo = self.cosmo
# Simulate the gravitational waves that the CMB generates
a, b = self.cmb_min * u.m, self.cmb_max * u.m
wavelengths = np.linspace(a, b , 1000)
# Planck radiance energy, converted to J/m3
density = np.array(
[planck_radiance_law_wavelength(x.value) / c.c.value
for x in wavelengths])
# convert to total CMB energy in universe
energy = density * (self.volume_of_universe << u.m**3).value
# Get Schwartzchild radii
sc = np.array([
(2 * c.G * c.h / (x * c.c**3)).value for x in wavelengths])
# radius of visible universe in m
R = cosmo.hubble_distance.to(u.m).value
# total mass of CMB
M = cosmo.Ogamma0 * self.mass_of_universe
# Mass of single photon of each wavelength
mass_per_photon = (c.h * c.c / wavelengths).to(
u.kg, equivalencies=u.mass_energy())
total_photons = energy / mass_per_photon
# 4/3 pi R3 * energy -- rho 4 pi r2
# The effect integrated over volume of the universe 3M / 2R
# the wave decays with 1/r, and we have to divide energy
# by the volume
energy = (sc * 3 * energy / (2 * R))
# can also calculate energy from density
# energy2 = sc * 2 * math.pi * density * R * R
ax = await self.get()
ax.set_title('Mass per photon')
ax.plot(wavelengths, mass_per_photon)
ax.show()
ax = await self.get()
ax.set_title('Total photons')
ax.plot(wavelengths, total_photons)
ax.show()
ax = await self.get()
ax.set_title('planck radiance energy')
ax.plot(wavelengths, energy)
ax.show()
ax = await self.get()
ax.set_title("Schwartzschild Radius")
ax.plot(wavelengths, sc)
ax.show()
ax = await self.get()
ax.plot(wavelengths, energy * self.fudge)
ax.set_title('Scaled energy')
ax.show()
# simulate the waves
points, wave = waves(
[x for x in energy],
[x.value for x in wavelengths],
period=self.period)
ax = await self.get()
ax.set_title('CMB generated gravity waves')
ax.plot(points, wave)
ax.show()
ax = await self.get()
ax.set_title('frequency v energy')
ax.plot([c.c.value/x for x in points], energy)
ax.show()
await self.htension(wavelengths, energy)
async def htension(self, wavelengths, energy):
# Now take the ratio of this amplitude to the wavelength
# and see how it compares to the Hubble constant in Hz.
# Why should there be a relation? We are multiplying
# by the Hubble time, to get the size of the effect
# over the Hubble time - it appears commensurate
# with the Hubble tension.
h0 = self.cosmo.H0.to(1/u.s).value
hubble_tension = np.array([
x.value/h0 for x in self.fudge * energy / wavelengths])
# Plot some stuff
ax = await self.get()
ax.set_title("Hubble Tension")
ax.plot(wavelengths, (1+hubble_tension) * self.cosmo.H0)
ax.plot(wavelengths, [self.cosmo.H0.value] * len(wavelengths))
ax.show()
def waves(amplitudes, wavelengths, phases=None, period=1, n=1000):
wave = np.zeros(n)
if phases is None:
phases = 2 * math.pi * np.random.random(n)
points = np.linspace(0, period, n)
for amplitude, wavelength, phase in zip(
amplitudes, wavelengths, phases):
wave += amplitude * np.sin((points/wavelength) + phase)
return points, wave
[docs]
class Spiral(magic.Ball):
""" Model a spiral galaxy
Or any rotating mass?
Want to convert this to use astropy units.
"""
def __init__(self):
super().__init__()
# set up an initial mode
self.cosmo = Cosmo()
# need a way to calculate the density of the medium
# and the temperature -- the two are linked.
# this is for z and Eddington sphere calculations.
self.n = 1e5 # density in protons / m**3 1e2 - 1e12
self.T = 3000
self.galaxy()
self.details = True
[docs]
def galaxy(self):
""" Set parameters for a galaxy """
# A = K * \omega_0. K = M for Sciama principle
# note A is half the asymptotic tangential velocity.
self.A = (110 * u.km/u.s / c.c).decompose().value
# Apparent rate of precession of the roots of the spiral.
self.B = 0.00000015
# Central mass. Mass converted to Schwartzschild radius (in light years)
# Mass of 1 is approximately 3e12 solar masses.
self.set_mcent(0.03)
self.Mball = 0.
self.Mdisc = 0.
# magic constant determined by overall energy in the orbit
self.EE = -0.00000345
# constant, can be read from tangential velocity for small r
self.CC = -10
# range of radius r to consider, in light years
self.rmin = 5000
self.rmax = 100000
def set_mcent(self, mcent):
self.Mcent = mcent
self.K = self.Mcent
self.omega0 = self.A / self.K # angular velocity in radians per year
def print_parms(self):
print('omega0', self.omega0)
print('CC', self.CC)
print('A/K', self.A / self.K)
print('rmin_check', self.rmin_check())
[docs]
def rmin_check(self):
""" The length of the roots of the spirals
This can be used to set the B value.
Assume that the spiral roots end at radius r0
And assume the roots are moving with the inertial frame at that
radius.
The rate of precession will match that of the inertial frame at
that radius.
"""
return self.A / self.B
def find_cc(self, tangential_velocity, r=None):
# constant, can be read from tangential velocity for small r
r = r or self.rmin
self.CC = 0
vtan = self.v(r)
ccoverr = tangential_velocity - vtan
self.CC = ccoverr * r
return self.CC
[docs]
def v(self, r):
""" Velocity at radius r
A = 0.0005
K = Mcent
CC = -10
??
"""
A = self.A
K = self.K
CC = self.CC
return (2 * A) - (2 * K * A * math.log(1 + K) / r) + CC / r
[docs]
def vinert(self, r, v):
""" Inertial part of the velocity
Part of velocity relative to inertial frame.
Notes
-----
K is central mass. A = 0.0005
"""
return v - (self.A * r) / (self.K + r)
def rdoubledot(self, r, vinert):
rdd = ((vinert ** 2) / r) - (self.Mcent/(r**2))
# if we have Mdisc of Mball, adjust as appropriate?
rdd -= self.Mdisc/(self.rmax ** 2)
rdd -= self.Mball * r /(self.rmax ** 3)
return rdd
def energy(self, r):
CC = self.CC
Mcent = self.Mcent
Mdisc = self.Mdisc
Mball = self.Mball
rmax = self.rmax
EE = self.EE
K = self.K
A = self.A
Log = math.log
# ok this deserves an explanation!
energy = (-CC**2/(2*r**2) + (Mcent - 2*A*CC)/r -
# adjustment for a uniform disk radius rmax, mass Mdisc
Mdisc*r/rmax**2 +
# adjustment for a spherical mass
Mball*r**2/(2*rmax**3) +
#
A**2*K/(K + r) +
A**2*Log(K + r) +
2 * A*K * (CC + 2*A*r) * Log(1 + r/K)/(r**2)
- (2 * A*K*Log(1 + r/K)/r)**2 + EE);
return energy
[docs]
def spheres3(self):
""" Bondi, Eddington and Schwartzchild
We have mass in lightyears, the Schwartzchild radius.
To get the mass in kilograms, S = 2GM/c^2, M = Sc^2/2G
"""
T = self.T or 3000.
S = self.Mcent # mass in lightyears, we need mass in kg
M = self.lightyear_to_kg()
msuns = M / c.M_sun
print("mass in suns:", msuns)
print("Bondi:", self.bondi())
print("Acretion:", self.accretion() << u.solMass/u.year)
print('Eddington:', self.eddington())
print('Schwartzchild:', self.schwartzchild())
self.schwartzchild()
[docs]
def density(self):
""" Return density of Mcent if it is a black hole """
S = self.schwartzchild() << u.m
M = S * c.c * c.c/ (2*c.G)
M = M.to(u.kg, equivalencies=u.mass_energy())
mpcm = (M / ((4/3) * math.pi * (S**3))) << u.kg / u.m**3
return mpcm
[docs]
def amplitude(self, R=None):
"""amplitude of the gravitational wave wave
Returns the amplitude of the wave if there was
just one such object in the sphere radius R. Scale
by the actual density to get useful numbers.
Suppose rho is the number of such objects per unit volume,
the the full wave will be::
a = 4 * pi rho Sum M r^2 dr / r for 0 < r < R
Which simplifies to::
a = 4 * pi * rho * M Sum r dr for 0 < r < R
or:
a = 2 * pi * rho * S * R ** 2 [1]
Now::
N = rho * 4/3 * pi * R**3
So if N==1, rho = 4 * pi / (3* R**3), substituting in [1]
a = 3 S / 2 R
Where S be the Schwartzchild radius of Mcent.
R is the radius out to which bodies can contribute to the sum.
To get the full wave, scale by the number of objects in
the region of radius R.
"""
if R is None:
R = self.cosmo.hubble_distance
S = self.schwartzchild()
amp = 3 * S / (2 * R)
return amp
[docs]
def accretion(self):
""" Accretion rate """
bb = self.bondi()
aa = 2 * bb * bb * (self.n/u.m**3) * (
2 * math.pi * c.k_B * self.T * u.K * c.m_p)**0.5
return aa.decompose()
[docs]
def bondi(self):
""" The Bondi Sphere
Equate the root mean square velocity of hydrogen atoms in the
medium, sqrt(3kT/m_h) with the escape velocity sqrt(2GM/B).
"""
# convert central mass to kg.
# Mcent is actually schwartzchild radius in light years
M = self.lightyear_to_kg(self.Mcent)
B = 2 * c.G * M * c.m_p / (3 * c.k_B * self.T * u.K)
#to(u.m*83/u.s**2, equivalencies=u.mass_energy())}")
print('Bondi Sphere:', B.decompose() << u.lightyear)
return B.decompose() << u.lightyear
[docs]
def eddington(self, z=None):
""" Eddington Sphere
Easiest to calculate if we know the redshift.
"""
z = z or self.z_calc()
return self.schwartzchild() / (1 - (1+z)**-2)
[docs]
def schwartzchild(self):
""" Schwartzchild radius in lightyears """
return self.Mcent * u.lightyear
def lightyear_to_kg(self, m=None):
m = m or self.Mcent
# self.Mcent = 2 * M * c.G / (c.c * c.c)
mass = m * u.lightyear * c.c * c.c / (2 * c.G)
return mass.decompose()
[docs]
def z_calc(self):
""" 1.27e7 theta / M
M mass in suns
"""
m = self.lightyear_to_kg() / c.M_sun
theta = (self.T**1.5)/ self.n
return 1.27e7 * theta / m
[docs]
def critical_radius(self):
""" Radius of Mcent based on critical density
Convert Mcent to the radius of the sphere it
would occupy based on critical density.
"""
mass = self.lightyear_to_kg()
cd = self.cosmo.critical_density0 << u.kg/u.m**3
# M / ((4/3) pi r**3) = cd
r = (3 * mass / (4*math.pi*cd))**(1/3)
return r << u.lightyear
[docs]
def tofu(self, u):
""" Work in natural units, c=G=Hubble Distance=1
Time in Hubble times.
"""
return tofu(u, self.theta, self.phi)
def e2tofu(self, u):
return e2tofu(u, self.theta, self.phi)
def tstar(self):
a = cosh(self.phi)
d = cos(self.theta)
# work with U = e**u and T=e**t
A = D = (a - d)/2
B = C = (a + d)/2
# first time visible is tstar given by
tstar = log(sqrt(A/B))
return tstar
[docs]
def umax(self):
""" the last time u that the source can be seen
notice e^umax = 1/(e**tstar = e**(-tstar)
so umax = -1 * tstar.
t for umax is infinite and u for tstar is -infinity
"""
return -1 * self.tstar()
[docs]
def tb(self):
""" blue shift period """
a = cosh(self.phi)
d = cos(self.theta)
tb = log((sqrt(a+1) + sqrt(1-d))/sqrt(a-d))
return tb
def zandx(self, t=None, u=None):
if u is None: u = self.uoft(t)
if t is None: t = self.tofu(u)
return zandx(t, u, self.theta, self.phi)
def uoft(self, t):
return uoft(t, self.theta, self.phi)
def e2uoft(self, t):
return e2uoft(t, self.theta, self.phi)
def rho(self, t):
a = cosh(self.phi)
d = cos(self.theta)
u = self.uoft(t)
return a * sinh(t) * cosh(u) - d * cosh(t) * sinh(u)
async def wits(self):
rr = np.linspace(self.rmin, self.rmax, 100)
nn = 100
ax = await self.get()
ax.projection('polar')
data = []
for r in rr:
r, self.v(r), random.random() * pi
#vv = [self.v(r) for r in rr]
vv = self.v(rr)
async def run(self):
self.spheres3()
rr = np.linspace(self.rmin, self.rmax, 100)
#vv = [self.v(r) for r in rr]
vv = self.v(rr)
ii = self.vinert(rr, vv)
rdd = self.rdoubledot(rr, ii)
energy = np.array([self.energy(r) for r in rr])
#ii = [self.vinert(r, v) for (r, v) in zip(rr, vv)]
#rdd = [self.rdoubledot(r, v) for (r, v) in zip(rr, ii)]
rdot = np.sqrt(2 * energy)
print('energy', max(energy), min(energy))
#print('spiral', len(rr), len(rdot))
if self.details:
ax = await self.get()
print(min(vv), max(vv),
statistics.mean(vv), statistics.variance(vv))
ax.plot(rr, vv, label='velocity')
ax.plot(rr, ii, label='vinert')
ax.plot(rr, rdot, label='rdot')
#ax.plot(rr, energy, label='energy')
ax.legend(loc=0)
ax.show()
ax = await self.get()
ax.plot(rr, rdd, label='rdoubledot')
ax.legend(loc=0)
ax.show()
await self.wits()
thetadot = vv/rr;
dthetabydr = thetadot/rdot
dtbydr = 1/rdot
NIntegrate = integrate.cumtrapz
thetaValues = NIntegrate(dthetabydr, rr, initial=0.)
tvalues = NIntegrate(dtbydr, rr, initial=0.)
B = self.B
ax = await self.get()
# fixme - want polar projections - wibni this worked?
ax.projection('polar')
ax.plot(thetaValues - (B * tvalues), rr)
ax.plot(thetaValues - (B * tvalues) + math.pi, rr)
ax.axis('off')
ax.show()
ax = await self.get()
ax.set_title('r v t')
ax.plot(rr, tvalues)
ax.show()
def zandx(t, u, theta, phi):
a = cosh(phi)
d = cos(theta)
# work with U = e**u and T=e**t
U = e ** u
T = e ** t
A = D = (a - d)/2
B = C = (a + d)/2
z = ((d * tanh(u) - a * tanh(t)) / (a * tanh(u) - d * tanh(t))) - 1
try:
sinhu = sinh(u)
except OverflowError:
#print("sinh overflow for", u)
sinhu = math.log(sys.float_info.max)
try:
coshu = cosh(u)
except OverflowError:
#print("cosh overflow for", u)
coshu = math.log(sys.float_info.max)
x = 1 - (((B/T) - (A*T)) * U)
return z, x
[docs]
def e2tofu(u, theta, phi):
""" Work in natural units, c=G=Hubble Distance=1
Time in Hubble times.
"""
a = cosh(phi)
d = cos(theta)
# work with U = e**u and T=e**t
A = D = (a - d)/2
B = C = (a + d)/2
# first time visible is tstar given by
tstar = log(sqrt(A/B))
# the last time u that the source can be seen
# notice e^umax = 1/(e**tstar = e**(-tstar)
# so umax = -1 * tstar.
# t for umax is infinite and u for tstar is -infinity
umax = log(sqrt(B/A))
U = e**u
T = ((U + sqrt(A*B + ((1 - B*B - A*A) * U*U) + A * B* U*U*U*U))
/ (B - A * U*U))
return T
def tofu(u, theta, phi):
T = e2tofu(u, theta, phi)
t = log(T)
# z = (d * tanh(u) - a * tanh(t)) / (a * tanh(u) - d * tanh(t))
# use distance as t0, adjust time for t0
#hd = self.cosmo.hubble_distance
#t0 = [(-1 * ball.distance/hd) for ball in self.balls]
return t
[docs]
def e2uoft(t, theta, phi):
""" Work in natural units, c=G=Hubble Distance=1
Time in Hubble times.
"""
a = cosh(phi)
d = cos(theta)
# work with U = e**u and T=e**t
A = D = (a - d)/2
B = C = (a + d)/2
# first time visible is tstar given by
ustar = log(sqrt(A/C))
# the last time u that the source can be seen
# notice e^umax = 1/(e**tstar = e**(-tstar)
# so umax = -1 * tstar.
# t for umax is infinite and u for tstar is -infinity
tmax = log(sqrt(C/A))
T = e**t
U = ((T - sqrt(C*D + ((1 - B*C - A*D) * T*T) + A * B * T*T*T*T))
/ (C - A * T*T))
return U
def uoft(t, theta, phi):
#print('ABCDU', A,B,C,D,U)
U = e2uoft(t, theta, phi)
try:
u = log(U)
except:
u = -sys.float_info.max
# z = (d * tanh(u) - a * tanh(t)) / (a * tanh(u) - d * tanh(t))
# use distance as t0, adjust time for t0
#hd = self.cosmo.hubble_distance
#t0 = [(-1 * ball.distance/hd) for ball in self.balls]
return u
[docs]
class Sun(Spiral):
def __init__(self):
super().__init__()
solar_angular_velocity = 2 * math.pi * 365 / 25 # radians per year
# Central mass. Mass converted to Schwartzschild radius (in light years)
# Mass of 1 is approximately 3e12 solar masses.
self.Mcent = (schwartzchild(c.M_sun) << u.lightyear).value
self.Mball = 0.
self.Mdisc = 0.
self.omega0 = solar_angular_velocity # radians per year
# astronomical unit in light years
# au = 1 / 63241.08 ### can't remember how I calculated this
auasly = u.au.to(u.lightyear)
self.rmin = 0.1 * auasly
self.rmax = 50 * auasly
self.K = self.Mcent
# solar wind goes from 30 km/s at 3 AU to 500 km/s at 40 AU
# so set A to 2 * 500 km/s in our units
# ??
# A = K * \omega_0. K = M for Sciama principle
# want 2 * A to be the asymptotic tangential velocity
self.A = self.K * solar_angular_velocity
# or go with asymptotic tangential velocity of 0.6km/h
self.A = (((60 / 3600) * u.m/u.s) / c.c).value
# magic constant determined by overall energy in the orbit
self.EE = 5000
# constant, can be read from tangential velocity for small r
self.find_cc(tangential_velocity=self.rmin * self.omega0)
#self.CC = -0.1
# Apparent rate of precession of the roots of the spiral.
self.B = self.A / self.rmin
self.omega0 = self.A / self.K # angular velocity in radians per year
[docs]
class SpanishDancer(Spiral):
""" Spanish Dancer """
def __init__(self):
super().__init__()
# TODO: import values from sd module?
[docs]
class LLPegasi(Spiral):
""" Beautiful spiral with an 810 year period.
Dust around a binary system showing a clear spiral structure.
The gap between bands and their rate of expansion are consistent with
the pair's 810 year orbit based on their angular separation.
Distance is 1300 * u.pc
temperature 1800K
"""
def __init__(self):
super().__init__()
def pick(x, v, vmin, vmax):
n = len(v)
loc = n * (x - vmin) / vmax
return v[loc]
[docs]
def cpr():
""" Started as Mathematica code from the new paradigm.
adapted to python over time.
See spiral class for more information over time.
"""
Plot = plt.plot
Log = np.log
Sqrt = np.sqrt
NIntegrate = integrate.cumtrapz
A = 0.0005; Mcent = .03; EE = -.00000345; CC = -10;
B = .00000015; Mball = 0; Mdisc = 0; K = Mcent;
rmin = 5000; rmax = 50000; iterate = 1000;
step = (rmax - rmin)/(iterate - 1)
r = np.arange(rmin, rmax)
ax = plt.subplot(121)
v = 2*A - 2*K*A*Log(1 + r/K)/r + CC/r
inert = v - A*r/(K + r);
ax.plot(r, v, label='velocity')
ax.plot(r, inert, label='vinert')
rdoubledot = inert**2/r - Mcent/r**2 - Mdisc/rmax**2 - Mball*r/rmax**3
ax.plot(r, rdoubledot, label='rdoubledot')
energy = (-CC**2/(2*r**2) + (Mcent - 2*A*CC)/r -
Mdisc*r/rmax**2 +
Mball*r**2/(2*rmax**3) +
A**2*K/(K + r) +
A**2*Log(K + r) +
2 * A*K * (CC + 2*A*r) * Log(1 + r/K)/(r**2)
- (2 * A*K*Log(1 + r/K)/r)**2 + EE);
#Plot(energy, r, label='energy')
rdot = Sqrt(2*energy)
ax.plot(r, rdot, label='rdot')
ax.legend(loc=0)
thetadot = v/r;
dthetabydr = thetadot/rdot
dtbydr = 1/rdot
thetaValues = NIntegrate(dthetabydr, r, initial=0.)
print(thetaValues)
print(len(thetaValues))
tvalues = NIntegrate(dtbydr, r, initial=0.)
#thetavalues = Table(
# NIntegrate(dthetabydr, rmin, rmax), ivalue, i, iterate))
#tvalues = Table(
# NIntegrate(dtbydr, r, ivalue, i, iterate))
#ListPolarPlot[{ Table[{thetavalues[[i]] - B*tvalues[[i]], ivalue},
#{i, iterate}] ,
#Table[{thetavalues[[i]] - B*tvalues[[i]] + Pi, ivalue},
#{i, iterate}] }]
print('theta', thetaValues[:5])
ax = plt.subplot(122, projection='polar')
ax.plot(thetaValues - (B * tvalues), r)
ax.plot(thetaValues - (B * tvalues) + math.pi, r)
values = (thetaValues - (B * tvalues))
print(min(values), max(values))
return rdot, inert, v, values
[docs]
def near_galaxies():
""" parse galaxy.txt from
https://heasarc.gsfc.nasa.gov/w3browse/all/neargalcat.html
"""
from astroquery import heasarc
from astropy import coordinates, table
cached = Path('~/karmapi/heasarc.fits').expanduser()
if cached.exists():
table = table.Table.read(cached)
else:
heas = heasarc.Heasarc()
coord = coordinates.SkyCoord('12h00m00 + 0d0m0s', frame='gcrs')
table = heas.query_region(
fields='All',
mission='heasarc_neargalcat',
position=coord,
radius='360 deg')
# write the cache if karmapi exists
if cached.parent.exists():
table.write(cached)
for item in table:
yield item
[docs]
def schwartzchild(m):
""" Mass in kg to schwartzchild radius in lightyears """
return 2 * m * c.G / (c.c * c.c)
[docs]
def from_heasarc(kwargs):
""" Create a Spiral from a Heasarc record """
galaxy = Spiral()
# set mass to M26 mass
m26 = 10 ** kwargs['LOG_MASS_26']
mh1 = 10 ** kwargs['LOG_H1_MASS']
print(m26)
galaxy.Mstellar = m26
galaxy.Mh1 = mh1
galaxy.distance = kwargs['DISTANCE'] * u.Mpc << u.lightyear
galaxy.dec = kwargs['DEC']
galaxy.ra = kwargs['RA']
galaxy.name = kwargs['NAME']
return galaxy
def sample_galaxies(n=1000, fudge=42, t0=0, delta_t=0, min_phi=None, cosmo=None):
# distribution of stellar mass based on heasarc catalog
cosmo = cosmo or Cosmo()
hd = cosmo.hubble_distance
stellar = statistics.NormalDist(
mu=8.459767827529022,
sigma=1.358427349161494)
#h1 = NormalDist(mu=6.368135788262371, sigma=2.9859747769592895)
min_phi = min_phi or sqrt(fudge/2)
random_phi = RandomPhi(max_phi=sqrt(fudge), min_phi=min_phi)
for mass in stellar.samples(n):
galaxy = Spiral()
galaxy.Mstellar = 10**mass
# distance will also get used as t0: time of first contact
galaxy.distance = random.random() * hd * math.sqrt(fudge)
galaxy.dec = (random.random() * 180.0) - 90.0
galaxy.ra = random.random() * 360.0
# magic to get theta distributed as sin(theta)
galaxy.theta = math.acos((2*random.random()-1))
# want phi distributed as sinh(phi) ** 2
galaxy.phi = random_phi()
galaxy.name = f'galaxy{n}'
galaxy.t0 = t0 + random.random() * delta_t
galaxy.xdistance = sqrt(fudge)
yield galaxy
class RandomPhi:
def __init__(self, max_phi=5, min_phi=1e-3, n=10000):
""" return random phi distributed as sinh(phi)**2
How to do this?
scipy.integrate and random.sample
"""
self.max_phi = max_phi or math.sqrt(42)
self.min_phi = min_phi or 1/1000
self.n = n or 10000
self.start()
def start(self):
def sinh2(x):
return math.sinh(x) ** 2
sample = []
values = np.linspace(self.min_phi, self.max_phi, self.n)
self.values = list(values)
lastx = values[0]
for x in values[1:]:
result = integrate.quad(sinh2, lastx, x)
sample.append(result[0])
lastx = x
zval = sample[0]
self.counts = [int(s/zval) for s in sample]
def __call__(self):
return random.sample(self.values[1:], 1, counts=self.counts)[0]
def parse_radec(value):
d, m, s = [float(s) for s in value.split()]
scale = 1
if d < 0:
d *= -1
scale = -1
d += m / 60.
d += s / 3600.
return d * scale
def cleanse(data):
clean = {}
for key in data.keys():
value = data[key]
try:
value = float(value)
except:
pass
if key.lower() in ('ra', 'dec'):
print(key, value)
value = parse_radec(value)
clean[key] = value
return clean
async def run(args):
farm = fm.Farm()
spiral = Spiral()
farm.add(spiral)
if args.skymap:
galaxies = None
if args.heasarc:
galaxies = list(from_heasarc(x) for x in near_galaxies())
skymap = SkyMap(galaxies, n=args.n, fudge=args.hitchhiker)
farm.add(skymap)
await farm.start()
print('about to run farm')
await farm.run()
def main(args=None):
parser = argparse.ArgumentParser()
parser.add_argument('--skymap', action='store_true')
parser.add_argument('--hitchhiker', type=float, default=42)
parser.add_argument('--heasarc', action='store_true')
parser.add_argument('-n', type=int, default=1000)
args = parser.parse_args(args)
magic.run(run(args))
if __name__ == '__main__':
#cpr()
#plt.show()
main()