v1.1.0 release

changelog.txt

@@ -1,8 +1,14 @@
-v1.0.1
+v1.1.0
  - edited README.md
  - made changes to docstrings to improve readthedocs
  - added model.atan() and model.a_gal_sph() functions to calculate proper accelerations
  - added convenience.dot and convenience.mags helper functions
+ - added CoxGomezSpiralArm potential to Models
+ - changed UniformAccel Model to UniformAlos
+ - added Uniform3DAccel Model
+ - fixed bug in models.a_gal_sph() by changing np.put -> np.place, which can take conditions instead of indices
+ - fixed bug in fit chi^2 statistic calculation, is now correct
+ - fixed bug in fix_arrays() wrapper that caused things to break when not running the docs generator
 
 v1.0.0
  - put out public repo release on Pypi

docs/conf.py

@@ -57,7 +57,7 @@ def setup(app):
 project = 'peebee'
 copyright = '2024, Tom Donlon'
 author = 'Tom Donlon'
-release = '1.0.0' #TODO: read this from a file
+release = '1.1.0' #TODO: read this from a file
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

peebee/convenience.py

@@ -20,7 +20,7 @@ def fix_arrays(func):
     """:meta private:"""
     @wraps(func)
     def wrapper(*args, **kwargs):
-        self.__doc__ = args[0].__doc__
+        #self.__doc__ = args[0].__doc__ #I removed this because it was breaking things, hopefully @wraps fixes the documentation?
 
         use_array = True
 
@@ -216,7 +216,7 @@ def dm_over_bary_alos(l, b, d, model_bary, model_dm, frame='gal'):
 
 #I can never be bothered to write out a pandas df or whatever so I wrote this instead
 def write_to_csv(path, *args, titles=None):
-	"""text autodoc3"""
+	"""text autodoc"""
 
 	if len(args) == 0: 
 		raise Exception("Have to provide at least one array-like in addition to the file path")

peebee/fitter.py

