def earth_distance(lat_lng1, lat_lng2):
"""
Compute the distance (in km) along earth between two latitude and longitude pairs
Parameters
----------
lat_lng1: tuple
the first latitude and longitude pair
lat_lng2: tuple
the second latitude and longitude pair
Returns
-------
float
the distance along earth in km
"""
lat1, lng1 = [l*pi/180 for l in lat_lng1]
lat2, lng2 = [l*pi/180 for l in lat_lng2]
dlat, dlng = lat1-lat2, lng1-lng2
ds = 2 * asin(sqrt(sin(dlat/2.0) ** 2 + cos(lat1) * cos(lat2) * sin(dlng/2.0) ** 2))
return 6371.01 * ds # spherical earth...
python类sin()的实例源码
def lonlat_to_pixel(self, lonlat, zoom):
"Converts a longitude, latitude coordinate pair for the given zoom level."
# Setting up, unpacking the longitude, latitude values and getting the
# number of pixels for the given zoom level.
lon, lat = self.get_lon_lat(lonlat)
npix = self._npix[zoom]
# Calculating the pixel x coordinate by multiplying the longitude value
# with the number of degrees/pixel at the given zoom level.
px_x = round(npix + (lon * self._degpp[zoom]))
# Creating the factor, and ensuring that 1 or -1 is not passed in as the
# base to the logarithm. Here's why:
# if fac = -1, we'll get log(0) which is undefined;
# if fac = 1, our logarithm base will be divided by 0, also undefined.
fac = min(max(sin(DTOR * lat), -0.9999), 0.9999)
# Calculating the pixel y coordinate.
px_y = round(npix + (0.5 * log((1 + fac) / (1 - fac)) * (-1.0 * self._radpp[zoom])))
# Returning the pixel x, y to the caller of the function.
return (px_x, px_y)
def sinox(x):
'''sin(x)/x'''
if(abs(x) > 1e-2):
return math.sin(x) / x
else:
return 1 - x*x / 6 + x*x*x*x / 120
def sinox3(x):
'''(x-sin(x))/(x*x*x)'''
if (abs(x) > 1e-2):
return (x - math.sin(x)) / (x*x*x)
else:
return 1./6. - x*x / 120. + x*x*x*x / 5040.
def specialFun1(x):
'''(x*sin(x) - 2.*(1.-cos(x)))/(x*x*x*x)'''
if (abs(x) > 1e-2):
return (x*math.sin(x) - 2.*(1. - math.cos(x))) / (x*x*x*x)
else:
return -1./12. + x*x / 180. - x*x*x*x / 6720.
def specialFun2(x):
'''(2.*(1.-cos(x)) - x*sin(x))/(2.*x*x*(1.-cos(x)))'''
if (abs(x) > 1e-2):
return (2.*(1. - math.cos(x)) - x*math.sin(x)) / (2.*x*x*(1. - math.cos(x)))
else:
return 1./12. + x*x / 720. + x*x*x*x / 30240.
def specialFun3(x):
'''(-2.*x + 3.*sin(x) - x*cos(x))/(x*x*x*x*x)'''
if (abs(x) > 1e-2):
return (-2.*x + 3.*math.sin(x) - x*math.cos(x)) / (x*x*x*x*x)
else:
return - 1./60. + x*x / 1260. - x*x*x*x / 60480.
def hexagon_generator(edge_length, offset):
"""Generator for coordinates in a hexagon."""
x, y = offset
for angle in range(0, 360, 60):
x += math.cos(math.radians(angle)) * edge_length
y += math.sin(math.radians(angle)) * edge_length
yield x, y
def hexagon_generator(self, edge_length, offset):
"""Generator for coordinates in a hexagon."""
x, y = offset
for angle in range(0, 360, 60):
x += math.cos(math.radians(angle)) * edge_length
y += math.sin(math.radians(angle)) * edge_length
yield x, y
def calc_helix_points(turtle, rad, pitch):
""" calculates required points to produce helix bezier curve with given radius and pitch in direction of turtle"""
# alpha = radians(90)
# pit = pitch/(2*pi)
# a_x = rad*cos(alpha)
# a_y = rad*sin(alpha)
# a = pit*alpha*(rad - a_x)*(3*rad - a_x)/(a_y*(4*rad - a_x)*tan(alpha))
# b_0 = Vector([a_x, -a_y, -alpha*pit])
# b_1 = Vector([(4*rad - a_x)/3, -(rad - a_x)*(3*rad - a_x)/(3*a_y), -a])
# b_2 = Vector([(4*rad - a_x)/3, (rad - a_x)*(3*rad - a_x)/(3*a_y), a])
# b_3 = Vector([a_x, a_y, alpha*pit])
# axis = Vector([0, 0, 1])
# simplifies greatly for case inc_angle = 90
points = [Vector([0, -rad, -pitch / 4]),
Vector([(4 * rad) / 3, -rad, 0]),
Vector([(4 * rad) / 3, rad, 0]),
Vector([0, rad, pitch / 4])]
# align helix points to turtle direction and randomize rotation around axis
trf = turtle.dir.to_track_quat('Z', 'Y')
spin_ang = rand_in_range(0, 2 * pi)
for p in points:
p.rotate(Quaternion(Vector([0, 0, 1]), spin_ang))
p.rotate(trf)
return points[1] - points[0], points[2] - points[0], points[3] - points[0], turtle.dir.copy()
def points_for_floor_split(self):
"""Calculate Poissonly distributed points for stem start points"""
array = []
# calculate approx spacing radius for dummy stem
self.tree_scale = self.param.g_scale + self.param.g_scale_v
stem = Stem(0, None)
stem.length = self.calc_stem_length(stem)
rad = 2.5 * self.calc_stem_radius(stem)
# generate points
for _ in range(self.param.floor_splits + 1):
point_ok = False
while not point_ok:
# distance from center proportional for number of splits, tree scale and stem radius
dis = sqrt(rand_in_range(0, 1) * self.param.floor_splits / 2.5 * self.param.g_scale * self.param.ratio)
# angle random in circle
theta = rand_in_range(0, 2 * pi)
pos = Vector([dis * cos(theta), dis * sin(theta), 0])
# test point against those already in array to ensure it will not intersect
point_m_ok = True
for point in array:
if (point[0] - pos).magnitude < rad:
point_m_ok = False
break
if point_m_ok:
point_ok = True
array.append((pos, theta))
return array
def shape_ratio(self, shape, ratio):
"""Calculate shape ratio as defined in paper"""
if shape == 1: # spherical
result = 0.2 + 0.8 * sin(pi * ratio)
elif shape == 2: # hemispherical
result = 0.2 + 0.8 * sin(0.5 * pi * ratio)
elif shape == 3: # cylindrical
result = 1.0
elif shape == 4: # tapered cylindrical
result = 0.5 + 0.5 * ratio
elif shape == 5: # flame
if ratio <= 0.7:
result = ratio / 0.7
else:
result = (1.0 - ratio) / 0.3
elif shape == 6: # inverse conical
result = 1.0 - 0.8 * ratio
elif shape == 7: # tend flame
if ratio <= 0.7:
result = 0.5 + 0.5 * ratio / 0.7
else:
result = 0.5 + 0.5 * (1.0 - ratio) / 0.3
elif shape == 8: # envelope
if ratio < 0 or ratio > 1:
result = 0.0
elif ratio < 1 - self.param.prune_width_peak:
result = pow(ratio / (1 - self.param.prune_width_peak),
self.param.prune_power_high)
else:
result = pow((1 - ratio) / (1 - self.param.prune_width_peak),
self.param.prune_power_low)
else: # conical (0)
result = 0.2 + 0.8 * ratio
return result
def q_prod(sym):
"""Production rule for Q"""
prop_off = sym.parameters["t"] / __t_max__
if prop_off < 1:
res = [LSymbol("!", {"w": 0.85 + 0.15 * sin(sym.parameters["t"])}),
LSymbol("^", {"a": random() - 0.65})]
if prop_off > __p_max__:
d_ang = 1 / (1 - __p_max__) * (1 - prop_off) * 110 + 15
res.extend([LSymbol("!", {"w": 0.1})])
for ind in range(int(random() * 2 + 5)):
r_ang = sym.parameters["t"] * 10 + ind * (random() * 50 + 40)
e_d_ang = d_ang * (random() * 0.4 + 0.8)
res.extend([LSymbol("/", {"a": r_ang}),
LSymbol("&", {"a": e_d_ang}),
LSymbol("["),
LSymbol("A"),
LSymbol("]"),
LSymbol("^", {"a": e_d_ang}),
LSymbol("\\", {"a": r_ang})],)
res.append(LSymbol("F", {"l": 0.05}))
else:
res.append(LSymbol("F", {"l": 0.15}))
res.append(LSymbol("Q", {"t": sym.parameters["t"] + __d_t__}))
else:
res = [LSymbol("!", {"w": 0}),
LSymbol("F", {"l": 0.15})]
return res
def calculate_cut_off(self):
"""
Calculate cut-off radius as Rc = L/2 from a given MOF object.
"""
width_a = self.ucv / (self.uc_size[1] * self.uc_size[2] / math.sin(math.radians(self.uc_angle[0])))
width_b = self.ucv / (self.uc_size[0] * self.uc_size[2] / math.sin(math.radians(self.uc_angle[1])))
width_c = self.ucv / (self.uc_size[0] * self.uc_size[1] / math.sin(math.radians(self.uc_angle[2])))
self.cut_off = min(width_a / 2, width_b / 2, width_c / 2)
def pbc_parameters(self):
"""
Calculates constants used in periodic boundary conditions.
"""
uc_cos = [math.cos(math.radians(a)) for a in self.uc_angle]
uc_sin = [math.sin(math.radians(a)) for a in self.uc_angle]
a, b, c = self.uc_size
v = self.frac_ucv
xf1 = 1 / a
xf2 = - uc_cos[2] / (a * uc_sin[2])
xf3 = (uc_cos[0] * uc_cos[2] - uc_cos[1]) / (a * v * uc_sin[2])
yf1 = 1 / (b * uc_sin[2])
yf2 = (uc_cos[1] * uc_cos[2] - uc_cos[0]) / (b * v * uc_sin[2])
zf1 = uc_sin[2] / (c * v)
self.to_frac = [xf1, xf2, xf3, yf1, yf2, zf1]
xc1 = a
xc2 = b * uc_cos[2]
xc3 = c * uc_cos[1]
yc1 = b * uc_sin[2]
yc2 = c * (uc_cos[0] - uc_cos[1] * uc_cos[2]) / uc_sin[2]
zc1 = c * v / uc_sin[2]
self.to_car = [xc1, xc2, xc3, yc1, yc2, zc1]
def uc_vectors(cls, uc_size, uc_angle):
"""
Calculate unit cell vectors for given unit cell size and angles
"""
a = uc_size[0]
b = uc_size[1]
c = uc_size[2]
alpha = math.radians(uc_angle[0])
beta = math.radians(uc_angle[1])
gamma = math.radians(uc_angle[2])
x_v = [a, 0, 0]
y_v = [b * math.cos(gamma), b * math.sin(gamma), 0]
z_v = [0.0] * 3
z_v[0] = c * math.cos(beta)
z_v[1] = (c * b * math.cos(alpha) - y_v[0] * z_v[0]) / y_v[1]
z_v[2] = math.sqrt(c * c - z_v[0] * z_v[0] - z_v[1] * z_v[1])
uc_vectors = [x_v, y_v, z_v]
return uc_vectors
def earth_Rreal(latrad):
return (1.0 / (((cos(latrad)) / earth_Rmax) ** 2 + ((sin(latrad)) / earth_Rmin) ** 2)) ** 0.5
def get_distance(location1, location2):
lat1, lng1 = location1
lat2, lng2 = location2
lat1, lng1, lat2, lng2 = map(radians, (lat1, lng1, lat2, lng2))
d = sin(0.5*(lat2 - lat1)) ** 2 + cos(lat1) * cos(lat2) * sin(0.5*(lng2 - lng1)) ** 2
return 2 * earth_Rrect * asin(sqrt(d))
def init_grid(self):
grid_all = []
lats = self.init_lats()
c = 2 * pi / (3 ** 0.5 * self.r_sight * self.safety) * self.earth_R
even_lng = True
strip_amount = int(ceil(c))
grid_all.append((0, strip_amount, even_lng))
ind_lat = 2
while ind_lat < len(lats):
amount = int(ceil(c * cos(lats[ind_lat])))
if amount < strip_amount - (sin(lats[ind_lat]*2)*self.param_shift+self.param_stretch):
ind_lat -= 1
strip_amount = int(ceil(c * cos(lats[ind_lat])))
else:
even_lng = not even_lng
if ind_lat + 1 < len(lats):
lat = lats[ind_lat + 1] * 180 / pi
grid_all.append((lat, strip_amount, even_lng))
ind_lat += 3
grid_all.append((90.0, 1, True)) # pole
return grid_all
def dist_cmp(self, location1, location2):
return sin(0.5 * (location2[0] - location1[0])) ** 2 + cos(location2[0]) * cos(location1[0]) * sin(0.5 * (location2[1] - location1[1])) ** 2
def getEarthRadius(latrad):
return (1.0 / (((math.cos(latrad)) / EARTH_Rmax) ** 2 + ((math.sin(latrad)) / EARTH_Rmin) ** 2)) ** (1.0 / 2)
def getCircularBounds(fitCloud=None,width=64,height=64,smoothing=0.01):
circumference = 2*(width+height)
if not fitCloud is None:
cx = np.mean(fitCloud[:,0])
cy = np.mean(fitCloud[:,1])
r = 0.5* max( np.max(fitCloud[:,0])- np.min(fitCloud[:,0]),np.max(fitCloud[:,1])- np.min(fitCloud[:,1]))
else:
r = circumference /(2.0*math.pi)
cx = cy = r
perimeterPoints = np.zeros((circumference,2),dtype=float)
for i in range(circumference):
angle = (2.0*math.pi)*float(i) / circumference - math.pi * 0.5
perimeterPoints[i][0] = cx + r * math.cos(angle)
perimeterPoints[i][1] = cy + r * math.sin(angle)
bounds = {'top':perimeterPoints[0:width],
'right':perimeterPoints[width-1:width+height-1],
'bottom':perimeterPoints[width+height-2:2*width+height-2],
'left':perimeterPoints[2*width+height-3:]}
bounds['s_top'],u = interpolate.splprep([bounds['top'][:,0], bounds['top'][:,1]],s=smoothing)
bounds['s_right'],u = interpolate.splprep([bounds['right'][:,0],bounds['right'][:,1]],s=smoothing)
bounds['s_bottom'],u = interpolate.splprep([bounds['bottom'][:,0],bounds['bottom'][:,1]],s=smoothing)
bounds['s_left'],u = interpolate.splprep([bounds['left'][:,0],bounds['left'][:,1]],s=smoothing)
return bounds
def geometry(self, frame=0):
t = frame / self.frames
Rot = Matrix.Rotation(0.5*pi, 4, 'Y')
bm = bmesh.new()
for i in range(self.n):
t0 = i / self.n
r0, theta = t0*self.r0, i*goldenAngle - frame*goldenAngle + t*self.offset
x = r0*cos(theta)
y = r0*sin(theta)
z = self.h0/2 - (self.h0 / (self.r0*self.r0))*r0*r0
p0 = Vector((x, y, z))
T0, N0, B0 = getTNBfromVector(p0)
M0 = Matrix([T0, B0, N0]).to_4x4().transposed()
for j in range(self.m):
t1 = j / self.m
t2 = 0.4 + 0.6*t0
r1, theta = t2*t1*self.r1, j*goldenAngle #- frame*goldenAngle + t*self.offset
x = r1*cos(theta)
y = r1*sin(theta)
z = self.h1 - (self.h1 / (self.r1*self.r1))*r1*r1
p1 = Vector((x, y, z))
T1, N1, B1 = getTNBfromVector(p1)
M1 = Matrix([T1, B1, N1]).to_4x4().transposed()
p = p0 + M0*p1
r2 = t2*t1*self.r2
T = Matrix.Translation(p)
bmesh.ops.create_cone(bm,
cap_ends=True, segments=6,
diameter1=r2, diameter2=r2,
depth=0.1*r2, matrix=T*M0*M1*Rot)
return bm
def torusSurface(r0, r1):
def surface(u, v):
point = ((r0 + r1*cos(TAU*v))*cos(TAU*u), \
(r0 + r1*cos(TAU*v))*sin(TAU*u), \
r1*sin(TAU*v))
return point
return surface
# Create an object from a surface parameterization
def rainbowLights(r=5, n=100, freq=2, energy=0.1):
for i in range(n):
t = float(i)/float(n)
pos = (r*sin(tau*t), r*cos(tau*t), r*sin(freq*tau*t))
# Create lamp
bpy.ops.object.add(type='LAMP', location=pos)
obj = bpy.context.object
obj.data.type = 'POINT'
# Apply gamma correction for Blender
color = tuple(pow(c, 2.2) for c in colorsys.hsv_to_rgb(t, 0.6, 1))
# Set HSV color and lamp energy
obj.data.color = color
obj.data.energy = energy
def radial(cls, r, theta):
"""Provide a radial acceleration.
Arguments:
r (float): speed in pixels per second (per second)
theta (float): angle in degrees (0 = +X axis, 90 = +Y axis)
"""
radians = math.radians(theta)
ax = r * math.cos(radians)
ay = r * math.sin(radians)
return cls(ax=ax, ay=ay)
def create_matrix(self,parent_joint,parent_joint_matrix):
# The calculation of the local matrix is an optimized version of
# local_matrix = T*IPS*R*S if ignore_parent_scale else T*R*S
# where S, R and T is the scale, rotation and translation matrix
# respectively and IPS is the inverse parent scale matrix.
cx = cos(radians(self.rotation_x))
sx = sin(radians(self.rotation_x))
cy = cos(radians(self.rotation_y))
sy = sin(radians(self.rotation_y))
cz = cos(radians(self.rotation_z))
sz = sin(radians(self.rotation_z))
if self.ignore_parent_scale:
ips_x = 1/parent_joint.scale_x
ips_y = 1/parent_joint.scale_y
ips_z = 1/parent_joint.scale_z
else:
ips_x = 1
ips_y = 1
ips_z = 1
local_matrix = numpy.empty((3,4),numpy.float32)
local_matrix[0,0] = cy*cz*self.scale_x*ips_x
local_matrix[1,0] = cy*sz*self.scale_x*ips_y
local_matrix[2,0] = -sy*self.scale_x*ips_z
local_matrix[0,1] = (sx*sy*cz - cx*sz)*self.scale_y*ips_x
local_matrix[1,1] = (sx*sy*sz + cx*cz)*self.scale_y*ips_y
local_matrix[2,1] = sx*cy*self.scale_y*ips_z
local_matrix[0,2] = (cx*sy*cz + sx*sz)*self.scale_z*ips_x
local_matrix[1,2] = (cx*sy*sz - sx*cz)*self.scale_z*ips_y
local_matrix[2,2] = cx*cy*self.scale_z*ips_z
local_matrix[0,3] = self.translation_x
local_matrix[1,3] = self.translation_y
local_matrix[2,3] = self.translation_z
return matrix3x4_multiply(parent_joint_matrix,local_matrix)
def create_matrix(self):
c = cos(radians(self.rotation))
s = sin(radians(self.rotation))
R = numpy.matrix([[c,-s,0],[s,c,0],[0,0,1]])
S = numpy.matrix([[self.scale_s,0,0],[0,self.scale_t,0],[0,0,1]])
C = numpy.matrix([[1,0,self.center_s],[0,1,self.center_t],[0,0,1]])
T = numpy.matrix([[1,0,self.translation_s],[0,1,self.translation_t],[0,0,1]])
# Only types 0x00, 0x06, 0x07, 0x08 and 0x09 have been tested
if self.matrix_type in {0x00,0x02,0x0A,0x0B,0x80}:
P = numpy.matrix([[1,0,0,0],[0,1,0,0],[0,0,0,1]])
elif self.matrix_type == 0x06:
P = numpy.matrix([[0.5,0,0,0.5],[0,-0.5,0,0.5],[0,0,0,1]])
elif self.matrix_type == 0x07:
P = numpy.matrix([[0.5,0,0.5,0],[0,-0.5,0.5,0],[0,0,1,0]])
elif self.matrix_type in {0x08,0x09}:
P = numpy.matrix([[0.5,0,0.5,0],[0,-0.5,0.5,0],[0,0,1,0]])*numpy.matrix(self.projection_matrix)
else:
raise ValueError('invalid texture matrix type')
M = T*C*S*R*C.I*P
if self.shape == gx.TG_MTX2x4:
return M[:2,:]
elif self.shape == gx.TG_MTX3x4:
return M
else:
raise ValueError('invalid texture matrix shape')
def update(self,time):
scale_x = self.scale_x.interpolate(time)
scale_y = self.scale_y.interpolate(time)
scale_z = self.scale_z.interpolate(time)
rotation_x = self.rotation_x.interpolate(time)
rotation_y = self.rotation_y.interpolate(time)
rotation_z = self.rotation_z.interpolate(time)
translation_x = self.translation_x.interpolate(time)
translation_y = self.translation_y.interpolate(time)
translation_z = self.translation_z.interpolate(time)
cx = cos(radians(rotation_x))
sx = sin(radians(rotation_x))
cy = cos(radians(rotation_y))
sy = sin(radians(rotation_y))
cz = cos(radians(rotation_z))
sz = sin(radians(rotation_z))
R = numpy.matrix([[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,1.0]]) #<-?
R[0,0] = cy*cz
R[0,1] = (sx*sy*cz - cx*sz)
R[0,2] = (cx*sy*cz + sx*sz)
R[1,0] = cy*sz
R[1,1] = (sx*sy*sz + cx*cz)
R[1,2] = (cx*sy*sz - sx*cz)
R[2,0] = -sy
R[2,1] = sx*cy
R[2,2] = cx*cy
S = numpy.matrix([[scale_x,0,0,0],[0,scale_y,0,0],[0,0,scale_z,0],[0,0,0,1]])
C = numpy.matrix([[1,0,0,self.center_x],[0,1,0,self.center_y],[0,0,1,self.center_z],[0,0,0,1]])
T = numpy.matrix([[1,0,0,translation_x],[0,1,0,translation_y],[0,0,1,translation_z],[0,0,0,1]])
self.texture_matrix[:] = (T*C*S*R*C.I)[:self.row_count,:]
def rotation(axis_x,axis_y,axis_z,angle):
s = sin(angle/2)
return Quarternion(cos(angle/2),s*axis_x,s*axis_y,s*axis_z)
