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
python类radians()的实例源码
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
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
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
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))
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]])
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))
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]])
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())
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'
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))
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' ]]
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)
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)
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
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]])
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)
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
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
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))
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)]])
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
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
########################################
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
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
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
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
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
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
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
