Skip to content
Snippets Groups Projects
Commit 5e81b103 authored by Folkert Nobels's avatar Folkert Nobels
Browse files

Formatting the code using black

parent b341ce73
No related branches found
No related tags found
1 merge request!681Add example of isolated galaxy and change hernquist potential configuration
Hide whitespace changes
Inline Side-by-side
examples/IsolatedGalaxy_dmparticles/angularmomentum.py
+ 65
− 53
View file @ 5e81b103
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
......@@ -7,11 +25,10 @@ from matplotlib.ticker import MaxNLocator
import scipy.optimize as sco
Nmax = 2001
steps = 10
iterarray = np.arange(0,Nmax+1,steps)
iterarray = np.arange(0, Nmax + 1, steps)
Lxtot = np.zeros(len(iterarray))
Lytot = np.zeros(len(iterarray))
Lztot = np.zeros(len(iterarray))
......@@ -19,59 +36,54 @@ Ltot = np.zeros(len(iterarray))
time_array = np.zeros(len(iterarray))
for i in np.arange(0,Nmax+1,steps):
for i in np.arange(0, Nmax + 1, steps):
print(i)
f = h5py.File('output_%04d.hdf5'%i, 'r')
boxsize = f['Header'].attrs['BoxSize']/2.
time_array[int(i/steps)] = f['Header'].attrs['Time']
particles = f['PartType4']
coordinates = particles['Coordinates'][:,:]
velocities = particles['Velocities'][:,:]
masses = particles['Masses'][:]
R = ((coordinates[:,0]-boxsize[0])**2 + (coordinates[:,1]-boxsize[1])**2)**.5
X = np.abs(coordinates[:,0]-boxsize[0])
Y = np.abs(coordinates[:,1]-boxsize[1])
Z = np.abs(coordinates[:,2]-boxsize[2])
vx = velocities[:,0]
vy = velocities[:,1]
vz = velocities[:,2]
Lx = (Y*vz - Z*vy)*masses
Ly = (Z*vx - X*vz)*masses
Lz = (X*vy - Y*vx)*masses
L = (Lx**2 + Ly**2 + Lz**2)**.5
Lxtot[int(i/steps)] = np.sum(Lx)
Lytot[int(i/steps)] = np.sum(Ly)
Lztot[int(i/steps)] = np.sum(Lz)
Ltot[int(i/steps)] = np.sum(L)
time_array[-1]=2.
plt.plot(time_array,Lxtot/Lxtot[0]-1,label='Lx total')
plt.plot(time_array,Lytot/Lytot[0]-1,label='Ly total')
plt.plot(time_array,Lztot/Lztot[0]-1,label='Lz total')
plt.plot(time_array,Ltot/Ltot[0]-1,label='L total')
plt.xlabel('Time')
plt.ylabel('ratio between current and zero angular momentum')
f = h5py.File("output_%04d.hdf5" % i, "r")
boxsize = f["Header"].attrs["BoxSize"] / 2.0
time_array[int(i / steps)] = f["Header"].attrs["Time"]
particles = f["PartType4"]
coordinates = particles["Coordinates"][:, :]
velocities = particles["Velocities"][:, :]
masses = particles["Masses"][:]
R = (
(coordinates[:, 0] - boxsize[0]) ** 2 + (coordinates[:, 1] - boxsize[1]) ** 2
) ** 0.5
X = np.abs(coordinates[:, 0] - boxsize[0])
Y = np.abs(coordinates[:, 1] - boxsize[1])
Z = np.abs(coordinates[:, 2] - boxsize[2])
vx = velocities[:, 0]
vy = velocities[:, 1]
vz = velocities[:, 2]
Lx = (Y * vz - Z * vy) * masses
Ly = (Z * vx - X * vz) * masses
Lz = (X * vy - Y * vx) * masses
L = (Lx ** 2 + Ly ** 2 + Lz ** 2) ** 0.5
Lxtot[int(i / steps)] = np.sum(Lx)
Lytot[int(i / steps)] = np.sum(Ly)
Lztot[int(i / steps)] = np.sum(Lz)
Ltot[int(i / steps)] = np.sum(L)
time_array[-1] = 2.0
plt.plot(time_array, Lxtot / Lxtot[0] - 1, label="Lx total")
plt.plot(time_array, Lytot / Lytot[0] - 1, label="Ly total")
plt.plot(time_array, Lztot / Lztot[0] - 1, label="Lz total")
plt.plot(time_array, Ltot / Ltot[0] - 1, label="L total")
plt.xlabel("Time")
plt.ylabel("ratio between current and zero angular momentum")
plt.legend()
plt.show()
plt.show()
plt.semilogy(time_array, np.absolute(Ltot/Ltot[0]-1))
plt.xlabel('Time (Gyr)')
plt.ylabel('Fractional change of total angular momentum')
plt.savefig('Total_angular_momentum.png')
plt.semilogy(time_array, np.absolute(Ltot / Ltot[0] - 1))
plt.xlabel("Time (Gyr)")
plt.ylabel("Fractional change of total angular momentum")
plt.savefig("Total_angular_momentum.png")
plt.show()
plt.close()
examples/IsolatedGalaxy_dmparticles/profilefit.py
+ 163
− 93
View file @ 5e81b103
#!/usr/bin/env python
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
......@@ -7,40 +25,46 @@ from matplotlib.ticker import MaxNLocator
import scipy.optimize as sco
import os
def linearfunc(x,a,b):
return a*x+b
def radialfunc(r,h,A):
return A * np.exp(-r/h)*r
def linearfunc(x, a, b):
return a * x + b
def verticalfunc(z,A,z0,zoff):
return 2*A*np.exp(-(z-zoff)/z0)
def verticalfunc2(z,A,z0):
return 2*A*np.exp(-(z)/z0)
def verticalfunc3(z,A,z0,zoff,b):
return 2*A*np.exp(-(z-zoff)/z0)+b
def radialfunc(r, h, A):
return A * np.exp(-r / h) * r
def verticalfunc(z, A, z0, zoff):
return 2 * A * np.exp(-(z - zoff) / z0)
def verticalfunc2(z, A, z0):
return 2 * A * np.exp(-(z) / z0)
def verticalfunc3(z, A, z0, zoff, b):
return 2 * A * np.exp(-(z - zoff) / z0) + b
Nmax = 2001
steps = 10
storefits=False
storefits = False
logfit = True
normalfit = False
# if the user wants to store the indivudal fits
if storefits:
if not os.path.exists('radial'):
os.mkdir('radial')
os.mkdir('vertical')
os.mkdir('histsnap')
if not os.path.exists("radial"):
os.mkdir("radial")
os.mkdir("vertical")
os.mkdir("histsnap")
# Initialize the arrays
R_ideal = np.linspace(0,40,100)
Z_ideal = np.linspace(0,10,100)
R_ideal = np.linspace(0, 40, 100)
Z_ideal = np.linspace(0, 10, 100)
iterarray = np.arange(0,Nmax+1,steps)
iterarray = np.arange(0, Nmax + 1, steps)
Z0t = np.zeros(len(iterarray))
Z0terr = np.zeros(len(iterarray))
......@@ -60,112 +84,158 @@ az = np.zeros(len(iterarray))
azerr = np.zeros(len(iterarray))
bz = np.zeros(len(iterarray))
bzerr = np.zeros(len(iterarray))
eps=1e-6
eps = 1e-6
for i in np.arange(0,Nmax+1,steps):
for i in np.arange(0, Nmax + 1, steps):
# Getting the data from the snapshots
f = h5py.File('output_%04d.hdf5'%i, 'r')
f = h5py.File("output_%04d.hdf5" % i, "r")
boxsize = f["Header"].attrs["BoxSize"] / 2.0
boxsize = f['Header'].attrs['BoxSize']/2.
time_array[int(i/steps)] = f['Header'].attrs['Time']
time_array[int(i / steps)] = f["Header"].attrs["Time"]
particles = f['PartType4']
coordinates = particles['Coordinates'][:,:]
masses = particles['Masses'][:]
particles = f["PartType4"]
coordinates = particles["Coordinates"][:, :]
masses = particles["Masses"][:]
R = ((coordinates[:,0]-boxsize[0])**2 + (coordinates[:,1]-boxsize[1])**2)**.5
Z = np.abs(coordinates[:,1]-boxsize[1])
R = (
(coordinates[:, 0] - boxsize[0]) ** 2 + (coordinates[:, 1] - boxsize[1]) ** 2
) ** 0.5
Z = np.abs(coordinates[:, 1] - boxsize[1])
# Bin the coordinates to make them suitable for fitting
Rhist = np.histogram(R,bins=100,range=[0,40],normed=True)
Zhist = np.histogram(Z,bins=100,range=[0,10.0],normed=True)
Rhist = np.histogram(R, bins=100, range=[0, 40], normed=True)
Zhist = np.histogram(Z, bins=100, range=[0, 10.0], normed=True)
# Create correct variables for fitting
Ry = Rhist[0]
Rx = (Rhist[1][1:]+Rhist[1][:len(Rhist[0])])/2.
Rx = (Rhist[1][1:] + Rhist[1][: len(Rhist[0])]) / 2.0
Zy = Zhist[0]
Zx = (Zhist[1][1:]+Zhist[1][:len(Zhist[0])])/2.
Zx = (Zhist[1][1:] + Zhist[1][: len(Zhist[0])]) / 2.0
# Fit with two methods: non-linear LSQ and linear LSQ in log space
bestsolR = sco.curve_fit(radialfunc, Rx[10:], Ry[10:], p0=[2.,.2])
bestsolR = sco.curve_fit(radialfunc, Rx[10:], Ry[10:], p0=[2.0, 0.2])
bestsolZ = sco.curve_fit(verticalfunc, Zx[40:], Zy[40:])
bestsolRlog = sco.curve_fit(linearfunc,Rx[10:],np.log10(Ry[10:]+eps))
bestsolZlog = sco.curve_fit(linearfunc,Zx[40:],np.log10(Zy[40:]+eps))
bestsolRlog = sco.curve_fit(linearfunc, Rx[10:], np.log10(Ry[10:] + eps))
bestsolZlog = sco.curve_fit(linearfunc, Zx[40:], np.log10(Zy[40:] + eps))
# Store variables
h0t[int(i/steps)] = bestsolR[0][0]
Z0t[int(i/steps)] = bestsolZ[0][1]
Ar[int(i/steps)] = bestsolR[0][1]
Az[int(i/steps)] = bestsolZ[0][0]
Z0terr[int(i/steps)] = (bestsolZ[1][1,1])**.5
h0terr[int(i/steps)] = (bestsolR[1][0,0])**.5
Arerr[int(i/steps)] = (bestsolR[1][1,1])**.5
Azerr[int(i/steps)] = (bestsolZ[1][0,0])**.5
ar[int(i/steps)] = bestsolRlog[0][0]
arerr[int(i/steps)] = (bestsolRlog[1][0,0])**.5
br[int(i/steps)] = bestsolRlog[0][1]
brerr[int(i/steps)] = (bestsolRlog[1][1,1])**.5
az[int(i/steps)] = bestsolZlog[0][0]
azerr[int(i/steps)] = (bestsolZlog[1][0,0])**.5
bz[int(i/steps)] = bestsolZlog[0][1]
bzerr[int(i/steps)] = (bestsolZlog[1][1,1])**.5
h0t[int(i / steps)] = bestsolR[0][0]
Z0t[int(i / steps)] = bestsolZ[0][1]
Ar[int(i / steps)] = bestsolR[0][1]
Az[int(i / steps)] = bestsolZ[0][0]
Z0terr[int(i / steps)] = (bestsolZ[1][1, 1]) ** 0.5
h0terr[int(i / steps)] = (bestsolR[1][0, 0]) ** 0.5
Arerr[int(i / steps)] = (bestsolR[1][1, 1]) ** 0.5
Azerr[int(i / steps)] = (bestsolZ[1][0, 0]) ** 0.5
ar[int(i / steps)] = bestsolRlog[0][0]
arerr[int(i / steps)] = (bestsolRlog[1][0, 0]) ** 0.5
br[int(i / steps)] = bestsolRlog[0][1]
brerr[int(i / steps)] = (bestsolRlog[1][1, 1]) ** 0.5
az[int(i / steps)] = bestsolZlog[0][0]
azerr[int(i / steps)] = (bestsolZlog[1][0, 0]) ** 0.5
bz[int(i / steps)] = bestsolZlog[0][1]
bzerr[int(i / steps)] = (bestsolZlog[1][1, 1]) ** 0.5
if storefits:
plt.step(Rx, Ry)
plt.plot(R_ideal,radialfunc(R_ideal,bestsolR[0][0],bestsolR[0][1]),label='Non linear LSQ')
plt.plot(R_ideal,10**(linearfunc(R_ideal,bestsolRlog[0][0],bestsolRlog[0][1])),label='Linear LSQ')
plt.xlim(0,40)
plt.ylim(0,0.25)
plt.xlabel('R (kpc)')
plt.ylabel('Probability')
plt.savefig('./radial/radialsnap%04d.png'%i)
plt.plot(
R_ideal,
radialfunc(R_ideal, bestsolR[0][0], bestsolR[0][1]),
label="Non linear LSQ",
)
plt.plot(
R_ideal,
10 ** (linearfunc(R_ideal, bestsolRlog[0][0], bestsolRlog[0][1])),
label="Linear LSQ",
)
plt.xlim(0, 40)
plt.ylim(0, 0.25)
plt.xlabel("R (kpc)")
plt.ylabel("Probability")
plt.savefig("./radial/radialsnap%04d.png" % i)
plt.close()
plt.step(Zx,Zy)
plt.plot(Z_ideal,verticalfunc(Z_ideal,bestsolZ[0][0],bestsolZ[0][1],bestsolZ[0][2]),label='Non linear LSQ')
plt.plot(Z_ideal,10**(linearfunc(Z_ideal,bestsolZlog[0][0],bestsolZlog[0][1])),label='Linear LSQ')
plt.xlim(0,10.0)
plt.ylim(0,0.6)
plt.xlabel('z (kpc)')
plt.ylabel('Probability')
plt.savefig('./vertical/verticalsnap%04d.png'%i)
plt.step(Zx, Zy)
plt.plot(
Z_ideal,
verticalfunc(Z_ideal, bestsolZ[0][0], bestsolZ[0][1], bestsolZ[0][2]),
label="Non linear LSQ",
)
plt.plot(
Z_ideal,
10 ** (linearfunc(Z_ideal, bestsolZlog[0][0], bestsolZlog[0][1])),
label="Linear LSQ",
)
plt.xlim(0, 10.0)
plt.ylim(0, 0.6)
plt.xlabel("z (kpc)")
plt.ylabel("Probability")
plt.savefig("./vertical/verticalsnap%04d.png" % i)
plt.close()
time_array[-1]=2.
time_array[-1] = 2.0
ax = plt.subplot(111)
ax.set_yscale('log')
ax.set_yscale("log")
if logfit:
plt.errorbar(time_array,np.absolute(az/(az[0])-1), yerr=azerr/(az[0]),label='z0 scale height (Log space)')
plt.errorbar(time_array,np.absolute(ar/(ar[0])-1), yerr=arerr/(ar[0]),label='h scale lenght (Log space)')
plt.errorbar(
time_array,
np.absolute(az / (az[0]) - 1),
yerr=azerr / (az[0]),
label="z0 scale height (Log space)",
)
plt.errorbar(
time_array,
np.absolute(ar / (ar[0]) - 1),
yerr=arerr / (ar[0]),
label="h scale lenght (Log space)",
)
if normalfit:
plt.errorbar(time_array,np.absolute(Z0t/(Z0t[0])-1), yerr=Z0terr/(Z0t[0]),label='z0 scale height (normal space)')
plt.errorbar(time_array,np.absolute(h0t/(h0t[0])-1), yerr=h0terr/(h0t[0]),label='h scale height (normal space)')
ax.set_xlabel('Time (Gyr)')
ax.set_ylabel('Fractional difference')
plt.errorbar(
time_array,
np.absolute(Z0t / (Z0t[0]) - 1),
yerr=Z0terr / (Z0t[0]),
label="z0 scale height (normal space)",
)
plt.errorbar(
time_array,
np.absolute(h0t / (h0t[0]) - 1),
yerr=h0terr / (h0t[0]),
label="h scale height (normal space)",
)
ax.set_xlabel("Time (Gyr)")
ax.set_ylabel("Fractional difference")
plt.legend()
plt.savefig('Fitdifference-witherror.pdf')
plt.savefig("Fitdifference-witherror.pdf")
plt.close()
ax = plt.subplot(111)
ax.set_yscale('log')
ax.set_yscale("log")
if logfit:
plt.plot(time_array,np.absolute(az/(az[0])-1), label='z0 scale height (Log space)')
plt.plot(time_array,np.absolute(ar/(ar[0])-1), label='h scale lenght (Log space)')
plt.plot(
time_array, np.absolute(az / (az[0]) - 1), label="z0 scale height (Log space)"
)
plt.plot(
time_array, np.absolute(ar / (ar[0]) - 1), label="h scale lenght (Log space)"
)
if normalfit:
plt.plot(time_array,np.absolute(Z0t/(Z0t[0])-1),label='z0 scale height (normal space)')
plt.plot(time_array,np.absolute(h0t/(h0t[0])-1),label='h scale height (normal space)')
ax.set_xlabel('Time (Gyr)')
ax.set_ylabel('Fractional difference')
plt.plot(
time_array,
np.absolute(Z0t / (Z0t[0]) - 1),
label="z0 scale height (normal space)",
)
plt.plot(
time_array,
np.absolute(h0t / (h0t[0]) - 1),
label="h scale height (normal space)",
)
ax.set_xlabel("Time (Gyr)")
ax.set_ylabel("Fractional difference")
plt.legend()
plt.savefig('Fitdifference.pdf')
plt.savefig("Fitdifference.pdf")
plt.show()
examples/IsolatedGalaxy_potential/angularmomentum.py
+ 65
− 53
View file @ 5e81b103
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
......@@ -7,11 +25,10 @@ from matplotlib.ticker import MaxNLocator
import scipy.optimize as sco
Nmax = 2001
steps = 10
iterarray = np.arange(0,Nmax+1,steps)
iterarray = np.arange(0, Nmax + 1, steps)
Lxtot = np.zeros(len(iterarray))
Lytot = np.zeros(len(iterarray))
Lztot = np.zeros(len(iterarray))
......@@ -19,59 +36,54 @@ Ltot = np.zeros(len(iterarray))
time_array = np.zeros(len(iterarray))
for i in np.arange(0,Nmax+1,steps):
for i in np.arange(0, Nmax + 1, steps):
print(i)
f = h5py.File('output_%04d.hdf5'%i, 'r')
boxsize = f['Header'].attrs['BoxSize']/2.
time_array[int(i/steps)] = f['Header'].attrs['Time']
particles = f['PartType4']
coordinates = particles['Coordinates'][:,:]
velocities = particles['Velocities'][:,:]
masses = particles['Masses'][:]
R = ((coordinates[:,0]-boxsize[0])**2 + (coordinates[:,1]-boxsize[1])**2)**.5
X = np.abs(coordinates[:,0]-boxsize[0])
Y = np.abs(coordinates[:,1]-boxsize[1])
Z = np.abs(coordinates[:,2]-boxsize[2])
vx = velocities[:,0]
vy = velocities[:,1]
vz = velocities[:,2]
Lx = (Y*vz - Z*vy)*masses
Ly = (Z*vx - X*vz)*masses
Lz = (X*vy - Y*vx)*masses
L = (Lx**2 + Ly**2 + Lz**2)**.5
Lxtot[int(i/steps)] = np.sum(Lx)
Lytot[int(i/steps)] = np.sum(Ly)
Lztot[int(i/steps)] = np.sum(Lz)
Ltot[int(i/steps)] = np.sum(L)
time_array[-1]=2.
plt.plot(time_array,Lxtot/Lxtot[0]-1,label='Lx total')
plt.plot(time_array,Lytot/Lytot[0]-1,label='Ly total')
plt.plot(time_array,Lztot/Lztot[0]-1,label='Lz total')
plt.plot(time_array,Ltot/Ltot[0]-1,label='L total')
plt.xlabel('Time')
plt.ylabel('ratio between current and zero angular momentum')
f = h5py.File("output_%04d.hdf5" % i, "r")
boxsize = f["Header"].attrs["BoxSize"] / 2.0
time_array[int(i / steps)] = f["Header"].attrs["Time"]
particles = f["PartType4"]
coordinates = particles["Coordinates"][:, :]
velocities = particles["Velocities"][:, :]
masses = particles["Masses"][:]
R = (
(coordinates[:, 0] - boxsize[0]) ** 2 + (coordinates[:, 1] - boxsize[1]) ** 2
) ** 0.5
X = np.abs(coordinates[:, 0] - boxsize[0])
Y = np.abs(coordinates[:, 1] - boxsize[1])
Z = np.abs(coordinates[:, 2] - boxsize[2])
vx = velocities[:, 0]
vy = velocities[:, 1]
vz = velocities[:, 2]
Lx = (Y * vz - Z * vy) * masses
Ly = (Z * vx - X * vz) * masses
Lz = (X * vy - Y * vx) * masses
L = (Lx ** 2 + Ly ** 2 + Lz ** 2) ** 0.5
Lxtot[int(i / steps)] = np.sum(Lx)
Lytot[int(i / steps)] = np.sum(Ly)
Lztot[int(i / steps)] = np.sum(Lz)
Ltot[int(i / steps)] = np.sum(L)
time_array[-1] = 2.0
plt.plot(time_array, Lxtot / Lxtot[0] - 1, label="Lx total")
plt.plot(time_array, Lytot / Lytot[0] - 1, label="Ly total")
plt.plot(time_array, Lztot / Lztot[0] - 1, label="Lz total")
plt.plot(time_array, Ltot / Ltot[0] - 1, label="L total")
plt.xlabel("Time")
plt.ylabel("ratio between current and zero angular momentum")
plt.legend()
plt.show()
plt.show()
plt.semilogy(time_array, np.absolute(Ltot/Ltot[0]-1))
plt.xlabel('Time (Gyr)')
plt.ylabel('Fractional change of total angular momentum')
plt.savefig('Total_angular_momentum.png')
plt.semilogy(time_array, np.absolute(Ltot / Ltot[0] - 1))
plt.xlabel("Time (Gyr)")
plt.ylabel("Fractional change of total angular momentum")
plt.savefig("Total_angular_momentum.png")
plt.show()
plt.close()
examples/IsolatedGalaxy_potential/profilefit.py
+ 163
− 93
View file @ 5e81b103
#!/usr/bin/env python
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
......@@ -7,40 +25,46 @@ from matplotlib.ticker import MaxNLocator
import scipy.optimize as sco
import os
def linearfunc(x,a,b):
return a*x+b
def radialfunc(r,h,A):
return A * np.exp(-r/h)*r
def linearfunc(x, a, b):
return a * x + b
def verticalfunc(z,A,z0,zoff):
return 2*A*np.exp(-(z-zoff)/z0)
def verticalfunc2(z,A,z0):
return 2*A*np.exp(-(z)/z0)
def verticalfunc3(z,A,z0,zoff,b):
return 2*A*np.exp(-(z-zoff)/z0)+b
def radialfunc(r, h, A):
return A * np.exp(-r / h) * r
def verticalfunc(z, A, z0, zoff):
return 2 * A * np.exp(-(z - zoff) / z0)
def verticalfunc2(z, A, z0):
return 2 * A * np.exp(-(z) / z0)
def verticalfunc3(z, A, z0, zoff, b):
return 2 * A * np.exp(-(z - zoff) / z0) + b
Nmax = 2001
steps = 10
storefits=False
storefits = False
logfit = True
normalfit = False
# if the user wants to store the indivudal fits
if storefits:
if not os.path.exists('radial'):
os.mkdir('radial')
os.mkdir('vertical')
os.mkdir('histsnap')
if not os.path.exists("radial"):
os.mkdir("radial")
os.mkdir("vertical")
os.mkdir("histsnap")
# Initialize the arrays
R_ideal = np.linspace(0,40,100)
Z_ideal = np.linspace(0,10,100)
R_ideal = np.linspace(0, 40, 100)
Z_ideal = np.linspace(0, 10, 100)
iterarray = np.arange(0,Nmax+1,steps)
iterarray = np.arange(0, Nmax + 1, steps)
Z0t = np.zeros(len(iterarray))
Z0terr = np.zeros(len(iterarray))
......@@ -60,112 +84,158 @@ az = np.zeros(len(iterarray))
azerr = np.zeros(len(iterarray))
bz = np.zeros(len(iterarray))
bzerr = np.zeros(len(iterarray))
eps=1e-6
eps = 1e-6
for i in np.arange(0,Nmax+1,steps):
for i in np.arange(0, Nmax + 1, steps):
# Getting the data from the snapshots
f = h5py.File('output_%04d.hdf5'%i, 'r')
f = h5py.File("output_%04d.hdf5" % i, "r")
boxsize = f["Header"].attrs["BoxSize"] / 2.0
boxsize = f['Header'].attrs['BoxSize']/2.
time_array[int(i/steps)] = f['Header'].attrs['Time']
time_array[int(i / steps)] = f["Header"].attrs["Time"]
particles = f['PartType4']
coordinates = particles['Coordinates'][:,:]
masses = particles['Masses'][:]
particles = f["PartType4"]
coordinates = particles["Coordinates"][:, :]
masses = particles["Masses"][:]
R = ((coordinates[:,0]-boxsize[0])**2 + (coordinates[:,1]-boxsize[1])**2)**.5
Z = np.abs(coordinates[:,1]-boxsize[1])
R = (
(coordinates[:, 0] - boxsize[0]) ** 2 + (coordinates[:, 1] - boxsize[1]) ** 2
) ** 0.5
Z = np.abs(coordinates[:, 1] - boxsize[1])
# Bin the coordinates to make them suitable for fitting
Rhist = np.histogram(R,bins=100,range=[0,40],normed=True)
Zhist = np.histogram(Z,bins=100,range=[0,10.0],normed=True)
Rhist = np.histogram(R, bins=100, range=[0, 40], normed=True)
Zhist = np.histogram(Z, bins=100, range=[0, 10.0], normed=True)
# Create correct variables for fitting
Ry = Rhist[0]
Rx = (Rhist[1][1:]+Rhist[1][:len(Rhist[0])])/2.
Rx = (Rhist[1][1:] + Rhist[1][: len(Rhist[0])]) / 2.0
Zy = Zhist[0]
Zx = (Zhist[1][1:]+Zhist[1][:len(Zhist[0])])/2.
Zx = (Zhist[1][1:] + Zhist[1][: len(Zhist[0])]) / 2.0
# Fit with two methods: non-linear LSQ and linear LSQ in log space
bestsolR = sco.curve_fit(radialfunc, Rx[10:], Ry[10:], p0=[2.,.2])
bestsolR = sco.curve_fit(radialfunc, Rx[10:], Ry[10:], p0=[2.0, 0.2])
bestsolZ = sco.curve_fit(verticalfunc, Zx[40:], Zy[40:])
bestsolRlog = sco.curve_fit(linearfunc,Rx[10:],np.log10(Ry[10:]+eps))
bestsolZlog = sco.curve_fit(linearfunc,Zx[40:],np.log10(Zy[40:]+eps))
bestsolRlog = sco.curve_fit(linearfunc, Rx[10:], np.log10(Ry[10:] + eps))
bestsolZlog = sco.curve_fit(linearfunc, Zx[40:], np.log10(Zy[40:] + eps))
# Store variables
h0t[int(i/steps)] = bestsolR[0][0]
Z0t[int(i/steps)] = bestsolZ[0][1]
Ar[int(i/steps)] = bestsolR[0][1]
Az[int(i/steps)] = bestsolZ[0][0]
Z0terr[int(i/steps)] = (bestsolZ[1][1,1])**.5
h0terr[int(i/steps)] = (bestsolR[1][0,0])**.5
Arerr[int(i/steps)] = (bestsolR[1][1,1])**.5
Azerr[int(i/steps)] = (bestsolZ[1][0,0])**.5
ar[int(i/steps)] = bestsolRlog[0][0]
arerr[int(i/steps)] = (bestsolRlog[1][0,0])**.5
br[int(i/steps)] = bestsolRlog[0][1]
brerr[int(i/steps)] = (bestsolRlog[1][1,1])**.5
az[int(i/steps)] = bestsolZlog[0][0]
azerr[int(i/steps)] = (bestsolZlog[1][0,0])**.5
bz[int(i/steps)] = bestsolZlog[0][1]
bzerr[int(i/steps)] = (bestsolZlog[1][1,1])**.5
h0t[int(i / steps)] = bestsolR[0][0]
Z0t[int(i / steps)] = bestsolZ[0][1]
Ar[int(i / steps)] = bestsolR[0][1]
Az[int(i / steps)] = bestsolZ[0][0]
Z0terr[int(i / steps)] = (bestsolZ[1][1, 1]) ** 0.5
h0terr[int(i / steps)] = (bestsolR[1][0, 0]) ** 0.5
Arerr[int(i / steps)] = (bestsolR[1][1, 1]) ** 0.5
Azerr[int(i / steps)] = (bestsolZ[1][0, 0]) ** 0.5
ar[int(i / steps)] = bestsolRlog[0][0]
arerr[int(i / steps)] = (bestsolRlog[1][0, 0]) ** 0.5
br[int(i / steps)] = bestsolRlog[0][1]
brerr[int(i / steps)] = (bestsolRlog[1][1, 1]) ** 0.5
az[int(i / steps)] = bestsolZlog[0][0]
azerr[int(i / steps)] = (bestsolZlog[1][0, 0]) ** 0.5
bz[int(i / steps)] = bestsolZlog[0][1]
bzerr[int(i / steps)] = (bestsolZlog[1][1, 1]) ** 0.5
if storefits:
plt.step(Rx, Ry)
plt.plot(R_ideal,radialfunc(R_ideal,bestsolR[0][0],bestsolR[0][1]),label='Non linear LSQ')
plt.plot(R_ideal,10**(linearfunc(R_ideal,bestsolRlog[0][0],bestsolRlog[0][1])),label='Linear LSQ')
plt.xlim(0,40)
plt.ylim(0,0.25)
plt.xlabel('R (kpc)')
plt.ylabel('Probability')
plt.savefig('./radial/radialsnap%04d.png'%i)
plt.plot(
R_ideal,
radialfunc(R_ideal, bestsolR[0][0], bestsolR[0][1]),
label="Non linear LSQ",
)
plt.plot(
R_ideal,
10 ** (linearfunc(R_ideal, bestsolRlog[0][0], bestsolRlog[0][1])),
label="Linear LSQ",
)
plt.xlim(0, 40)
plt.ylim(0, 0.25)
plt.xlabel("R (kpc)")
plt.ylabel("Probability")
plt.savefig("./radial/radialsnap%04d.png" % i)
plt.close()
plt.step(Zx,Zy)
plt.plot(Z_ideal,verticalfunc(Z_ideal,bestsolZ[0][0],bestsolZ[0][1],bestsolZ[0][2]),label='Non linear LSQ')
plt.plot(Z_ideal,10**(linearfunc(Z_ideal,bestsolZlog[0][0],bestsolZlog[0][1])),label='Linear LSQ')
plt.xlim(0,10.0)
plt.ylim(0,0.6)
plt.xlabel('z (kpc)')
plt.ylabel('Probability')
plt.savefig('./vertical/verticalsnap%04d.png'%i)
plt.step(Zx, Zy)
plt.plot(
Z_ideal,
verticalfunc(Z_ideal, bestsolZ[0][0], bestsolZ[0][1], bestsolZ[0][2]),
label="Non linear LSQ",
)
plt.plot(
Z_ideal,
10 ** (linearfunc(Z_ideal, bestsolZlog[0][0], bestsolZlog[0][1])),
label="Linear LSQ",
)
plt.xlim(0, 10.0)
plt.ylim(0, 0.6)
plt.xlabel("z (kpc)")
plt.ylabel("Probability")
plt.savefig("./vertical/verticalsnap%04d.png" % i)
plt.close()
time_array[-1]=2.
time_array[-1] = 2.0
ax = plt.subplot(111)
ax.set_yscale('log')
ax.set_yscale("log")
if logfit:
plt.errorbar(time_array,np.absolute(az/(az[0])-1), yerr=azerr/(az[0]),label='z0 scale height (Log space)')
plt.errorbar(time_array,np.absolute(ar/(ar[0])-1), yerr=arerr/(ar[0]),label='h scale lenght (Log space)')
plt.errorbar(
time_array,
np.absolute(az / (az[0]) - 1),
yerr=azerr / (az[0]),
label="z0 scale height (Log space)",
)
plt.errorbar(
time_array,
np.absolute(ar / (ar[0]) - 1),
yerr=arerr / (ar[0]),
label="h scale lenght (Log space)",
)
if normalfit:
plt.errorbar(time_array,np.absolute(Z0t/(Z0t[0])-1), yerr=Z0terr/(Z0t[0]),label='z0 scale height (normal space)')
plt.errorbar(time_array,np.absolute(h0t/(h0t[0])-1), yerr=h0terr/(h0t[0]),label='h scale height (normal space)')
ax.set_xlabel('Time (Gyr)')
ax.set_ylabel('Fractional difference')
plt.errorbar(
time_array,
np.absolute(Z0t / (Z0t[0]) - 1),
yerr=Z0terr / (Z0t[0]),
label="z0 scale height (normal space)",
)
plt.errorbar(
time_array,
np.absolute(h0t / (h0t[0]) - 1),
yerr=h0terr / (h0t[0]),
label="h scale height (normal space)",
)
ax.set_xlabel("Time (Gyr)")
ax.set_ylabel("Fractional difference")
plt.legend()
plt.savefig('Fitdifference-witherror.pdf')
plt.savefig("Fitdifference-witherror.pdf")
plt.close()
ax = plt.subplot(111)
ax.set_yscale('log')
ax.set_yscale("log")
if logfit:
plt.plot(time_array,np.absolute(az/(az[0])-1), label='z0 scale height (Log space)')
plt.plot(time_array,np.absolute(ar/(ar[0])-1), label='h scale lenght (Log space)')
plt.plot(
time_array, np.absolute(az / (az[0]) - 1), label="z0 scale height (Log space)"
)
plt.plot(
time_array, np.absolute(ar / (ar[0]) - 1), label="h scale lenght (Log space)"
)
if normalfit:
plt.plot(time_array,np.absolute(Z0t/(Z0t[0])-1),label='z0 scale height (normal space)')
plt.plot(time_array,np.absolute(h0t/(h0t[0])-1),label='h scale height (normal space)')
ax.set_xlabel('Time (Gyr)')
ax.set_ylabel('Fractional difference')
plt.plot(
time_array,
np.absolute(Z0t / (Z0t[0]) - 1),
label="z0 scale height (normal space)",
)
plt.plot(
time_array,
np.absolute(h0t / (h0t[0]) - 1),
label="h scale height (normal space)",
)
ax.set_xlabel("Time (Gyr)")
ax.set_ylabel("Fractional difference")
plt.legend()
plt.savefig('Fitdifference.pdf')
plt.savefig("Fitdifference.pdf")
plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment