Skip to content

Commit

Permalink
complete example script.
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbwhit
jbwhit committed Feb 1, 2013
1 parent 3b3a2de commit 62d0def
Show file tree
Hide file tree
Showing 2 changed files with 148 additions and 43 deletions.
115 changes: 99 additions & 16 deletions dipole_error.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,23 @@
#!/usr/bin/env python
"""Calculates the expected $\Delta \alpha/\alpha$ from the King, et al. (2012) and error estimate. Section 5.3."""
import angles
import uncertainties
from uncertainties.umath import *
"""Calculates the expected $\Delta \alpha/\alpha$ from the King, et al. (2012) and error estimate. Section 5.3.
King, J.A., Webb, J.K., Murphy, M.T., Flambaum, V.V., Carswell, R.F., Bainbridge, M.B., Wilczynska, M.R., Koch, F.E.,
Spatial variation in the fine-structure constant - new results from VLT/UVES,
Monthly Notices of the Royal Astronomical Society, 422, 3370-3414. (2012)
"""
try:
import angles
except:
print """python package ''angles'' is required. Try running:
$ sudo pip install -U angles"""
raise
try:
import uncertainties
except:
print """python package ''uncertainties'' is required. Try running:
$ sudo pip install -U uncertainties"""
raise
from uncertainties import umath