@@ -188,7 +188,7 @@ def fit_model(l, b, d, alos, alos_err, model, bounds, frame='gal', mode='gd', sc
 		scale = kwargs['scale']
 
 	params = result.x
-	red_chi2 = result.fun/(len(alos) - len(params))
+	red_chi2 = result.fun/(len(alos) - len(params) - 1)
 	aic = 2*len(params) + result.fun
 	if print_out:
 		rss(params, model, data, scale, print_out=print_out)

peebee/models.py

@@ -183,7 +183,7 @@ def atan(self, l, b, d, frame='gal', sun_pos=(r_sun, 0., 0.), angular=True): #in
 		tan_accels = np.sum((accels - los_accels*los_vecs)**2, axis=1)**0.5 #magnitude of tangential accels
 
 		if angular: #kpc/s^2 -> radians/s^2
-			tan_accels /=d 
+			tan_accels /= d #TODO: not sure this is correct because it doesn't include cos(b) for one component
 
 		return tan_accels
 
@@ -218,22 +218,22 @@ def a_gal_sph(self, l, b, d, frame='gal', sun_pos=(r_sun, 0., 0.), angular=True)
 		atan_vec = accels - alos*rhats
 
 		bhats = np.array([np.cos(l)*np.sin(b), -np.sin(l)*np.sin(b), np.cos(b)])
-	    lhats = np.array([np.sin(l), np.cos(l), 0])
+		lhats = np.array([np.sin(l), np.cos(l), 0])
 
-	    atan_b = np.dot(atan_vec, bhats)
-	    atan_l = np.dot(atan_vec, lhats)
+		atan_b = np.dot(atan_vec, bhats)
+		atan_l = np.dot(atan_vec, lhats)
 
-	    #fix polar axis
-	    polar_indx = (b == 90.) | (b == -90.)
-	    np.put(atan_l, polar_indx, 0.)
-	    np.put(atan_b, (b == 90.), -mags(atan_vec[(b == 90.)]))
-	    np.put(atan_b, (b == -90.), mags(atan_vec[(b == -90.)]))
+		#fix polar axis
+		polar_indx = (b == 90.) | (b == -90.)
+		np.place(atan_l, polar_indx, 0.)
+		np.place(atan_b, (b == 90.), -mags(atan_vec[(b == 90.)]))
+		np.place(atan_b, (b == -90.), mags(atan_vec[(b == -90.)]))
 
 		if angular: #kpc/s^2 -> radians/s^2
-			atan_b /=d 
-			atan_l /= d
+			atan_b /= d 
+			atan_l = atan_l / d * np.cos(b)
 
-	    return alos, atan_l, atan_b
+		return alos, atan_l, atan_b
 
 	#model1 + model2 returns a new CompositeModel
 	def __add__(self, model2):
@@ -648,7 +648,7 @@ class SinusoidalDisk(Model):
 	def __init__(self, **kwargs):
 		super().__init__()
 		self.name = 'Anharmonic Disk'
-		self.param_names = ['alpha', 'amp', 'lambda', 'phi']
+		self.param_names = ['alpha', 'amp', 'lam', 'phi']
 		self.param_defaults = [None, None, None, None] #None if required param
 		self._finish_init_model(**kwargs)
 
@@ -790,30 +790,119 @@ def accel(self, x, y, z):
 		return ax, ay, az
 
 #---------------------------------------------------------------------------
-# Uniform Acceleration 
+# Uniform Line-of-Sight Acceleration 
 # useful for things like globular clusters where the MW acceleration field is ~ constant
 #---------------------------------------------------------------------------
-class UniformAccel(Model):
+class UniformAlos(Model):
 
 	def __init__(self, **kwargs):
 		super().__init__()
-		self.name = 'Uniform Acceleration'
+		self.name = 'Uniform Line-of-Sight Acceleration'
 		self.param_names = ['alos']
-		self.param_defaults = [None] #None if required param
+		self.param_defaults = [0.] #None if required param
 		self._finish_init_model(**kwargs)
 
 	def accel(self, x, y, z, **kwargs): #should catch everything?
-		raise NotImplementedError('Uniform Acceleration model has no acc() method, it is meant to be used with alos().')
+		raise NotImplementedError('UniformAlos model has no acc() method, it is meant to be used with alos().')
 
 	def alos(self, l, b, d, **kwargs): #should catch everything?
 		return np.zeros(len(l)) + self.alos
 
+#---------------------------------------------------------------------------
+# Uniform 3D Acceleration 
+# useful for things like globular clusters where the MW acceleration field is ~ constant
+#---------------------------------------------------------------------------
+class Uniform3DAccel(Model):
+
+	def __init__(self, **kwargs):
+		super().__init__()
+		self.name = 'Uniform Acceleration'
+		self.param_names = ['ax', 'ay', 'az']
+		self.param_defaults = [0., 0., 0.] #None if required param
+		self._finish_init_model(**kwargs)
+
+	def accel(self, x, y, z, **kwargs): #should catch everything?
+		ax, ay, az = self.log_corr_params()
+		return np.zeros(len(x)) + ax, np.zeros(len(x)) + ay, np.zeros(len(x)) + az
+
+#-----------------------------------------
+# Cox & Gomez (2002) Spiral Arm Potential
+# Rather than do out all the derivatives I just numerically evaluate the acceleration at each point
+#    this is slower than if I had evaluated things out, but much better for my sanity
+#    might eventually add in functionality for linear combinations eventually, but for now we're just going to use 1 term in the sum (n=1)
+#    if speed is enough of a problem that the linear combination inside the model becomes important we're also better off analytically computing the derivative
+#-----------------------------------------
+class CoxGomezSpiralArm(Model):
+
+	def __init__(self, **kwargs):
+		super().__init__()
+		self.name = 'Cox-Gomez Spiral Arm'
+		self.param_names = ['N', 'alpha', 'rs', 'rh', 'phi0', 'rho0']
+		#N: number of spiral arms
+		#alpha: the pitch angle (in degrees)
+		#rs: scale length, kpc (i.e. double-exp disk)
+		#h: scale height, kpc (i.e. double-exp disk)
+		#r0: fiducial radius, kpc
+		#phi0: fiducial rotation phase (rad)
+		#rho_0: fiducial density at r0 and phi0, Msun/kpc^3
+
+		self.param_defaults = [2., 15., 7., 0.18, 0., 3.123e7] #None if required param; these are from Cox & Gomez (2002) as their fiducial parameters
+		self._finish_init_model(**kwargs)
+
+	def set_number_arms(self, N):
+		#this is a really dumb helper function that helps with optimization. However, it is really symptom of a greater problem that you can't currently lock parameters when fitting.
+		self.params['N'] = N
+
+	#helper function to evlauate the potential, since this is done multiple times below
+	def pot(self, x, y, z):
+		N, alpha, rs, rh, phi0, rho0 = self.log_corr_params()
+
+		r = (x**2 + y**2 + z**2)**0.5
+		phi = np.arctan2(y,-x) #left-handed, in radians
+
+		Kn = N/(r*np.sin(alpha*np.pi/180))
+		Bn = Kn*rh*(1 + 0.4*Kn*rh)
+		Dn = (1 + Kn*rh + 0.3*(Kn*rh)**2)/(1 + 0.3*Kn*rh)
+		gamma = N*(phi - phi0 - np.log(r/r_sun)/np.tan(alpha*np.pi/180)) #something like radians
+
+		out = -4*np.pi*G*rh*rho0 * np.exp(-(r - r_sun)/rs) * (1 / (Kn * Dn)) * np.cos(gamma) * (np.cosh(Kn*z/Bn))**(-Bn)
+
+		return out 
+
+	def accel(self, x, y, z, sun_pos=(r_sun, 0., 0.)): #this is probably pretty slow, fyi
+
+		#compute accelerations of each object
+		dr = 0.001 #kpc
+
+		pot_base    = self.pot(x, y, z)
+		pot_plus_dx = self.pot(x+dr, y,    z   )
+		pot_plus_dy = self.pot(x,    y+dr, z   )
+		pot_plus_dz = self.pot(x,    y,    z+dr)
+
+		dpot_dx = -(pot_plus_dx - pot_base)/dr
+		dpot_dy = -(pot_plus_dy - pot_base)/dr
+		dpot_dz = -(pot_plus_dz - pot_base)/dr
+
+		#compute Solar acceleration and subtract off
+		pot_base_sun    = self.pot(sun_pos[0],    sun_pos[1],    sun_pos[2]   )
+		pot_plus_dx_sun = self.pot(sun_pos[0]+dr, sun_pos[1],    sun_pos[2]   )
+		pot_plus_dy_sun = self.pot(sun_pos[0],    sun_pos[1]+dr, sun_pos[2]   )
+		pot_plus_dz_sun = self.pot(sun_pos[0],    sun_pos[1],    sun_pos[2]+dr)
+
+		dpot_dx_sun = -(pot_plus_dx_sun - pot_base_sun)/dr
+		dpot_dy_sun = -(pot_plus_dy_sun - pot_base_sun)/dr
+		dpot_dz_sun = -(pot_plus_dz_sun - pot_base_sun)/dr
+
+		ax = dpot_dx - dpot_dx_sun
+		ay = dpot_dy - dpot_dy_sun
+		az = dpot_dz - dpot_dz_sun
+
+		return ax, ay, az
+
 #-------------------------------------------------
 # Generic Gala Potential Wrapper
 # allows us to easily alos for any Gala potential
-# TODO: make this more homogenous
-#       i.e. act the same as the other models
-# CANNOT UPDATE PARAMS RIGHT NOW
+# TODO: cannot update params right now
 #-------------------------------------------------
 class GalaPotential(Model):
 
@@ -843,7 +932,7 @@ def accel(self, x, y, z):
 #-------------------------------------------------
 # Generic Galpy Potential Wrapper
 # allows us to easily alos for any Galpy potential
-# TODO: CANNOT UPDATE PARAMS RIGHT NOW
+# TODO: cannot update params right now
 #-------------------------------------------------
 class GalpyPotential(Model):
 

peebee/transforms.py

@@ -137,4 +137,4 @@ def wrapper(*args, **kwargs):
 
 			return func(*args, **kwargs)
 		return wrapper
-	return internal_decorator
+	return internal_decorator

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="peebee",
-    version="1.0.0",
+    version="1.1.0",
     author="Tom Donlon",
     author_email="thomas.donlon@uah.edu",
     description="A python package for the intersection of pulsar accelerations and Galactic structure",