The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def angle_wrap(angle,radians=False):
'''
Wraps the input angle to 360.0 degrees.
if radians is True: input is assumed to be in radians, output is also in
radians
'''
if radians:
wrapped = angle % (2.0*PI)
if wrapped < 0.0:
wrapped = 2.0*PI + wrapped
else:
wrapped = angle % 360.0
if wrapped < 0.0:
wrapped = 360.0 + wrapped
return wrapped Example 2
def great_circle_dist(p1, p2):
"""Return the distance (in km) between two points in
geographical coordinates.
"""
lon0, lat0 = p1
lon1, lat1 = p2
EARTH_R = 6372.8
lat0 = np.radians(float(lat0))
lon0 = np.radians(float(lon0))
lat1 = np.radians(float(lat1))
lon1 = np.radians(float(lon1))
dlon = lon0 - lon1
y = np.sqrt(
(np.cos(lat1) * np.sin(dlon)) ** 2
+ (np.cos(lat0) * np.sin(lat1)
- np.sin(lat0) * np.cos(lat1) * np.cos(dlon)) ** 2)
x = np.sin(lat0) * np.sin(lat1) + \
np.cos(lat0) * np.cos(lat1) * np.cos(dlon)
c = np.arctan2(y, x)
return EARTH_R * c Example 3
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data Example 4
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data Example 5
def plotArc(start_angle, stop_angle, radius, width, **kwargs):
""" write a docstring for this function"""
numsegments = 100
theta = np.radians(np.linspace(start_angle+90, stop_angle+90, numsegments))
centerx = 0
centery = 0
x1 = -np.cos(theta) * (radius)
y1 = np.sin(theta) * (radius)
stack1 = np.column_stack([x1, y1])
x2 = -np.cos(theta) * (radius + width)
y2 = np.sin(theta) * (radius + width)
stack2 = np.column_stack([np.flip(x2, axis=0), np.flip(y2,axis=0)])
#add the first values from the first set to close the polygon
np.append(stack2, [[x1[0],y1[0]]], axis=0)
arcArray = np.concatenate((stack1,stack2), axis=0)
return patches.Polygon(arcArray, True, **kwargs), ((x1, y1), (x2, y2)) Example 6
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 7
def computeRotMatrix(self,Phi=False):
#######################################
# COMPUTE ROTATION MATRIX SUCH THAT m(t) = A*L(t)*A'*Hp
# Default set such that phi1,phi2 = 0 is UXO pointed towards North
if Phi is False:
phi1 = np.radians(self.phi[0])
phi2 = np.radians(self.phi[1])
phi3 = np.radians(self.phi[2])
else:
phi1 = np.radians(Phi[0]) # Roll (CCW)
phi2 = np.radians(Phi[1]) # Inclination (+ve is nose pointing down)
phi3 = np.radians(Phi[2]) # Declination (degrees CW from North)
# A1 = np.r_[np.c_[np.cos(phi1),-np.sin(phi1),0.],np.c_[np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis
# A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
# A3 = np.r_[np.c_[np.cos(phi3),-np.sin(phi3),0.],np.c_[np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis (direction of head of object)
A1 = np.r_[np.c_[np.cos(phi1),np.sin(phi1),0.],np.c_[-np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis
A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
A3 = np.r_[np.c_[np.cos(phi3),np.sin(phi3),0.],np.c_[-np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis (direction of head of object)
return np.dot(A3,np.dot(A2,A1)) Example 8
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 9
def test_projection(self):
"""Tests the electric field projection."""
projection = self.field.projection
# Top-right quadrant
a = radians(45)
self.assertTrue(isclose(projection([0, 0], a), -4*cos(a)))
self.assertTrue(isclose(projection([3, 0], a), 0.375*cos(a)))
self.assertTrue(isclose(projection([0, 1], a), -sqrt(2)*cos(a)))
self.assertTrue(isclose(projection([[0, 0], [3, 0], [0, 1]], a),
array([-4, 0.375, -sqrt(2)])*cos(a)).all())
# Bottom-left quadrant
a1 = radians(-135)
a2 = radians(45)
self.assertTrue(isclose(projection([0, 0], a1), 4*cos(a2)))
self.assertTrue(isclose(projection([3, 0], a1), -0.375*cos(a2)))
self.assertTrue(isclose(projection([0, 1], a1),
sqrt(2)*cos(a2)))
self.assertTrue(isclose(projection([[0, 0], [3, 0], [0, 1]], a1),
array([4, -0.375, sqrt(2)])*cos(a2)).all()) Example 10
def cie_relative_luminance(sky_elevation, sky_azimuth=None, sun_elevation=None,
sun_azimuth=None, type='soc'):
""" cie relative luminance of a sky element relative to the luminance
at zenith
angle in radians
type is one of 'soc' (standard overcast sky), 'uoc' (uniform radiance)
or 'clear_sky' (standard clear sky low turbidity)
"""
if type == 'clear_sky' and (
sun_elevation is None or sun_azimuth is None or sky_azimuth is None):
raise ValueError, 'Clear sky requires sun position'
if type == 'soc':
return cie_luminance_gradation(sky_elevation, 4, -0.7)
elif type == 'uoc':
return cie_luminance_gradation(sky_elevation, 0, -1)
elif type == 'clear_sky':
return cie_luminance_gradation(sky_elevation, -1,
-0.32) * cie_scattering_indicatrix(
sun_azimuth, sun_elevation, sky_azimuth, sky_elevation, 10, -3,
0.45)
else:
raise ValueError, 'Unknown sky type' Example 11
def ecliptic_longitude(hUTC, dayofyear, year):
""" Ecliptic longitude
Args:
hUTC: fractional hour (UTC time)
dayofyear (int):
year (int):
Returns:
(float) the ecliptic longitude (degrees)
Details:
World Meteorological Organization (2006).Guide to meteorological
instruments and methods of observation. Geneva, Switzerland.
"""
jd = julian_date(hUTC, dayofyear, year)
n = jd - 2451545
# mean longitude (deg)
L = numpy.mod(280.46 + 0.9856474 * n, 360)
# mean anomaly (deg)
g = numpy.mod(357.528 + 0.9856003 * n, 360)
return L + 1.915 * numpy.sin(numpy.radians(g)) + 0.02 * numpy.sin(
numpy.radians(2 * g)) Example 12
def actual_sky_irradiances(dates, ghi, dhi=None, Tdew=None, longitude=_longitude, latitude=_latitude, altitude=_altitude, method='dirint'):
""" derive a sky irradiance dataframe from actual weather data"""
df = sun_position(dates=dates, latitude=latitude, longitude=longitude, altitude=altitude, filter_night=False)
df['am'] = air_mass(df['zenith'], altitude)
df['dni_extra'] = sun_extraradiation(df.index)
if dhi is None:
pressure = pvlib.atmosphere.alt2pres(altitude)
dhi = diffuse_horizontal_irradiance(ghi, df['elevation'], dates, pressure=pressure, temp_dew=Tdew, method=method)
df['ghi'] = ghi
df['dhi'] = dhi
el = numpy.radians(df['elevation'])
df['dni'] = (df['ghi'] - df['dhi']) / numpy.sin(el)
df['brightness'] = brightness(df['am'], df['dhi'], df['dni_extra'])
df['clearness'] = clearness(df['dni'], df['dhi'], df['zenith'])
return df.loc[(df['elevation'] > 0) & (df['ghi'] > 0) , ['azimuth', 'zenith', 'elevation', 'am', 'dni_extra', 'clearness', 'brightness', 'ghi', 'dni', 'dhi' ]] Example 13
def rotation_matrix(alpha, beta, gamma):
"""
Return the rotation matrix used to rotate a set of cartesian
coordinates by alpha radians about the z-axis, then beta radians
about the y'-axis and then gamma radians about the z''-axis.
"""
ALPHA = np.array([[np.cos(alpha), -np.sin(alpha), 0],
[np.sin(alpha), np.cos(alpha), 0],
[0, 0, 1]])
BETA = np.array([[np.cos(beta), 0, np.sin(beta)],
[0, 1, 0],
[-np.sin(beta), 0, np.cos(beta)]])
GAMMA = np.array([[np.cos(gamma), -np.sin(gamma), 0],
[np.sin(gamma), np.cos(gamma), 0],
[0, 0, 1]])
R = ALPHA.dot(BETA).dot(GAMMA)
return(R) Example 14
def __init__(self, lat0, lon0, depth0, nlat, nlon, ndepth, dlat, dlon, ddepth):
# NOTE: Origin of spherical coordinate system and geographic coordinate
# system is not the same!
# Geographic coordinate system
self.lat0, self.lon0, self.depth0 =\
seispy.coords.as_geographic([lat0, lon0, depth0])
self.nlat, self.nlon, self.ndepth = nlat, nlon, ndepth
self.dlat, self.dlon, self.ddepth = dlat, dlon, ddepth
# Spherical/Pseudospherical coordinate systems
self.nrho = self.ndepth
self.ntheta = self.nlambda = self.nlat
self.nphi = self.nlon
self.drho = self.ddepth
self.dtheta = self.dlambda = np.radians(self.dlat)
self.dphi = np.radians(self.dlon)
self.rho0 = seispy.constants.EARTH_RADIUS\
- (self.depth0 + (self.ndepth - 1) * self.ddepth)
self.lambda0 = np.radians(self.lat0)
self.theta0 = ?/2 - (self.lambda0 + (self.nlambda - 1) * self.dlambda)
self.phi0 = np.radians(self.lon0) Example 15
def _add_triangular_sides(self, xy_mask, angle, y_top_right, y_bot_left,
x_top_right, x_bot_left, n_material):
angle = np.radians(angle)
trap_len = (y_top_right - y_bot_left) / np.tan(angle)
num_x_iterations = round(trap_len/self.x_step)
y_per_iteration = round(self.y_pts / num_x_iterations)
lhs_x_start_index = int(x_bot_left/ self.x_step + 0.5)
rhs_x_stop_index = int(x_top_right/ self.x_step + 1 + 0.5)
for i, _ in enumerate(xy_mask):
xy_mask[i][:lhs_x_start_index] = False
xy_mask[i][lhs_x_start_index:rhs_x_stop_index] = True
if i % y_per_iteration == 0:
lhs_x_start_index -= 1
rhs_x_stop_index += 1
self.n[xy_mask] = n_material
return self.n Example 16
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 17
def semiMajAmp(m1,m2,inc,ecc,p):
"""
K = [(2*pi*G)/p]^(1/3) [m2*sin(i)/m2^(2/3)*sqrt(1-e^2)]
units:
K [m/s]
m1 [Msun]
m2 [Mj]
p [yrs]
inc [deg]
"""
pSecs = p*sec_per_year
m1KG = m1*const.M_sun.value
A = ((2.0*np.pi*const.G.value)/pSecs)**(1.0/3.0)
B = (m2*const.M_jup.value*np.sin(np.radians(inc)))
C = m1KG**(2.0/3.0)*np.sqrt(1.0-ecc**2.0)
print('Resulting K is '+repr(A*(B/C)))
#return A*(B/C) Example 18
def inc_prior_fn(self, inc):
ret = 1.0
if (self.inc_prior==False) or (self.inc_prior=="uniform"):
ret = self.uniform_fn(inc,7)
else:
if inc not in [0.0,90.0,180.0]:
mn = self.mins_ary[7]
mx = self.maxs_ary[7]
# Account for case where min = -(max), as it causes error in denominator
if mn == -1*mx:
mn = mn-0.1
mn_rad = np.radians(mn)
mx_rad = np.radians(mx)
inc_rad = np.radians(inc)
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
ret = np.abs(np.sin(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
elif self.inc_prior == 'cos':
ret = np.abs(np.cos(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
#if ret==0: ret=-np.inf
return ret Example 19
def inc_prior_fn(self, inc):
ret = 1.0
if (self.inc_prior==False) or (self.inc_prior=="uniform"):
ret = self.uniform_fn(inc,7)
else:
if inc not in [0.0,90.0,180.0]:
mn = self.mins_ary[7]
mx = self.maxs_ary[7]
# Account for case where min = -(max), as it causes error in denominator
if mn == -1*mx:
mn = mn-0.1
mn_rad = np.radians(mn)
mx_rad = np.radians(mx)
inc_rad = np.radians(inc)
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
ret = np.abs(np.sin(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
elif self.inc_prior == 'cos':
ret = np.abs(np.cos(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
#if ret==0: ret=-np.inf
return ret Example 20
def hav_dist(locs1, locs2):
"""
Return a distance matrix between two set of coordinates.
Use geometric distance (default) or haversine distance (if longlat=True).
Parameters
----------
locs1 : numpy.array
The first set of coordinates as [(long, lat), (long, lat)].
locs2 : numpy.array
The second set of coordinates as [(long, lat), (long, lat)].
Returns
-------
mat_dist : numpy.array
The distance matrix between locs1 and locs2
"""
locs1 = np.radians(locs1)
locs2 = np.radians(locs2)
cos_lat1 = np.cos(locs1[..., 0])
cos_lat2 = np.cos(locs2[..., 0])
cos_lat_d = np.cos(locs1[..., 0] - locs2[..., 0])
cos_lon_d = np.cos(locs1[..., 1] - locs2[..., 1])
return 6367000 * np.arccos(
cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d)) Example 21
def _get_rotation_matrix(axis, angle):
"""
Helper function to generate a rotation matrix for an x, y, or z axis
Used in get_major_angles
"""
cos = np.cos
sin = np.sin
angle = np.radians(angle)
if axis == 2:
# z axis
return np.array([[cos(angle), -sin(angle), 0], [sin(angle), cos(angle), 0], [0, 0, 1]])
if axis == 1:
# y axis
return np.array([[cos(angle), 0, sin(angle)], [0, 1, 0], [-sin(angle), 0, cos(angle)]])
else:
# x axis
return np.array([[1, 0, 0], [0, cos(angle), -sin(angle)], [0, sin(angle), cos(angle)]]) Example 22
def euler(xyz, order='xyz', units='deg'):
if not hasattr(xyz, '__iter__'):
xyz = [xyz]
if units == 'deg':
xyz = np.radians(xyz)
r = np.eye(3)
for theta, axis in zip(xyz, order):
c = np.cos(theta)
s = np.sin(theta)
if axis == 'x':
r = np.dot(np.array([[1, 0, 0], [0, c, -s], [0, s, c]]), r)
if axis == 'y':
r = np.dot(np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]]), r)
if axis == 'z':
r = np.dot(np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]]), r)
return r Example 23
def get_q_per_pixel(self):
'''Gets the delta-q associated with a single pixel. This is computed in
the small-angle limit, so it should only be considered approximate.
For instance, wide-angle detectors will have different delta-q across
the detector face.'''
if self.q_per_pixel is not None:
return self.q_per_pixel
c = (self.pixel_size_um/1e6)/self.distance_m
twotheta = np.arctan(c) # radians
self.q_per_pixel = 2.0*self.get_k()*np.sin(twotheta/2.0)
return self.q_per_pixel
# Maps
######################################## Example 24
def add_motion_spikes(motion_mat,frequency,severity,TR): time = motion_mat[:,0] max_translation = 5 * severity / 1000 * np.sqrt(2*np.pi)#Max translations in m, factor of sqrt(2*pi) accounts for normalisation factor in norm.pdf later on max_rotation = np.radians(5 * severity) *np.sqrt(2*np.pi) #Max rotation in rad time_blocks = np.floor(time[-1]/TR).astype(np.int32) for i in range(time_blocks): if np.random.uniform(0,1) < frequency: #Decide whether to add spike for j in range(1,4): if np.random.uniform(0,1) < 1/6: motion_mat[:,j] = motion_mat[:,j] \ + max_translation * random.uniform(-1,1) \ * norm.pdf(time,loc = (i+0.5)*TR,scale = TR/5) for j in range(4,7): if np.random.uniform(0,1) < 1/6: motion_mat[:,j] = motion_mat[:,j] \ + max_rotation * random.uniform(-1,1) \ * norm.pdf(time,loc = (i+0.5 + np.random.uniform(-0.25,-.25))*TR,scale = TR/5) return motion_mat
Example 25
def tiltFactor(self, midpointdepth=None,
printAvAngle=False):
'''
get tilt factor from inverse distance law
https://en.wikipedia.org/wiki/Inverse-square_law
'''
# TODO: can also be only def. with FOV, rot, tilt
beta2 = self.viewAngle(midpointdepth=midpointdepth)
try:
angles, vals = getattr(
emissivity_vs_angle, self.opts['material'])()
except AttributeError:
raise AttributeError("material[%s] is not in list of know materials: %s" % (
self.opts['material'], [o[0] for o in getmembers(emissivity_vs_angle)
if isfunction(o[1])]))
if printAvAngle:
avg_angle = beta2[self.foreground()].mean()
print('angle: %s DEG' % np.degrees(avg_angle))
# use averaged angle instead of beta2 to not overemphasize correction
normEmissivity = np.clip(
InterpolatedUnivariateSpline(
np.radians(angles), vals)(beta2), 0, 1)
return normEmissivity Example 26
def setPose(self, obj_center=None, distance=None,
rotation=None, tilt=None, pitch=None):
tvec, rvec = self.pose()
if distance is not None:
tvec[2, 0] = distance
if obj_center is not None:
tvec[0, 0] = obj_center[0]
tvec[1, 0] = obj_center[1]
if rotation is None and tilt is None:
return rvec
r, t, p = rvec2euler(rvec)
if rotation is not None:
r = np.radians(rotation)
if tilt is not None:
t = np.radians(tilt)
if pitch is not None:
p = np.radians(pitch)
rvec = euler2rvec(r, t, p)
self._pose = tvec, rvec Example 27
def hav(alpha):
""" Formula for haversine
Parameters
----------
alpha : (float)
Angle in radians
Returns
--------
hav_alpha : (float)
Haversine of alpha, equal to the square of the sine of half-alpha
"""
hav_alpha = np.sin(alpha * 0.5)**2
return hav_alpha Example 28
def archav(hav):
""" Formula for the inverse haversine
Parameters
-----------
hav : (float)
Haversine of an angle
Returns
---------
alpha : (float)
Angle in radians
"""
alpha = 2.0 * np.arcsin(np.sqrt(hav))
return alpha Example 29
def spherical_to_cartesian(s,degrees=True,normalize=False):
'''
Takes a vector in spherical coordinates and converts it to cartesian.
Assumes the input vector is given as [radius,colat,lon]
'''
if degrees:
s[1] = np.radians(s[1])
s[2] = np.radians(s[2])
x1 = s[0]*np.sin(s[1])*np.cos(s[2])
x2 = s[0]*np.sin(s[1])*np.sin(s[2])
x3 = s[0]*np.cos(s[1])
x = [x1,x2,x3]
if normalize:
x /= np.linalg.norm(x)
return x Example 30
def rotate_delays(lat_r,lon_r,lon_0=0.0,lat_0=0.0,degrees=0): ''' Rotates the source and receiver of a trace object around an arbitrary axis. ''' alpha = np.radians(degrees) colat_r = 90.0-lat_r colat_0 = 90.0-lat_0 x_r = lon_r - lon_0 y_r = colat_0 - colat_r #rotate receivers lat_rotated = 90.0-colat_0+x_r*np.sin(alpha) + y_r*np.cos(alpha) lon_rotated = lon_0+x_r*np.cos(alpha) - y_r*np.sin(alpha) return lat_rotated, lon_rotated
Example 31
def rotate_delays(lat_r,lon_r,lon_0=0.0,lat_0=0.0,degrees=0): ''' Rotates the source and receiver of a trace object around an arbitrary axis. ''' alpha = np.radians(degrees) colat_r = 90.0-lat_r colat_0 = 90.0-lat_0 x_r = lon_r - lon_0 y_r = colat_0 - colat_r #rotate receivers lat_rotated = 90.0-colat_0+x_r*np.sin(alpha) + y_r*np.cos(alpha) lon_rotated = lon_0+x_r*np.cos(alpha) - y_r*np.sin(alpha) return lat_rotated, lon_rotated
Example 32
def rotate_setup(self,lon_0=60.0,colat_0=90.0,degrees=0):
alpha = np.radians(degrees)
lon_s = self.sy
lon_r = self.ry
colat_s = self.sx
colat_r = self.rx
x_s = lon_s - lon_0
y_s = colat_0 - colat_s
x_r = lon_r - lon_0
y_r = colat_0 - colat_r
#rotate receiver
self.rx = colat_0+x_r*np.sin(alpha) + y_r*np.cos(alpha)
self.ry = lon_0+x_r*np.cos(alpha) - y_r*np.sin(alpha)
#rotate source
self.sx = colat_0+x_s*np.sin(alpha) + y_s*np.cos(alpha)
self.sy = lon_0+x_s*np.cos(alpha) - y_s*np.sin(alpha)
#########################################################################
# Plot map of earthquake and station
######################################################################### Example 33
def _calcArea_(self, v1, v2):
"""
Private method to calculate the area covered by a spherical
quadrilateral with one corner defined by the normal vectors
of the two intersecting great circles.
INPUTS:
v1, v2: float array(3), the normal vectors
RETURNS:
area: float, the area given in square radians
"""
angle = self.calcAngle(v1, v2)
area = (4*angle - 2*np.math.pi)
return area Example 34
def _rotVector_(self, v, angle, axis):
"""
Rotate a vector by an angle around an axis
INPUTS:
v: 3-dim float array
angle: float, the rotation angle in radians
axis: string, 'x', 'y', or 'z'
RETURNS:
float array(3): the rotated vector
"""
# axisd = {'x':[1,0,0], 'y':[0,1,0], 'z':[0,0,1]}
# construct quaternion and rotate...
rot = cgt.quat()
rot.fromAngleAxis(angle, axis)
return list(rot.rotateVec(v)) Example 35
def _geodetic_to_cartesian(lat, lon, alt):
"""Conversion from latitude, longitue and altitude coordinates to
cartesian with respect to an ellipsoid
Args:
lat (float): Latitude in radians
lon (float): Longitue in radians
alt (float): Altitude to sea level in meters
Return:
numpy.array: 3D element (in meters)
"""
C = Earth.r / np.sqrt(1 - (Earth.e * np.sin(lat)) ** 2)
S = Earth.r * (1 - Earth.e ** 2) / np.sqrt(1 - (Earth.e * np.sin(lat)) ** 2)
r_d = (C + alt) * np.cos(lat)
r_k = (S + alt) * np.sin(lat)
norm = np.sqrt(r_d ** 2 + r_k ** 2)
return norm * np.array([
np.cos(lat) * np.cos(lon),
np.cos(lat) * np.sin(lon),
np.sin(lat)
]) Example 36
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 37
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 38
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 39
def set_radmask(self, cluster, mpcscale):
"""
Assign mask (0/1) values to maskgals for a given cluster
parameters
----------
cluster: Cluster object
mpcscale: float
scaling to go from mpc to degrees (check units) at cluster redshift
results
-------
sets maskgals.mark
"""
# note this probably can be in the superclass, no?
ras = cluster.ra + self.maskgals.x/(mpcscale*SEC_PER_DEG)/np.cos(np.radians(cluster.dec))
decs = cluster.dec + self.maskgals.y/(mpcscale*SEC_PER_DEG)
self.maskgals.mark = self.compute_radmask(ras,decs) Example 40
def _calc_bkg_density(self, r, chisq, refmag):
"""
Internal method to compute background filter
parameters
----------
bkg: Background object
background
cosmo: Cosmology object
cosmology scaling info
returns
-------
bcounts: float array
b(x) for the cluster
"""
mpc_scale = np.radians(1.) * self.cosmo.Dl(0, self.z) / (1 + self.z)**2
sigma_g = self.bkg.sigma_g_lookup(self.z, chisq, refmag)
return 2 * np.pi * r * (sigma_g/mpc_scale**2) Example 41
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 42
def apply_grad_cartesian_tensor(grad_X, zmat_dist):
"""Apply the gradient for transformation to cartesian space onto zmat_dist.
Args:
grad_X (:class:`numpy.ndarray`): A ``(3, n, n, 3)`` array.
The mathematical details of the index layout is explained in
:meth:`~chemcoord.Cartesian.get_grad_zmat()`.
zmat_dist (:class:`~chemcoord.Zmat`):
Distortions in Zmatrix space.
Returns:
:class:`~chemcoord.Cartesian`: Distortions in cartesian space.
"""
columns = ['bond', 'angle', 'dihedral']
C_dist = zmat_dist.loc[:, columns].values.T
try:
C_dist = C_dist.astype('f8')
C_dist[[1, 2], :] = np.radians(C_dist[[1, 2], :])
except (TypeError, AttributeError):
C_dist[[1, 2], :] = sympy.rad(C_dist[[1, 2], :])
cart_dist = np.tensordot(grad_X, C_dist, axes=([3, 2], [0, 1])).T
from chemcoord.cartesian_coordinates.cartesian_class_main import Cartesian
return Cartesian(atoms=zmat_dist['atom'],
coords=cart_dist, index=zmat_dist.index) Example 43
def drawSkymapCatalog(ax,lon,lat,**kwargs):
mapping = {
'ait':'aitoff',
'mol':'mollweide',
'lam':'lambert',
'ham':'hammer'
}
kwargs.setdefault('proj','aitoff')
kwargs.setdefault('s',2)
kwargs.setdefault('marker','.')
kwargs.setdefault('c','k')
proj = kwargs.pop('proj')
projection = mapping.get(proj,proj)
#ax.grid()
# Convert from
# [0. < lon < 360] -> [-pi < lon < pi]
# [-90 < lat < 90] -> [-pi/2 < lat < pi/2]
lon,lat= np.radians([lon-360.*(lon>180),lat])
ax.scatter(lon,lat,**kwargs) Example 44
def healpixMap(nside, lon, lat, fill_value=0., nest=False):
"""
Input (lon, lat) in degrees instead of (theta, phi) in radians.
Returns HEALPix map at the desired resolution
"""
lon_median, lat_median = np.median(lon), np.median(lat)
max_angsep = np.max(ugali.utils.projector.angsep(lon, lat, lon_median, lat_median))
pix = angToPix(nside, lon, lat, nest=nest)
if max_angsep < 10:
# More efficient histograming for small regions of sky
m = np.tile(fill_value, healpy.nside2npix(nside))
pix_subset = ugali.utils.healpix.angToDisc(nside, lon_median, lat_median, max_angsep, nest=nest)
bins = np.arange(np.min(pix_subset), np.max(pix_subset) + 1)
m_subset = np.histogram(pix, bins=bins - 0.5)[0].astype(float)
m[bins[0:-1]] = m_subset
else:
m = np.histogram(pix, np.arange(hp.nside2npix(nside) + 1))[0].astype(float)
if fill_value != 0.:
m[m == 0.] = fill_value
return m
############################################################ Example 45
def galToCel(ll, bb):
"""
Converts Galactic (deg) to Celestial J2000 (deg) coordinates
"""
bb = numpy.radians(bb)
sin_bb = numpy.sin(bb)
cos_bb = numpy.cos(bb)
ll = numpy.radians(ll)
ra_gp = numpy.radians(192.85948)
de_gp = numpy.radians(27.12825)
lcp = numpy.radians(122.932)
sin_lcp_ll = numpy.sin(lcp - ll)
cos_lcp_ll = numpy.cos(lcp - ll)
sin_d = (numpy.sin(de_gp) * sin_bb) \
+ (numpy.cos(de_gp) * cos_bb * cos_lcp_ll)
ramragp = numpy.arctan2(cos_bb * sin_lcp_ll,
(numpy.cos(de_gp) * sin_bb) \
- (numpy.sin(de_gp) * cos_bb * cos_lcp_ll))
dec = numpy.arcsin(sin_d)
ra = (ramragp + ra_gp + (2. * numpy.pi)) % (2. * numpy.pi)
return numpy.degrees(ra), numpy.degrees(dec) Example 46
def ang2const(lon,lat,coord='gal'):
import ephem
scalar = np.isscalar(lon)
lon = np.array(lon,copy=False,ndmin=1)
lat = np.array(lat,copy=False,ndmin=1)
if coord.lower() == 'cel':
ra,dec = lon,lat
elif coord.lower() == 'gal':
ra,dec = gal2cel(lon,lat)
else:
msg = "Unrecognized coordinate"
raise Exception(msg)
x,y = np.radians([ra,dec])
const = [ephem.constellation(coord) for coord in zip(x,y)]
if scalar: return const[0]
return const Example 47
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 48
def __init__(self, survey_features=None, condition_features=None, time_lag=10.,
filtername='r', twi_change=-18.):
"""
Paramters
---------
time_lag : float (10.)
If there is a gap between observations longer than this, let the filter change (minutes)
twi_change : float (-18.)
The sun altitude to consider twilight starting/ending
"""
self.time_lag = time_lag/60./24. # Convert to days
self.twi_change = np.radians(twi_change)
self.filtername = filtername
if condition_features is None:
self.condition_features = {}
self.condition_features['Current_filter'] = features.Current_filter()
self.condition_features['Sun_moon_alts'] = features.Sun_moon_alts()
self.condition_features['Current_mjd'] = features.Current_mjd()
if survey_features is None:
self.survey_features = {}
self.survey_features['Last_observation'] = features.Last_observation()
super(Strict_filter_basis_function, self).__init__(survey_features=self.survey_features,
condition_features=self.condition_features) Example 49
def treexyz(ra, dec):
"""
Utility to convert RA,dec postions in x,y,z space, useful for constructing KD-trees.
Parameters
----------
ra : float or array
RA in radians
dec : float or array
Dec in radians
Returns
-------
x,y,z : floats or arrays
The position of the given points on the unit sphere.
"""
# Note ra/dec can be arrays.
x = np.cos(dec) * np.cos(ra)
y = np.cos(dec) * np.sin(ra)
z = np.sin(dec)
return x, y, z Example 50
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix((10., 10., 10.), (90., 90., 90.))
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array((
( a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0),
(-a*sinb*co, b*sina, 0.0, 0.0),
( a*cosb, b*cosa, c, 0.0),
( 0.0, 0.0, 0.0, 1.0)),
dtype=numpy.float64) 