# Define exceptions
class DipoleError(Exception):
Expand Down Expand Up @@ -58,7 +73,7 @@ def basic_dipole_monopole(amplitude=DIP_AMPLITUDE, \
:returns: value of predicted dipole at a theta radians away from dipole.
:rtype: number
"""
return amplitude * uncertainties.umath.cos(theta) + monopole
return amplitude * umath.cos(theta) + monopole

wrap_dipole_monopole = uncertainties.wrap(basic_dipole_monopole)

Expand All @@ -84,10 +99,14 @@ def dipole_monopole(right_ascension=QSO_RA, \
:type amplitude: uncertainties.AffineScalarFunc
:param monopole: Monopole term (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type monopole: uncertainties.AffineScalarFunc
returns
"""
pos1 = angles.AngularPosition(alpha=right_ascension, delta=declination)
return wrap_dipole_monopole(amplitude, \
wrap_sep(wrap_radian_RA(dipole_ra), wrap_radian_DEC(dipole_dec), pos1.alpha.r, pos1.delta.r), \
wrap_sep(wrap_radian_RA(dipole_ra), \
wrap_radian_DEC(dipole_dec), \
pos1.alpha.r, \
pos1.delta.r), \
monopole)

# ======================
Expand Down Expand Up @@ -119,8 +138,22 @@ def basic_z_dipole_monopole(prefactor=Z_DIP_PREFACTOR, \
monopole=Z_DIP_MONOPOLE,\
*args,\
**kwargs):
"""docstring for basic_z_dipole_monopole"""
return prefactor * z_redshift ** beta * uncertainties.umath.cos(theta) + monopole
"""
Arguments:
:param prefactor: Amplitude of dipole.
:type prefactor: number
:param z_redshift: Redshift of absorber in sky.
:type z_redshift: number
:param beta: power law exponent.
:type beta: number
:param theta: angle in radians between the dipole and a positions on the sky.
:type theta: number
:param monopole: monopole term.
:type monopole: number
:returns: value of predicted dipole at a theta radians away from dipole.
:rtype: number
"""
return prefactor * z_redshift ** beta * umath.cos(theta) + monopole

wrap_z_dipole_monopole = uncertainties.wrap(basic_z_dipole_monopole)

Expand All @@ -132,13 +165,32 @@ def z_dipole_monopole(right_ascension=QSO_RA, \
z_redshift=REDSHIFT, \
beta=Z_DIP_BETA, \
monopole=Z_DIP_MONOPOLE):
"""docstring for z_dipole_monopole"""
"""Returns the predicted value of alpha from the z_dipole equation.
Arguments:
:param right_ascension:right ascension of point in sky under consideration.
:type right_ascension: number
:param declination: declination of point in sky under consideration.
:type declination: number
:param dipole_ra: RA of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type dipole_ra: uncertainties.AffineScalarFunc
:param dipole_dec: DEC of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type dipole_dec: uncertainties.AffineScalarFunc
:param prefactor: Prefactor term of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type prefactor: uncertainties.AffineScalarFunc
:param monopole: Monopole term (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type monopole: uncertainties.AffineScalarFunc
returns
"""
pos1 = angles.AngularPosition(alpha=right_ascension, delta=declination)
return wrap_z_dipole_monopole(prefactor, \
z_redshift, \
beta, \
wrap_sep(wrap_radian_RA(dipole_ra), wrap_radian_DEC(dipole_dec), \
pos1.alpha.r, pos1.delta.r), \
wrap_sep(wrap_radian_RA(dipole_ra), \
wrap_radian_DEC(dipole_dec), \
pos1.alpha.r, \
pos1.delta.r), \
monopole)

# ======================
Expand Down Expand Up @@ -167,8 +219,19 @@ def basic_r_dipole_monopole(amplitude=R_DIPOLE_AMPLITUDE,\
monopole=R_DIPOLE_MONOPOLE,\
*args,\
**kwargs):
"""docstring for basic_r_dipole_monopole"""
return amplitude * radial_distance * uncertainties.umath.cos(theta) + monopole
"""
Arguments:
:param amplitude: Amplitude term of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type amplitude: number
:param radial_distance: Radial distance in GLyr
:type radial_distance: number
:type theta: number
:param monopole: monopole term.
:type monopole: number
returns
"""
return amplitude * radial_distance * umath.cos(theta) + monopole

wrap_r_dipole_monopole = uncertainties.wrap(basic_r_dipole_monopole)

Expand All @@ -179,12 +242,32 @@ def r_dipole_monopole(right_ascension=QSO_RA, \
amplitude=R_DIPOLE_AMPLITUDE, \
radial_distance=RADIAL_DISTANCE, \
monopole=R_DIPOLE_MONOPOLE):
"""docstring for r_dipole_monopole"""
"""docstring for r_dipole_monopole
Arguments:
:param right_ascension:right ascension of point in sky under consideration.
:type right_ascension: number
:param declination: declination of point in sky under consideration.
:type declination: number
:param dipole_ra: RA of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type dipole_ra: uncertainties.AffineScalarFunc
:param dipole_dec: DEC of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type dipole_dec: uncertainties.AffineScalarFunc
:param amplitude: Amplitude term of dipole (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type amplitude: uncertainties.AffineScalarFunc
:param radial_distance: Radial distance in GLyr
:type radial_distance: uncertainties.AffineScalarFunc
:param monopole: Monopole term (optional with uncertainty via uncertainties.ufloat((value, error)) ).
:type monopole: uncertainties.AffineScalarFunc
returns
"""
pos1 = angles.AngularPosition(alpha=right_ascension, delta=declination)
return wrap_r_dipole_monopole(amplitude, \
radial_distance, \
wrap_sep(wrap_radian_RA(dipole_ra), wrap_radian_DEC(dipole_dec), \
pos1.alpha.r, pos1.delta.r), \
wrap_sep(wrap_radian_RA(dipole_ra), \
wrap_radian_DEC(dipole_dec), \
pos1.alpha.r, \
pos1.delta.r), \
monopole)


Expand Down
76 changes: 49 additions & 27 deletions calculation.py → example.py
View file
Original file line number Diff line number Diff line change
@@ -1,45 +1,67 @@
#!/usr/bin/env python
"""Calculates the expected $\Delta \alpha/\alpha$ from the King, et al. (2012) and error estimate. Section 5.3."""
import angles
import numpy as np
"""Example of dipole_error usage"""
import uncertainties
from dipole_error import *
import dipole_error


QSO_RA = "22h20m06.757" # RA
QSO_DEC = "-28d03m23.34" # DEC
THETA = 1.02 # radians
REDSHIFT = 1.5 # Just pick a redshift to start
RADIAL_DISTANCE = 1.3 # GLyr

# King et al. (2012) dipole location
DIP_RA = 17.3
DIP_RA_ERR = 1.0
DIP_DEC = -61.0
DIP_DEC_ERR = 10.0
DIP_AMPLITUDE = 0.97e-5
DIP_AMPLITUDE_ERR = 0.21e-5 # average of asymmetric errors
DIP_MONOPOLE = -0.178e-5
DIP_MONOPOLE_ERR = 0.084e-5

# Values and errors combined for uncertainties package.
DIPOLE_AMPLITUDE = uncertainties.ufloat((DIP_AMPLITUDE, DIP_AMPLITUDE_ERR))
MONOPOLE = uncertainties.ufloat((DIP_MONOPOLE, DIP_MONOPOLE_ERR))
DIPOLE_RA = uncertainties.ufloat((DIP_RA, DIP_RA_ERR))
DIPOLE_DEC = uncertainties.ufloat((DIP_DEC, DIP_DEC_ERR))

QSO_RA = "22h20m06.757" # RA
QSO_DEC = "-28d03m23.34" # DEC
# QSO_decimal_RA = 335.028153
# QSO_decimal_DEC = -28.056483
QSO_decimal_RA = 335.028153
QSO_decimal_DEC = -28.056483
DIP_decimal_RA = angles.h2d(DIP_RA)
DIP_decimal_DEC = angles.h2d(DIP_DEC)

z_dipole_monopole(right_ascension=QSO_RA, declination=QSO_DEC, dipole_ra=Z_DIPOLE_RA, dipole_dec=Z_DIPOLE_DEC, prefactor=Z_DIPOLE_PREFACTOR, z_redshift=REDSHIFT, beta=Z_DIPOLE_BETA, monopole=Z_DIPOLE_MONOPOLE)
print "Default: "
print dipole_error.dipole_monopole()

print
print "Different RA values: "
for ra_value in ["12h20m06.757", 17.2, "17h12m"]:
print "RA:", ra_value, "\n da/a:", dipole_error.dipole_monopole(right_ascension=ra_value, \
declination=QSO_DEC, \
dipole_ra=DIPOLE_RA, \
dipole_dec=DIPOLE_DEC, \
amplitude=DIPOLE_AMPLITUDE, \
monopole=MONOPOLE)

print
print "Different DEC values: "
for dec_value in ["-28d03m23.34", "-61d03m", 15.0]:
print "DEC:", dec_value, "\n da/a:", \
dipole_error.dipole_monopole(right_ascension=QSO_RA, \
declination=dec_value, \
dipole_ra=DIPOLE_RA, \
dipole_dec=DIPOLE_DEC, \
amplitude=DIPOLE_AMPLITUDE, \
monopole=MONOPOLE)

print
print "Different Amplitude/error values: "
for amplitude_error in [0.1e-5, 2.0e-5]:
print "Amplitude error:", amplitude_error, "\n da/a:", \
dipole_error.dipole_monopole(right_ascension=QSO_RA, \
declination=dec_value, \
dipole_ra=DIPOLE_RA, \
dipole_dec=DIPOLE_DEC, \
amplitude=uncertainties.ufloat((DIP_AMPLITUDE, amplitude_error)), \
monopole=MONOPOLE)


# np.degrees(jw_sep(QSO_decimal_RA, QSO_decimal_DEC, DIP_decimal_RA, DIP_decimal_DEC)) # 58.032267492266556

# def sky_position(right_ascension, declination):
# """docstring for qso_position"""
# return angles.AngularPosition(alpha=right_ascension, delta=declination)
#
# wrapped_sky_position = uncertainties.wrap(sky_position)

# QSO_POSITION = sky_position(QSO_RA, QSO_DEC)
# DIPOLE_POSITION = sky_position(DIP_RA, DIP_DEC)
# # print "QSO_POSITION: ", QSO_POSITION
# # print "DIPOLE_POSITION: ", sky_position(DIP_RA, DIP_DEC)
#
# def monte_carlo(value, error):
# """Return a random value around value +/- error"""
# if error == 0:
Expand Down

0 comments on commit 62d0def

Please sign in to comment.