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

'''

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
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
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "deg":
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
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "deg":
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
centerx = 0
centery = 0
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
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:
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
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

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())

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

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':
elif type == 'uoc':
elif type == 'clear_sky':
-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(

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)
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
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
"""
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.depth0 + (self.ndepth - 1) * self.ddepth)
self.theta0 = ?/2 - (self.lambda0 + (self.nlambda - 1) * self.dlambda)

Example 15

```def _add_triangular_sides(self, xy_mask, angle, y_top_right, y_bot_left,
x_top_right, x_bot_left, n_material):
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)

if i % y_per_iteration == 0:
lhs_x_start_index -= 1
rhs_x_stop_index += 1

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
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)
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
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
elif self.inc_prior == 'cos':
#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
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
elif self.inc_prior == 'cos':
#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
"""
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
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':
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

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
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(
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:
if tilt is not None:
if pitch is not None:
rvec = euler2rvec(r, t, p)

self._pose = tvec, rvec ```

Example 27

```def hav(alpha):
""" Formula for haversine

Parameters
----------
alpha : (float)

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)
"""
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:

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.
'''

colat_r = 90.0-lat_r
colat_0 = 90.0-lat_0

x_r = lon_r - lon_0
y_r = colat_0 - colat_r

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.
'''

colat_r = 90.0-lat_r
colat_0 = 90.0-lat_0

x_r = lon_r - lon_0
y_r = colat_0 - colat_r

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):

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

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:
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
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
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
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):
"""

parameters
----------
cluster: Cluster object
mpcscale: float
scaling to go from mpc to degrees (check units) at cluster redshift

results
-------

"""

# note this probably can be in the superclass, no?

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
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
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]
ax.scatter(lon,lat,**kwargs) ```

Example 44

```def healpixMap(nside, lon, lat, fill_value=0., nest=False):
"""
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
"""
sin_bb = numpy.sin(bb)
cos_bb = numpy.cos(bb)

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)

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
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.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
dec : float or array

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