我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.cos()。
def deriveKernel(self, params, i): self.checkParamsI(params, i) ell = np.exp(params[0]) p = np.exp(params[1]) #compute d2 if (self.K_sq is None): d2 = sq_dist(self.X_scaled.T / ell) #precompute squared distances else: d2 = self.K_sq / ell**2 #compute dp dp = self.dp/p K = np.exp(-d2 / 2.0) if (i==0): return d2*K*np.cos(2*np.pi*dp) elif (i==1): return 2*np.pi*dp*np.sin(2*np.pi*dp)*K else: raise Exception('invalid parameter index:' + str(i))
def getTrainTestKernel(self, params, Xtest): self.checkParams(params) ell = np.exp(params[0]) p = np.exp(params[1]) Xtest_scaled = Xtest/np.sqrt(Xtest.shape[1]) d2 = sq_dist(self.X_scaled.T/ell, Xtest_scaled.T/ell) #precompute squared distances #compute dp dp = np.zeros(d2.shape) for d in xrange(self.X_scaled.shape[1]): dp += (np.outer(self.X_scaled[:,d], np.ones((1, Xtest_scaled.shape[0]))) - np.outer(np.ones((self.X_scaled.shape[0], 1)), Xtest_scaled[:,d])) dp /= p K = np.exp(-d2 / 2.0) return np.cos(2*np.pi*dp)*K
def rotate_point_cloud(batch_data): """ Randomly rotate the point clouds to augument the dataset rotation is per shape based along up direction Input: BxNx3 array, original batch of point clouds Return: BxNx3 array, rotated batch of point clouds """ rotated_data = np.zeros(batch_data.shape, dtype=np.float32) for k in range(batch_data.shape[0]): rotation_angle = np.random.uniform() * 2 * np.pi cosval = np.cos(rotation_angle) sinval = np.sin(rotation_angle) rotation_matrix = np.array([[cosval, 0, sinval], [0, 1, 0], [-sinval, 0, cosval]]) shape_pc = batch_data[k, ...] rotated_data[k, ...] = np.dot(shape_pc.reshape((-1, 3)), rotation_matrix) return rotated_data
def rotate_point_cloud_by_angle(batch_data, rotation_angle): """ Rotate the point cloud along up direction with certain angle. Input: BxNx3 array, original batch of point clouds Return: BxNx3 array, rotated batch of point clouds """ rotated_data = np.zeros(batch_data.shape, dtype=np.float32) for k in range(batch_data.shape[0]): #rotation_angle = np.random.uniform() * 2 * np.pi cosval = np.cos(rotation_angle) sinval = np.sin(rotation_angle) rotation_matrix = np.array([[cosval, 0, sinval], [0, 1, 0], [-sinval, 0, cosval]]) shape_pc = batch_data[k, ...] rotated_data[k, ...] = np.dot(shape_pc.reshape((-1, 3)), rotation_matrix) return rotated_data
def _evalfull(self, x): fadd = self.fopt curshape, dim = self.shape_(x) # it is assumed x are row vectors if self.lastshape != curshape: self.initwithsize(curshape, dim) # BOUNDARY HANDLING # TRANSFORMATION IN SEARCH SPACE x = x - self.arrxopt x = monotoneTFosc(x) idx = (x > 0) x[idx] = x[idx] ** (1 + self.arrexpo[idx] * np.sqrt(x[idx])) x = self.arrscales * x # COMPUTATION core ftrue = 10 * (self.dim - np.sum(np.cos(2 * np.pi * x), -1)) + np.sum(x ** 2, -1) fval = self.noise(ftrue) # without noise # FINALIZE ftrue += fadd fval += fadd return fval, ftrue
def initwithsize(self, curshape, dim): # DIM-dependent initialization if self.dim != dim: if self.zerox: self.xopt = zeros(dim) else: self.xopt = compute_xopt(self.rseed, dim) self.rotation = compute_rotation(self.rseed + 1e6, dim) self.scales = (1. / self.condition ** .5) ** linspace(0, 1, dim) # CAVE? self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales)) # decouple scaling from function definition self.linearTF = dot(self.linearTF, self.rotation) K = np.arange(0, 12) self.aK = np.reshape(0.5 ** K, (1, 12)) self.bK = np.reshape(3. ** K, (1, 12)) self.f0 = np.sum(self.aK * np.cos(2 * np.pi * self.bK * 0.5)) # optimal value # DIM- and POPSI-dependent initialisations of DIM*POPSI matrices if self.lastshape != curshape: self.dim = dim self.lastshape = curshape self.arrxopt = resize(self.xopt, curshape)
def genplot(x, y, fit, xdata=None, ydata=None, maxpts=10000): bin_range = (0, 360) a = (np.arange(*bin_range)) f_a = nuth_func(a, fit[0], fit[1], fit[2]) nuth_func_str = r'$y=%0.2f*cos(%0.2f-x)+%0.2f$' % tuple(fit) if xdata.size > maxpts: import random idx = random.sample(list(range(xdata.size)), 10000) else: idx = np.arange(xdata.size) f, ax = plt.subplots() ax.set_xlabel('Aspect (deg)') ax.set_ylabel('dh/tan(slope) (m)') ax.plot(xdata[idx], ydata[idx], 'k.', label='Orig pixels') ax.plot(x, y, 'ro', label='Bin median') ax.axhline(color='k') ax.plot(a, f_a, 'b', label=nuth_func_str) ax.set_xlim(*bin_range) pad = 0.2 * np.max([np.abs(y.min()), np.abs(y.max())]) ax.set_ylim(y.min() - pad, y.max() + pad) ax.legend(prop={'size':8}) return f #Function copied from from openPIV pyprocess
def make_wafer(self,wafer_r,frame,label,labelloc,labelwidth): """ Generate wafer with primary flat on the left. From https://coresix.com/products/wafers/ I estimated that the angle defining the wafer flat to arctan(flat/2 / radius) """ angled = 18 angle = angled*np.pi/180 circ = cad.shapes.Circle((0,0), wafer_r, width=self.boxwidth, initial_angle=180+angled, final_angle=360+180-angled, layer=self.layer_box) flat = cad.core.Path([(-wafer_r*np.cos(angle),wafer_r*np.sin(angle)),(-wafer_r*np.cos(angle),-wafer_r*np.sin(angle))], width=self.boxwidth, layer=self.layer_box) date = time.strftime("%d/%m/%Y") if labelloc==(0,0): labelloc=(-2e3,wafer_r-1e3) # The label is added 100 um on top of the main cell label_grid_chip = cad.shapes.LineLabel( self.name + " " +\ date,500,position=labelloc, line_width=labelwidth, layer=self.layer_label) if frame==True: self.add(circ) self.add(flat) if label==True: self.add(label_grid_chip)
def ani_update(arg, ii=[0]): i = ii[0] # don't ask... if np.isclose(t_arr[i], np.around(t_arr[i], 1)): fig2.suptitle('Evolution (Time: {})'.format(t_arr[i]), fontsize=24) graphic_floor[0].set_data([-floor_lim*np.cos(incline_history[i]) + radius*np.sin(incline_history[i]), floor_lim*np.cos(incline_history[i]) + radius*np.sin(incline_history[i])], [-floor_lim*np.sin(incline_history[i])-radius*np.cos(incline_history[i]), floor_lim*np.sin(incline_history[i])-radius*np.cos(incline_history[i])]) graphic_wheel.center = (x_history[i], y_history[i]) graphic_ind[0].set_data([x_history[i], x_history[i] + radius*np.sin(w_history[i])], [y_history[i], y_history[i] + radius*np.cos(w_history[i])]) graphic_pend[0].set_data([x_history[i], x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2])], [y_history[i], y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2])]) graphic_dist[0].set_data([x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2]), x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2]) - pscale*p_history[i]*np.cos(q_history[i, 2])], [y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2]), y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2]) - pscale*p_history[i]*np.sin(q_history[i, 2])]) ii[0] += int(1 / (timestep * framerate)) if ii[0] >= len(t_arr): print("Resetting animation!") ii[0] = 0 return [graphic_floor, graphic_wheel, graphic_ind, graphic_pend, graphic_dist] # Run animation
def getZ2(el,recpos): """ sqrt( [a+h+s]^2 - [a*cos(el)]^2 ) - sqrt( [a+h]^2 - [a*cos(el)]^2 ) z(el) = --------------------------------------------------------------------- s a is height of observing station from earth center in km, h = 300km is height of ionosphere slab s = 200km is slab thickness """ a = np.linalg.norm(recpos)/1000 h=300 s = 200 Z = (np.sqrt((a+h+s)**2-(a*np.cos(np.radians(el[el>30])))**2) -np.sqrt((a+h)**2-(a*np.cos(np.radians(el[el>30])))**2))/s return Z
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 x_axis_rotation(theta): """Generates a 3x3 rotation matrix for a rotation of angle theta about the x axis. Parameters ---------- theta : float amount to rotate, in radians Returns ------- :obj:`numpy.ndarray` of float A random 3x3 rotation matrix. """ R = np.array([[1, 0, 0,], [0, np.cos(theta), -np.sin(theta)], [0, np.sin(theta), np.cos(theta)]]) return R
def y_axis_rotation(theta): """Generates a 3x3 rotation matrix for a rotation of angle theta about the y axis. Parameters ---------- theta : float amount to rotate, in radians Returns ------- :obj:`numpy.ndarray` of float A random 3x3 rotation matrix. """ R = np.array([[np.cos(theta), 0, np.sin(theta)], [0, 1, 0], [-np.sin(theta), 0, np.cos(theta)]]) return R
def z_axis_rotation(theta): """Generates a 3x3 rotation matrix for a rotation of angle theta about the z axis. Parameters ---------- theta : float amount to rotate, in radians Returns ------- :obj:`numpy.ndarray` of float A random 3x3 rotation matrix. """ R = np.array([[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]]) return R
def sph2cart(r, az, elev): """ Convert spherical to cartesian coordinates. Attributes ---------- r : float radius az : float aziumth (angle about z axis) elev : float elevation from xy plane Returns ------- float x-coordinate float y-coordinate float z-coordinate """ x = r * np.cos(az) * np.sin(elev) y = r * np.sin(az) * np.sin(elev) z = r * np.cos(elev) return x, y, z
def mdct(x, odd=True): """ Calculate modified discrete cosine transform of input signal in an inefficient pure-Python method. Use only for testing. Parameters ---------- X : array_like The input signal odd : boolean, optional Switch to oddly stacked transform. Defaults to :code:`True`. Returns ------- out : array_like The output signal """ return trans(x, func=numpy.cos, odd=odd) * numpy.sqrt(2)
def imdct(X, odd=True): """ Calculate inverse modified discrete cosine transform of input signal in an inefficient pure-Python method. Use only for testing. Parameters ---------- X : array_like The input signal odd : boolean, optional Switch to oddly stacked transform. Defaults to :code:`True`. Returns ------- out : array_like The output signal """ return itrans(X, func=numpy.cos, odd=odd) * numpy.sqrt(2)
def cmdct(x, odd=True): """ Calculate complex modified discrete cosine transform of input inefficient pure-Python method. Use only for testing. Parameters ---------- X : array_like The input signal odd : boolean, optional Switch to oddly stacked transform. Defaults to :code:`True`. Returns ------- out : array_like The output signal """ return trans(x, func=lambda x: numpy.cos(x) - 1j * numpy.sin(x), odd=odd)
def icmdct(X, odd=True): """ Calculate inverse complex modified discrete cosine transform of input signal in an inefficient pure-Python method. Use only for testing. Parameters ---------- X : array_like The input signal odd : boolean, optional Switch to oddly stacked transform. Defaults to :code:`True`. Returns ------- out : array_like The output signal """ return itrans(X, func=lambda x: numpy.cos(x) + 1j * numpy.sin(x), odd=odd)
def test_with_fake_log_prob(self): np.random.seed(42) def grad_log_prob(x): return -(x/2.0 + np.sin(x))*(1.0/2.0 + np.cos(x)) def fake_log_prob(x): return -(x/5.0 + np.sin(x) )**2.0/2.0 generator = mh_generator(log_density=fake_log_prob,x_start=1.0) tester = GaussianSteinTest(grad_log_prob,41) selector = SampleSelector(generator, sample_size=1000,thinning=20,tester=tester, max_iterations=5) data,converged = selector.points_from_stationary() assert converged is False
def test_with_ugly(self): np.random.seed(42) def grad_log_prob(x): return -(x/5.0 + np.sin(x))*(1.0/5.0 + np.cos(x)) def log_prob(x): return -(x/5.0 + np.sin(x) )**2.0/2.0 generator = mh_generator(log_density=log_prob,x_start=1.0) tester = GaussianSteinTest(grad_log_prob,41) selector = SampleSelector(generator, sample_size=1000,thinning=20,tester=tester, max_iterations=5) data,converged = selector.points_from_stationary() tester = GaussianSteinTest(grad_log_prob,41) assert tester.compute_pvalue(data)>0.05 assert converged is True
def random_points_in_circle(n,xx,yy,rr): """ get n random points in a circle. """ rnd = random(size=(n,3)) t = TWOPI*rnd[:,0] u = rnd[:,1:].sum(axis=1) r = zeros(n,'float') mask = u>1. xmask = logical_not(mask) r[mask] = 2.-u[mask] r[xmask] = u[xmask] xyp = reshape(rr*r,(n,1))*column_stack( (cos(t),sin(t)) ) dartsxy = xyp + array([xx,yy]) return dartsxy
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 genSphCoords(): """ Generates cartesian (x,y,z) and spherical (theta, phi) coordinates of a sphere Returns ------- coords : named tuple holds cartesian (x,y,z) and spherical (theta, phi) coordinates """ coords = namedtuple('coords', ['x', 'y', 'z', 'az', 'el']) az = _np.linspace(0, 2 * _np.pi, 360) el = _np.linspace(0, _np.pi, 181) coords.x = _np.outer(_np.cos(az), _np.sin(el)) coords.y = _np.outer(_np.sin(az), _np.sin(el)) coords.z = _np.outer(_np.ones(360), _np.cos(el)) coords.el, coords.az = _np.meshgrid(_np.linspace(0, _np.pi, 181), _np.linspace(0, 2 * _np.pi, 360)) return coords
def sph2cartMTX(vizMTX): """ Converts the spherical vizMTX data to named tuple contaibubg .xs/.ys/.zs Parameters ---------- vizMTX : array_like [180 x 360] matrix that hold amplitude information over phi and theta Returns ------- V : named_tuple Contains .xs, .ys, .zs cartesian coordinates """ rs = _np.abs(vizMTX.reshape((181, -1)).T) coords = genSphCoords() V = namedtuple('V', ['xs', 'ys', 'zs']) V.xs = rs * _np.sin(coords.el) * _np.cos(coords.az) V.ys = rs * _np.sin(coords.el) * _np.sin(coords.az) V.zs = rs * _np.cos(coords.el) return V
def convert_cof_mag2mass(t0, te, u0, alpha, s, q): """ function to convert from center of magnification to center of mass coordinates. Note that this function is for illustration only. It has not been tested and may have sign errors. """ if s <= 1.0: return t0, u0 else: delta = q / (1. + q) / s delta_u0 = delta * np.sin(alpha * np.pi / 180.) delta_tau = delta * np.cos(alpha * np.pi / 180.) t0_prime = t0 + delta_tau * te u0_prime = u0 + delta_u0 return t0_prime, u0_prime #Define model parameters in CoMAGN system
def _B_1_function(self, z, B_0=None): """ calculate B_1(z) function defined in: Gould A. 1994 ApJ 421L, 71 "Proper motions of MACHOs http://adsabs.harvard.edu/abs/1994ApJ...421L..71G Yoo J. et al. 2004 ApJ 603, 139 "OGLE-2003-BLG-262: Finite-Source Effects from a Point-Mass Lens" http://adsabs.harvard.edu/abs/2004ApJ...603..139Y """ if B_0 is None: B_0 = self._B_0_function(z) function = (lambda r, theta: r * np.sqrt(1.-r**2) / self.parameters.rho / np.sqrt(r**2+zz**2-2.*r*zz*np.cos(theta))) lim_0 = lambda x: 0 lim_1 = lambda x: 1 W_1 = 0. * z for (i, zz) in enumerate(z): W_1[i] = integrate.dblquad(function, 0., 2.*np.pi, lim_0, lim_1)[0] W_1 /= np.pi return B_0 - 1.5 * z * self.parameters.rho * W_1
def get_orientation_sector(dx,dy): # rotate (dx,dy) by pi/8 rotation = np.array([[np.cos(np.pi/8), -np.sin(np.pi/8)], [np.sin(np.pi/8), np.cos(np.pi/8)]]) rotated = np.dot(rotation, np.array([[dx], [dy]])) if rotated[1] < 0: rotated[0] = -rotated[0] rotated[1] = -rotated[1] s_theta = -1 if rotated[0] >= 0 and rotated[0] >= rotated[1]: s_theta = 0 elif rotated[0] >= 0 and rotated[0] < rotated[1]: s_theta = 1 elif rotated[0] < 0 and -rotated[0] < rotated[1]: s_theta = 2 elif rotated[0] < 0 and -rotated[0] >= rotated[1]: s_theta = 3 return s_theta
def to_radec(coords, xc=0, yc=0): """ Convert the generated coordinates to (ra, dec) (unit: degree). xc, yc: the center coordinate (ra, dec) """ results = [] for r, theta in coords: # FIXME: spherical algebra should be used!!! dx = r * np.cos(theta*np.pi/180) dy = r * np.sin(theta*np.pi/180) x = xc + dx y = yc + dy results.append((x, y)) if len(results) == 1: return results[0] else: return results
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 ct2lg(dX, dY, dZ, lat, lon): n = dX.size R = np.zeros((3, 3, n)) R[0, 0, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.cos(np.deg2rad(lon))) R[0, 1, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.sin(np.deg2rad(lon))) R[0, 2, :] = np.cos(np.deg2rad(lat)) R[1, 0, :] = -np.sin(np.deg2rad(lon)) R[1, 1, :] = np.cos(np.deg2rad(lon)) R[1, 2, :] = np.zeros((1, n)) R[2, 0, :] = np.multiply(np.cos(np.deg2rad(lat)), np.cos(np.deg2rad(lon))) R[2, 1, :] = np.multiply(np.cos(np.deg2rad(lat)), np.sin(np.deg2rad(lon))) R[2, 2, :] = np.sin(np.deg2rad(lat)) dxdydz = np.column_stack((np.column_stack((dX, dY)), dZ)) RR = np.reshape(R[0, :, :], (3, n)) dx = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0) RR = np.reshape(R[1, :, :], (3, n)) dy = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0) RR = np.reshape(R[2, :, :], (3, n)) dz = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0) return dx, dy, dz
def distance(self, lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1 = lon1*pi/180 lat1 = lat1*pi/180 lon2 = lon2*pi/180 lat2 = lat2*pi/180 # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2 c = 2 * np.arcsin(np.sqrt(a)) km = 6371 * c return km
def ct2lg(self, dX, dY, dZ, lat, lon): n = dX.size R = numpy.zeros((3, 3, n)) R[0, 0, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon))) R[0, 1, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon))) R[0, 2, :] = numpy.cos(numpy.deg2rad(lat)) R[1, 0, :] = -numpy.sin(numpy.deg2rad(lon)) R[1, 1, :] = numpy.cos(numpy.deg2rad(lon)) R[1, 2, :] = numpy.zeros((1, n)) R[2, 0, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon))) R[2, 1, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon))) R[2, 2, :] = numpy.sin(numpy.deg2rad(lat)) dxdydz = numpy.column_stack((numpy.column_stack((dX, dY)), dZ)) RR = numpy.reshape(R[0, :, :], (3, n)) dx = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0) RR = numpy.reshape(R[1, :, :], (3, n)) dy = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0) RR = numpy.reshape(R[2, :, :], (3, n)) dz = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0) return dx, dy, dz
def distance(self, lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1 = lon1*pi/180 lat1 = lat1*pi/180 lon2 = lon2*pi/180 lat2 = lat2*pi/180 # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2 c = 2 * numpy.arcsin(numpy.sqrt(a)) km = 6371 * c return km
def todictionary(self, time_series=False): # convert the ETM adjustment into a dirtionary # optionally, output the whole time series as well # start with the parameters etm = dict() etm['Linear'] = {'tref': self.Linear.tref, 'params': self.Linear.values.tolist()} etm['Jumps'] = [{'type':jump.type, 'year': jump.year, 'a': jump.a.tolist(), 'b': jump.b.tolist(), 'T': jump.T} for jump in self.Jumps.table] etm['Periodic'] = {'frequencies': self.Periodic.frequencies, 'sin': self.Periodic.sin.tolist(), 'cos': self.Periodic.cos.tolist()} if time_series: ts = dict() ts['t'] = self.ppp_soln.t.tolist() ts['x'] = self.ppp_soln.x.tolist() ts['y'] = self.ppp_soln.y.tolist() ts['z'] = self.ppp_soln.z.tolist() etm['time_series'] = ts return etm
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 plotPlaceTxRxSphereXY(Ax,xtx,ytx,xrx,yrx,x0,y0,a): Xlim = Ax.get_xlim() Ylim = Ax.get_ylim() FS = 20 Ax.scatter(xtx,ytx,s=100,color='k') Ax.text(xtx-0.75,ytx+1.5,'$\mathbf{Tx}$',fontsize=FS+6) Ax.scatter(xrx,yrx,s=100,color='k') Ax.text(xrx-0.75,yrx-4,'$\mathbf{Rx}$',fontsize=FS+6) xs = x0 + a*np.cos(np.linspace(0,2*np.pi,41)) ys = y0 + a*np.sin(np.linspace(0,2*np.pi,41)) Ax.plot(xs,ys,ls=':',color='k',linewidth=3) Ax.set_xbound(Xlim) Ax.set_ybound(Ylim) return Ax
def calc_IndCurrent_cos_range(self,f,t): """Induced current over a range of times""" Bpx = self.Bpx Bpz = self.Bpz a2 = self.a2 azm = np.pi*self.azm/180. R = self.R L = self.L w = 2*np.pi*f Ax = np.pi*a2**2*np.sin(azm) Az = np.pi*a2**2*np.cos(azm) Phi = (Ax*Bpx + Az*Bpz) phi = np.arctan(R/(w*L))-np.pi # This is the phase and not phase lag Is = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t + phi) Ire = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t)*np.cos(phi) Iim = (w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.sin(w*t)*np.sin(phi) return Ire,Iim,Is,phi
def calc_IndCurrent_FD_i(self,f): """Give FD EMF and current for single frequency""" #INITIALIZE ATTRIBUTES Bpx = self.Bpx Bpz = self.Bpz a2 = self.a2 azm = np.pi*self.azm/180. R = self.R L = self.L w = 2*np.pi*f Ax = np.pi*a2**2*np.sin(azm) Az = np.pi*a2**2*np.cos(azm) Phi = (Ax*Bpx + Az*Bpz) EMF = -1j*w*Phi Is = EMF/(R + 1j*w*L) return EMF,Is
def calc_IndCurrent_FD_spectrum(self): """Gives FD induced current spectrum""" #INITIALIZE ATTRIBUTES Bpx = self.Bpx Bpz = self.Bpz a2 = self.a2 azm = np.pi*self.azm/180. R = self.R L = self.L w = 2*np.pi*np.logspace(0,8,101) Ax = np.pi*a2**2*np.sin(azm) Az = np.pi*a2**2*np.cos(azm) Phi = (Ax*Bpx + Az*Bpz) EMF = -1j*w*Phi Is = EMF/(R + 1j*w*L) return EMF,Is
def calc_IndCurrent_TD_i(self,t): """Give FD EMF and current for single frequency""" #INITIALIZE ATTRIBUTES Bpx = self.Bpx Bpz = self.Bpz a2 = self.a2 azm = np.pi*self.azm/180. R = self.R L = self.L Ax = np.pi*a2**2*np.sin(azm) Az = np.pi*a2**2*np.cos(azm) Phi = (Ax*Bpx + Az*Bpz) Is = (Phi/L)*np.exp(-(R/L)*t) # V = (Phi*R/L)*np.exp(-(R/L)*t) - (Phi*R/L**2)*np.exp(-(R/L)*t) EMF = Phi return EMF,Is
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 save_hough(self, lines, clmap): """ :param lines: (rho, theta) pairs :param clmap: clusters assigned to lines :return: None """ height, width = self.image.shape ratio = 600. * (self.step+1) / min(height, width) temp = cv2.resize(self.image, None, fx=ratio, fy=ratio, interpolation=cv2.INTER_CUBIC) temp = cv2.cvtColor(temp, cv2.COLOR_GRAY2BGR) colors = [(0, 127, 255), (255, 0, 127)] for i in range(0, np.size(lines) / 2): rho = lines[i, 0] theta = lines[i, 1] color = colors[clmap[i, 0]] if theta < np.pi / 4 or theta > 3 * np.pi / 4: pt1 = (rho / np.cos(theta), 0) pt2 = (rho - height * np.sin(theta) / np.cos(theta), height) else: pt1 = (0, rho / np.sin(theta)) pt2 = (width, (rho - width * np.cos(theta)) / np.sin(theta)) pt1 = (int(pt1[0]), int(pt1[1])) pt2 = (int(pt2[0]), int(pt2[1])) cv2.line(temp, pt1, pt2, color, 5) self.save2image(temp)
def build_2D_cov_matrix(sigmax,sigmay,angle,verbose=True): """ Build a covariance matrix for a 2D multivariate Gaussian --- INPUT --- sigmax Standard deviation of the x-compoent of the multivariate Gaussian sigmay Standard deviation of the y-compoent of the multivariate Gaussian angle Angle to rotate matrix by in degrees (clockwise) to populate covariance cross terms verbose Toggle verbosity --- EXAMPLE OF USE --- import tdose_utilities as tu covmatrix = tu.build_2D_cov_matrix(3,1,35) """ if verbose: print ' - Build 2D covariance matrix with varinaces (x,y)=('+str(sigmax)+','+str(sigmay)+\ ') and then rotated '+str(angle)+' degrees' cov_orig = np.zeros([2,2]) cov_orig[0,0] = sigmay**2.0 cov_orig[1,1] = sigmax**2.0 angle_rad = (180.0-angle) * np.pi/180.0 # The (90-angle) makes sure the same convention as DS9 is used c, s = np.cos(angle_rad), np.sin(angle_rad) rotmatrix = np.matrix([[c, -s], [s, c]]) cov_rot = np.dot(np.dot(rotmatrix,cov_orig),np.transpose(rotmatrix)) # performing rot * cov * rot^T return cov_rot # = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
