skypyproject · Lucia-Fonseca · Feb 21, 2022 · Feb 21, 2022 · Feb 21, 2022 · Feb 22, 2022
diff --git a/examples/galaxies/plot_demographics.py b/examples/galaxies/plot_demographics.py
@@ -0,0 +1,239 @@
+"""
+Galaxy Demographics
+===================
+
+In this example we reproduce the results from de la Bella et al. 2021 [1]_.
+In that paper the authors distinguished between active galaxies
+(centrals and satellites) and quiescent galaxies (mass-quenched and satellite-quenched),
+describing the galaxy demographics with a set of continuity equations that are solved analytically.
+Such equations invoke two quenching mechanisms that transform star-forming galaxies into quiescent
+objects: mass quenching and satellite quenching.
+
+They showed that the combination of the two quenching mechanisms produces
+a double Schechter function for the quiescent population.
+They demonstrated that the satellite-quenched galaxies are indeed a subset
+of the active galaxies and that the mass-quenched galaxies have a different
+faint-end slope parameter by one unit. The connection between quenching and
+galaxy populations reduced significantly the parameter space of the simulations.
+
+This example uses de la Bella et al. model implemented in `skypy.galaxies.stellar_mass`
-This example uses de la Bella et al. model implemented in `skypy.galaxies.stellar_mass`
+This example uses de la Bella et al. 2021's model implemented in `skypy.galaxies.stellar_mass`
-This example uses de la Bella et al. model implemented in `skypy.galaxies.stellar_mass`
+This example uses de la Bella et al. 2021's model implemented in `skypy.galaxies.stellar_mass`
+and shows how to sample any type of galaxy population
+from a general Schechter mass function.
+"""
+
+# %%
+# Schechter Parameters
+# --------------------
+#
+# We generate their Figure 4 and use the best-fit values of the blue parameters
-# We generate their Figure 4 and use the best-fit values of the blue parameters
+# We generate Figure 3 in de la Bella et al. 2021 and use the best-fit values of the blue parameters
-# We generate their Figure 4 and use the best-fit values of the blue parameters
+# We generate Figure 3 in de la Bella et al. 2021 and use the best-fit values of the blue parameters
+# in Wiegel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
-# in Wiegel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
+# in Weigel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
-# in Wiegel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
+# in Wiegel et al. 2016 [2]_. We follow [1]_'s procedure
-# in Wiegel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
+# in Weigel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
-# in Wiegel et al. 2016 [2]_. We follow de la Bella et al.'s procedure
+# in Wiegel et al. 2016 [2]_. We follow [1]_'s procedure
+# to compute the fraction of satellite galaxies from [2]
-# to compute the fraction of satellite galaxies from [2]
+# to compute the fraction of satellite galaxies from [2]_
-# to compute the fraction of satellite galaxies from [2]
+# to compute the fraction of satellite galaxies from [2]_
+# and use their fixed value for the fraction of satellite-quenched galaxies
+# :math:`f_{\rho} = 0.5`.
+# 
+# Finally we generate the different galaxy samples by feeding the Schechter parameters
-# Finally we generate the different galaxy samples by feeding the Schechter parameters
+# Finally, we generate the different galaxy samples by feeding the Schechter parameters
-# Finally we generate the different galaxy samples by feeding the Schechter parameters
+# Finally, we generate the different galaxy samples by feeding the Schechter parameters
+# as inputs for the Schechter mass function.
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from astropy.cosmology import FlatLambdaCDM
+from astropy.table import Table
+from astropy.units import Quantity
+from skypy.galaxies import schechter_smf
+# from skypy.galaxies.stellar_mass import (schechter_smf_amplitude_centrals,
+                                  #    schechter_smf_amplitude_satellites,
+                                #    schechter_smf_amplitude_mass_quenched,
+                             #    schechter_smf_amplitude_satellite_quenched
+#                                          )
+
+# %%
+
+
+# Weigel et al. 2016 parameters for the active population
+phiblue = 10**-2.423
+mstarb = 10**10.60
+alphab = -1.21
+blue = (phiblue, alphab, mstarb)
+
+# Fraction of satellite-quenched galaxies and satellites
+frho = 0.5
+fsat = 0.4
+
+# Weigel+16 redshift and cosmology
+z_min, z_max = 0.02, 0.06
+z_range = np.linspace(z_min, z_max, 100)
+cosmology = FlatLambdaCDM(H0=70, Om0=0.3)
+
+# Sky area (SDSS DR7 8423 deg2)
+sky_area = Quantity(2000, "deg2")
+
+# %%
+# First we use the model implemented in `skypy.galaxies.stellar_mass`.
-# First we use the model implemented in `skypy.galaxies.stellar_mass`.
+# First, we use the model implemented in `skypy.galaxies.stellar_mass`.
-# First we use the model implemented in `skypy.galaxies.stellar_mass`.
+# First, we use the model implemented in `skypy.galaxies.stellar_mass`.
+
+
+# To be replaced by the SkyPy function once it's merged
+def schechter_smf_phi_centrals(phi_blue_total, fsatellite):
+    if np.ndim(phi_blue_total) == 1 and np.ndim(fsatellite) == 1:
+        phi_blue_total = phi_blue_total[:, np.newaxis]
+
+    sum_phics = (1 - fsatellite) * (1 - np.log(1 - fsatellite))
+    return (1 - fsatellite) * phi_blue_total / sum_phics
+
+
+def schechter_smf_phi_satellites(phi_centrals, fsatellite):
+    return phi_centrals * np.log(1 / (1 - fsatellite))
+
+
+def schechter_smf_phi_mass_quenched(phi_centrals, phi_satellites):
+    return phi_centrals + phi_satellites
+
+
+def schechter_smf_phi_satellite_quenched(phi_satellites, fenvironment):
+    return - np.log(1 - fenvironment) * phi_satellites
+
+
+# SkyPy amplitudes
+phic = schechter_smf_phi_centrals(phiblue, fsat)
+phis = schechter_smf_phi_satellites(phic, fsat)
+phimq = schechter_smf_phi_mass_quenched(phic, phis)
+phisq = schechter_smf_phi_satellite_quenched(phis, frho)
+
+
+# %%
+# Finally we simulate the galaxy populations.
-# Finally we simulate the galaxy populations.
+# Finally, we simulate the galaxy populations.
-# Finally we simulate the galaxy populations.
+# Finally, we simulate the galaxy populations.
+
+
+# Schechter mass functions
+z_centrals, m_centrals = schechter_smf(z_range, mstarb, phic, alphab, 1e9, 1e12, sky_area, cosmology)
+z_satellites, m_satellites = schechter_smf(z_range, mstarb, phis, alphab, 1e9, 1e12, sky_area, cosmology)
+z_massq, m_mass_quenched = schechter_smf(z_range, mstarb, phimq, alphab + 1, 1e9, 1e12, sky_area, cosmology)
+z_satq, m_satellite_quenched = schechter_smf(z_range, mstarb, phisq, alphab, 1e9, 1e12, sky_area, cosmology)
-z_centrals, m_centrals = schechter_smf(z_range, mstarb, phic, alphab, 1e9, 1e12, sky_area, cosmology)
-z_satellites, m_satellites = schechter_smf(z_range, mstarb, phis, alphab, 1e9, 1e12, sky_area, cosmology)
-z_massq, m_mass_quenched = schechter_smf(z_range, mstarb, phimq, alphab + 1, 1e9, 1e12, sky_area, cosmology)
-z_satq, m_satellite_quenched = schechter_smf(z_range, mstarb, phisq, alphab, 1e9, 1e12, sky_area, cosmology)
+z_centrals, m_centrals = schechter_smf(z_range, mstarb, phic, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_satellites, m_satellites = schechter_smf(z_range, mstarb, phis, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_massq, m_mass_quenched = schechter_smf(z_range, mstarb, phimq, alphab + 1, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_satq, m_satellite_quenched = schechter_smf(z_range, mstarb, phisq, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
-z_centrals, m_centrals = schechter_smf(z_range, mstarb, phic, alphab, 1e9, 1e12, sky_area, cosmology)
-z_satellites, m_satellites = schechter_smf(z_range, mstarb, phis, alphab, 1e9, 1e12, sky_area, cosmology)
-z_massq, m_mass_quenched = schechter_smf(z_range, mstarb, phimq, alphab + 1, 1e9, 1e12, sky_area, cosmology)
-z_satq, m_satellite_quenched = schechter_smf(z_range, mstarb, phisq, alphab, 1e9, 1e12, sky_area, cosmology)
+z_centrals, m_centrals = schechter_smf(z_range, mstarb, phic, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_satellites, m_satellites = schechter_smf(z_range, mstarb, phis, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_massq, m_mass_quenched = schechter_smf(z_range, mstarb, phimq, alphab + 1, 
+                                       1e9, 1e12, sky_area, cosmology)
+z_satq, m_satellite_quenched = schechter_smf(z_range, mstarb, phisq, alphab, 
+                                       1e9, 1e12, sky_area, cosmology)
+
+
+logm_centrals = np.log10(m_centrals)
+logm_satellites = np.log10(m_satellites)
+logm_massq = np.log10(m_mass_quenched)
+logm_satq = np.log10(m_satellite_quenched)
+
+# log Mass bins
+bins = np.linspace(9, 12, 35)
+
+# Sky volume
+dV_dz = (cosmology.differential_comoving_volume(z_range) * sky_area).to_value('Mpc3')
+dV = np.trapz(dV_dz, z_range)
+dlm = (np.max(bins)-np.min(bins)) / (np.size(bins)-1)
+
+# log distribution of masses
+logphi_centrals = np.histogram(logm_centrals, bins=bins)[0] / dV / dlm
+logphi_satellites = np.histogram(logm_satellites, bins=bins)[0] / dV / dlm
+logphi_massq = np.histogram(logm_massq, bins=bins)[0] / dV / dlm
+logphi_satq = np.histogram(logm_satq, bins=bins)[0] / dV / dlm
+
+logphi_active = logphi_centrals + logphi_satellites
+logphi_passive = logphi_massq + logphi_satq
+logphi_total = logphi_active + logphi_passive
+
+
+# %%
+# Validation against SSD DR7 data
+# -------------------------------
+#
+# Here we compare our sampled galaxies and
-# Here we compare our sampled galaxies and
+# Here, we compare our sampled galaxies and
-# Here we compare our sampled galaxies and
+# Here, we compare our sampled galaxies and
+# we validate the model using the results from the best fit to
+# SDSS DR7 data in Weigel et al. (2016) [2].
-# SDSS DR7 data in Weigel et al. (2016) [2].
+# SDSS DR7 data in Weigel et al. (2016) [2]_.
-# SDSS DR7 data in Weigel et al. (2016) [2].
+# SDSS DR7 data in Weigel et al. (2016) [2]_.
+# The authors presented a comprehensive method to determine
+# stellar mass functions and apply it to samples in the local universe,
+# in particular to SDSS DR7 data in the redshift range from 0.02 to 0.06.
+#
+# Their data is presented as logarithmic distribution, therefore,
-# Their data is presented as logarithmic distribution, therefore,
+# Their data is presented as a logarithmic distribution, therefore,
-# Their data is presented as logarithmic distribution, therefore,
+# Their data is presented as a logarithmic distribution, therefore,
+# to perform the comparison we need to apply a conversion factor.
+# This factor allows us to go from :math:`\phi` to a :math:`\log \phi` plot
+# and compare with Weigel et al 2016 best-fit model. 
+
+
+# Load the data
+wtotal = Table.read('weigel16_total.csv', format='csv')
+wred = Table.read('weigel16_quiescent.csv', format='csv')
+wblue = Table.read('weigel16_active.csv', format='csv')
+wcentral = Table.read('weigel16_central.csv', format='csv')
+wsatellite = Table.read('weigel16_satellite.csv', format='csv')
+
+# %%
+
+
+# Plot
+fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(16, 6), sharex=True, sharey=True)
+fig.suptitle('Galaxy Demographics Simulation', fontsize=26)
+
+ax1.plot(wblue['logm'], 10**wblue['logphi'], color='k', label='Weigel+16 active', lw=1)
+ax1.step(bins[:-1], logphi_active, where='post', label='SkyPy active', color='blue', zorder=3, lw=1)
+ax1.step(bins[:-1], logphi_centrals, where='post', label='SkyPy centrals', color='royalblue', zorder=3, lw=1)
+ax1.step(bins[:-1], logphi_satellites, where='post', label='SkyPy satellites', color='cyan', zorder=3, lw=1)
+
+
+ax2.plot(wred['logm'], 10**wred['logphi'], color='k', label='Weigel+16 passive', lw=1)
+ax2.fill_between(wred['logm'], 10**wred['upper_error'], 10**wred['lower_error'], color='salmon', alpha=0.1)
+ax2.step(bins[:-1], logphi_passive, where='post', label='SkyPy passive', color='red', zorder=3, lw=1)
+ax2.step(bins[:-1], logphi_massq, where='post', label='SkyPy mass-quenched', color='coral', zorder=3, lw=1)
+ax2.step(bins[:-1], logphi_satq, where='post', label='SkyPy sat-quenched', color='maroon', zorder=3, lw=1)
+
+ax3.plot(wtotal['logm'], 10**wtotal['logphi'], color='k', label='Weigel+16 total', lw=1)
+ax3.plot(wcentral['logm'], 10**wcentral['logphi'], '--', color='grey', label='Weigel+16 centrals', lw=1)
+ax3.plot(wsatellite['logm'], 10**wsatellite['logphi'], '--', color='grey', label='Weigel+16 satellites', lw=1)
+ax3.fill_between(wtotal['logm'], 10**wtotal['upper_error'], 10**wtotal['lower_error'], color='plum', alpha=0.1)
+ax3.step(bins[:-1], logphi_total, where='post', label='SkyPy total', color='purple', zorder=3, lw=1)
+
+for ax in [ax1, ax2, ax3]:
+    ax.legend(loc='lower left', fontsize='small', frameon=False)
+    ax.set_xlabel(r'Stellar mass, $log (M/M_{\odot})$', fontsize=18)
-    ax.set_xlabel(r'Stellar mass, $log (M/M_{\odot})$', fontsize=18)
+    ax.set_xlabel(r'Stellar mass, $\log (M/M_{\odot})$', fontsize=18)
-    ax.set_xlabel(r'Stellar mass, $log (M/M_{\odot})$', fontsize=18)
+    ax.set_xlabel(r'Stellar mass, $\log (M/M_{\odot})$', fontsize=18)
+    ax.set_xlim((9, 11.9))
+    ax.set_ylim((2e-6,5e-2))
+    ax.set_yscale('log')
+
+
+ax1.set_ylabel(r'$\phi /h^3 Mpc^{-3}$', fontsize=18)
-ax1.set_ylabel(r'$\phi /h^3 Mpc^{-3}$', fontsize=18)
+ax1.set_ylabel(r'$\phi /h^3$\,Mpc$^{-3}$', fontsize=18)
-ax1.set_ylabel(r'$\phi /h^3 Mpc^{-3}$', fontsize=18)
+ax1.set_ylabel(r'$\phi /h^3$\,Mpc$^{-3}$', fontsize=18)
+# plt.savefig('galaxy_simulation.pdf')
+plt.tight_layout()
+plt.show()
+
+
+# %%
+# Sonification
+# ------------
+# The sonification, or transformation of physical data via sound,
+# is becoming increasingly important to make astronomy accessible
+# for those who are visually impaired, and to enhance visualisations
+# and convey information that visualisation alone cannot.
+# In this work [1] the authors also made their main plot available
-# In this work [1] the authors also made their main plot available
+# In [1]_, the authors also made their main plot available
-# In this work [1] the authors also made their main plot available
+# In [1]_, the authors also made their main plot available
+# in sound format using the `STRAUSS`_ software (Trayford 2021) [3]_.
+
+import IPython
+# Active population
+IPython.display.Audio("skypy_sim_active.wav")
+
+# %%
+
+# Mass-quenched galaxies
+IPython.display.Audio("skypy_sim_massq.wav")
+
+# %%
+
+# Satellite-quenched galaxies
+IPython.display.Audio("skypy_sim_satq.wav")
+
+
+# %%
+# References
+# ----------
+#
+# .. [1] de la Bella et al. 2021, Quenching and Galaxy Demographics, arXiv 2112.11110.
+#
+# .. [2] Weigel A. K., Schawinski K., Bruderer C., 2016, Monthly Notices of the
+#   Royal Astronomical Society, 459, 2150
+#
+# .. [3] Trayford J., 2021, james-trayford/strauss: v0.1.0 Pre-release, `doi:10.5281/zenodo.5776280`_.
+#   
+# .. _doi:10.5281/zenodo.5776280: https://doi.org/10.5281/zenodo.5776280
+# 
+# .. _STRAUSS: https://strauss.readthedocs.io/en/latest/
diff --git a/examples/galaxies/skypy_passive.wav b/examples/galaxies/skypy_passive.wav
diff --git a/examples/galaxies/skypy_sim_active.wav b/examples/galaxies/skypy_sim_active.wav
diff --git a/examples/galaxies/skypy_sim_massq.wav b/examples/galaxies/skypy_sim_massq.wav
diff --git a/examples/galaxies/skypy_sim_satq.wav b/examples/galaxies/skypy_sim_satq.wav
diff --git a/examples/galaxies/weigel16_active.csv b/examples/galaxies/weigel16_active.csv
@@ -0,0 +1,52 @@
+logm,logphi,upper_error,lower_error
+9.013996889580094,-1.751662778641784,-1.7369399084154145,-1.757217038189045
+9.065318818040435,-1.7608151892172117,-1.7476177207534134,-1.7694202522896152
+9.139968895800934,-1.774127786417834,-1.7631490841541393,-1.787170381890445
+9.186625194401245,-1.7834562254139612,-1.773456225413961,-1.798456225413961
+9.237947122861586,-1.7963772756380938,-1.7863772756380938,-1.811377275638094
+9.349922239502332,-1.8288046334531025,-1.8188046334531025,-1.8438046334531026
+9.415241057542769,-1.8474828761457271,-1.8374828761457271,-1.863106449157871
+9.573872472783826,-1.8908196287349863,-1.8808196287349863,-1.9100545977962353
+9.629860031104199,-1.9081781284375836,-1.8981781284375836,-1.9281781284375838
+9.755832037325039,-1.9542039980275392,-1.9442039980275392,-1.9742039980275394
+9.8398133748056,-1.9858546585607249,-1.9758546585607246,-2.005854658560725
+9.979782270606531,-2.04031528347236,-2.03031528347236,-2.0603152834723604
+10.049766718506998,-2.07347316212013,-2.06347316212013,-2.094777226682359
+10.175738724727838,-2.1390466983313026,-2.1290466983313023,-2.1639937371274858
+10.250388802488336,-2.1904162473698667,-2.1758329521544244,-2.2108329521544245
+10.381026438569208,-2.283465526179369,-2.272888802488337,-2.303465526179369
+10.455676516329705,-2.3539982226171974,-2.343998222617197,-2.3739982226171974
+10.565318818040435,-2.4661782761725473,-2.456178276172547,-2.485920847998691
+10.614307931570762,-2.532390930670378,-2.522390930670378,-2.5493746418924457
+10.656298600311043,-2.5898083579992996,-2.5798083579992994,-2.60550895499925
+10.758942457231726,-2.7386088396127017,-2.728608839612702,-2.7649380424421808
+10.79860031104199,-2.79840435458787,-2.78840435458787,-2.8281041990668747
+10.847589424572316,-2.8962519440124406,-2.8862519440124403,-2.9224572317262822
+10.964230171073094,-3.164797822706063,-3.1547978227060627,-3.1897978227060633
+11.001555209953343,-3.2599934202655803,-3.24999342026558,-3.291659289388681
+11.043545878693624,-3.375161815892767,-3.3651618158927668,-3.419934088721028
+11.06687402799378,-3.455988298896544,-3.445988298896544,-3.497194697474638
+11.150855365474339,-3.75214858236631,-3.74214858236631,-3.783026677832275
+11.188180404354588,-3.899611197511665,-3.889611197511665,-3.9335614307931586
+11.220839813374806,-4.039066874027994,-4.029066874027994,-4.079875583203733
+11.253499222395023,-4.185854908059646,-4.175854908059646,-4.236147653462629
+11.276827371695179,-4.303319001006313,-4.293319001006314,-4.359484951056629
+11.332814930015552,-4.586227060653192,-4.576227060653192,-4.656449455676522
+11.356143079315707,-4.726009331259725,-4.716009331259724,-4.801223950233287
+11.384136858475895,-4.899372731985495,-4.889372731985494,-4.980578019699337
+11.40746500777605,-5.062452047693113,-5.052452047693113,-5.148649559357189
+11.426127527216174,-5.192915500259203,-5.1829155002592024,-5.283106791083465
+11.444790046656298,-5.323378952825292,-5.313378952825293,-5.417564022809741
+11.468118195956453,-5.492027993779154,-5.4820279937791545,-5.584939346811814
+11.519440124416796,-5.91852255054433,-5.883182348367036,-5.995
+11.542768273716952,-6.0,-5.995,-5.995
+11.5800933125972,-6.0,-5.995,-5.995
+11.608087091757387,-6.0,-5.995,-5.995
+11.636080870917574,-6.0,-5.995,-5.995
+11.673405909797822,-6.0,-5.995,-5.995
+11.706065318818041,-6.0,-5.995,-5.995
+11.724727838258165,-6.0,-5.995,-5.995
+11.74805598755832,-6.0,-5.995,-5.995
+11.766718506998444,-6.0,-5.995,-5.995
+11.780715396578538,-6.0,-5.995,-5.995
+11.808709175738725,-6.0,-5.995,-5.995
diff --git a/examples/galaxies/weigel16_central.csv b/examples/galaxies/weigel16_central.csv
@@ -0,0 +1,52 @@
+logm,logphi,upper_error,lower_error
+9.013996889580094,-1.901453957996769,-1.8045234248788353,-2.017770597738289
+9.065318818040435,-1.9084282085920479,-1.809745866849437,-2.017770597738289
+9.139968895800934,-1.9185761939804045,-1.8173566438795332,-2.017770597738289
+9.186625194401245,-1.9233041223271303,-1.8215202111860989,-2.0202411805566842
+9.237947122861586,-1.9275190501441357,-1.8257351390031042,-2.02445610837369
+9.349922239502332,-1.9367152562903298,-1.8349313451492981,-2.033652314519884
+9.415241057542769,-1.941285305619963,-1.8400319710462265,-2.0387510293810056
+9.573872472783826,-1.948729638681243,-1.8511986854074345,-2.0499177437422134
+9.629860031104199,-1.9528651328326536,-1.8551923497949228,-2.0531731908719033
+9.755832037325039,-1.9649773484286204,-1.8642763325112712,-2.0592290600337533
+9.8398133748056,-1.9711786711822379,-1.870332320988837,-2.0651434225686924
+9.979782270606531,-1.9779074146954043,-1.880425635118113,-2.078601439897544
+10.049766718506998,-1.9834370977518496,-1.8864991072917772,-2.0863337305089877
+10.175738724727838,-1.9952609050752397,-1.8983229146151672,-2.1011137739385646
+10.250388802488336,-2.006340782508102,-1.9093949697002188,-2.1119844017507146
+10.381026438569208,-2.0314643948854285,-1.934518582077545,-2.133967129423679
+10.455676516329705,-2.0516150864132165,-1.9525967494469734,-2.154453865139024
+10.565318818040435,-2.086802807443267,-1.982756519787154,-2.1921557032289236
+10.614307931570762,-2.1025249806694597,-1.9962321618540433,-2.2090012053541983
+10.656298600311043,-2.1246172383428195,-2.0179436199903673,-2.232383575046536
+10.758942457231726,-2.1833779030450664,-2.0767042846926143,-2.2944103380857555
+10.79860031104199,-2.212105873279827,-2.1054104671098126,-2.3236510706393863
+10.847589424572316,-2.2551202677917614,-2.14798633627989,-2.3658294825168618
+10.964230171073094,-2.3736469254314545,-2.2630240975821945,-2.4773785184139365
+11.001555209953343,-2.4140830538368827,-2.3024818187243987,-2.5165248918779795
+11.043545878693624,-2.4743431518967496,-2.3627419167842634,-2.5805559535406077
+11.06687402799378,-2.5114401646702804,-2.399819781817757,-2.618207487398778
+11.150855365474339,-2.649559916255992,-2.537939533403471,-2.756327238984492
+11.188180404354588,-2.717767141282523,-2.611211984586968,-2.824563257787379
+11.220839813374806,-2.791526329880765,-2.6864514919945846,-2.8952644416348776
+11.253499222395023,-2.8778767355095543,-2.7662768054628986,-2.975289651715326
+11.276827371695179,-2.942391759804674,-2.8319465420265875,-3.0421141940389838
+11.332814930015552,-3.0972278181129607,-2.98955390977944,-3.202493095615763
+11.356143079315707,-3.1643132107623333,-3.056839098475677,-3.2723935117850034
+11.384136858475895,-3.2518953137361692,-3.142110654673184,-3.3645976440080827
+11.40746500777605,-3.324880399547697,-3.213170284837771,-3.4414344208606473
+11.426127527216174,-3.4048657327766576,-3.2930013259969915,-3.5166434276997007
+11.444790046656298,-3.4854019994246777,-3.3735375926450115,-3.5921294895329563
+11.468118195956453,-3.5887935091801797,-3.4764797088355044,-3.691388145746976
+11.519440124416796,-3.832501920681778,-3.7167902901228724,-3.9384947548553852
+11.542768273716952,-3.9467760360038087,-3.8273056526447142,-4.053855682021509
+11.5800933125972,-4.141548879010175,-4.008592408869329,-4.248628525027869
+11.608087091757387,-4.3088210558961455,-4.170542654256135,-4.415998710978959
+11.636080870917574,-4.486293796788272,-4.34522819994518,-4.593471451871093
+11.673405909797822,-4.729953212353477,-4.606175961222472,-4.868885879439324
+11.706065318818041,-4.9637477231195515,-4.8341168340423035,-5.114926265519992
+11.724727838258165,-5.116771085496064,-4.970450378772958,-5.259629025528218
+11.74805598755832,-5.323471384715013,-5.15366956352008,-5.446522555966339
+11.766718506998444,-5.476882150930529,-5.315071620084174,-5.607924612530431
+11.780715396578538,-5.593426578386139,-5.436334462516161,-5.728603848101238
+11.808709175738725,-5.8786676506827655,-5.707588220947513,-5.976751605056448