我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.radians()。
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
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 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 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 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 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
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 #########################################################################
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
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))
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) ])
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)
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)
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)
def find_lines(img, acc_threshold=0.25, should_erode=True): if len(img.shape) == 3 and img.shape[2] == 3: # if it's color img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) img = cv2.GaussianBlur(img, (11, 11), 0) img = cv2.adaptiveThreshold( img, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 5, 2) img = cv2.bitwise_not(img) # thresh = 127 # edges = cv2.threshold(img, thresh, 255, cv2.THRESH_BINARY)[1] # edges = cv2.Canny(blur, 500, 500, apertureSize=3) if should_erode: element = cv2.getStructuringElement(cv2.MORPH_RECT, (4, 4)) img = cv2.erode(img, element) theta = np.pi/2000 angle_threshold = 2 horizontal = cv2.HoughLines( img, 1, theta, int(acc_threshold * img.shape[1]), min_theta=np.radians(90 - angle_threshold), max_theta=np.radians(90 + angle_threshold)) vertical = cv2.HoughLines( img, 1, theta, int(acc_threshold * img.shape[0]), min_theta=np.radians(-angle_threshold), max_theta=np.radians(angle_threshold), ) horizontal = list(horizontal) if horizontal is not None else [] vertical = list(vertical) if vertical is not None else [] horizontal = [line[0] for line in horizontal] vertical = [line[0] for line in vertical] horizontal = np.asarray(horizontal) vertical = np.asarray(vertical) return horizontal, vertical
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)
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 ############################################################
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)
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
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)
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
