def iterateRegio(self, pl, rwa, rwd, robl, rpoh, lon):
		
		okGa = okGd = True

		if pl.speculums[0][Planet.PMP] < 90.0 or (pl.speculums[0][Planet.PMP] >= 270.0 and pl.speculums[0][Planet.PMP] < 360.0):
			Ga = math.degrees(math.cos(rwa)*math.cos(robl)-math.sin(robl)*math.tan(rpoh))
			if Ga != 0.0:
				Fa = math.degrees(math.atan(math.sin(rwa)/(math.cos(rwa)*math.cos(robl)-math.sin(robl)*math.tan(rpoh))))

				if Fa >= 0.0 and Ga > 0.0:
					lon = Fa
				elif Fa < 0.0 and Ga > 0.0:
					lon = Fa+360.0
				elif Ga < 0.0:
					lon = Fa+180.0
			else:
				okGa = False
		else:
			Gd = math.degrees(math.cos(rwd)*math.cos(robl)+math.sin(robl)*math.tan(rpoh))
			if Gd != 0.0:
				Fd = math.degrees(math.atan(math.sin(rwd)/(math.cos(rwd)*math.cos(robl)+math.sin(robl)*math.tan(rpoh))))

				if Fd >= 0.0 and Gd > 0.0:
					lon = Fd
				elif Fd < 0.0 and Gd > 0.0:
					lon = Fd+360.0
				elif Gd < 0.0:
					lon = Fd+180.0
			else:
				okGd = False

		return okGa, okGd, lon

def ene2mots(energy, mat=None, hkl=None, r=None, alpha=None, pp=False):
    """calculates the real positions of the motors for a given energy (kev)

    Returns a list with positions: [atheh1, axeh1, dtheh1, dxeh1, dyeh1]
    """
    if mat is None:
        mat = CRYST_MAT
    if hkl is None:
        hkl = CRYST_HKL
    if r is None:
        r = CRYST_R
    if alpha is None:
        alpha = CRYST_ALPHA
    rthetab = theta_b(energy, get_dspacing(mat, hkl))
    ralpha = math.radians(alpha)
    p0 = r * math.sin(rthetab + ralpha)
    q0 = r * math.sin(rthetab - ralpha)

    atheh1 = math.degrees(rthetab)
    axeh1 = p0
    dtheh1 = 2 * math.degrees(rthetab)
    dxeh1 = p0 + q0 * math.cos(2 * rthetab)
    dyeh1 = q0 * math.sin(2 * rthetab)

    _mot_list = [atheh1, axeh1, dtheh1, dxeh1, dyeh1]

    if pp:
        # pretty print (= for humans)
        _tmpl_head = "MOT: {0:=^10} {1:=^10} {2:=^10} {3:=^10} {4:=^10}"
        _tmpl_data = "POS: {0:^ 10.4f} {1:^ 10.4f} {2:^ 10.4f} {3:^ 10.4f} {4:^ 10.4f}"
        print(_tmpl_head.format('ath', 'ax', 'dth', 'dx', 'dy'))
        print(_tmpl_data.format(*_mot_list))
    else:
        return _mot_list

def strokeOrientation(inStroke):
    "Input: List inPoints.  Returns the Angle of Orientation of a set of points (in degrees) where 0 is horizontal"

    cen = inStroke.Center
    inPoints = inStroke.Points
    pArea = area(inPoints)

    if (len(inPoints) <= 1):
        print "Warning: trying to get the Angle of Orientation of one or fewer points."
        return 0.0

    if pArea == 0:  #Perfect line        
        orientX = inPoints[-1].X - inPoints[0].X
        orientY = inPoints[-1].Y - inPoints[0].Y
        orientMag = vectorLength(orientX, orientY)
        if orientMag == 0:
            return 0.0
        orientX = orientX / orientMag
        orientY = orientY / orientMag
        return math.degrees( math.acos(orientX) )

    moment11 = momentOfOrder(cen, inPoints, 1, 1)
    if moment11 == 0:
        return 0    #There is no Moment of order 1,1.  We'd get a divide by zero.  Orientation is undefined; just return zero.

    angle = (.5 * math.atan((momentOfOrder(cen, inPoints, 0, 2) - momentOfOrder(cen, inPoints, 2, 0)) / (2 * moment11))) \
            + sign(moment11) * math.pi / 4

    return math.degrees(angle)

    def _calculate_position(self, radius, icon_size, index, children_count,
                            width, height, sin=math.sin, cos=math.cos):
        # tweak cos and sin in order to make the 'ring' into an equilateral
        # triangle.

        def cos_d(d):
            while d < -90:
                d += 360
            if d <= 30:
                return (d + 90) / 120.
            if d <= 90:
                return (90 - d) / 60.
            # mirror around 90
            return -cos_d(180 - d)

        sqrt_3 = math.sqrt(3)

        def sin_d(d):
            while d < -90:
                d += 360
            if d <= 30:
                return ((d + 90) / 120.) * sqrt_3 - 1
            if d <= 90:
                return sqrt_3 - 1
            # mirror around 90
            return sin_d(180 - d)

        cos = lambda r: cos_d(math.degrees(r))
        sin = lambda r: sin_d(math.degrees(r))

        return RingLayout._calculate_position(self, radius, icon_size, index,
                                              children_count, width, height,
                                              sin=sin, cos=cos)

    def getValues(self):        
        accx = self.twos_comp_combine(self.BUS.read_byte_data(self.LSM, self.ACC_X_MSB), self.BUS.read_byte_data(self.LSM, self.ACC_X_LSB))
        accy = self.twos_comp_combine(self.BUS.read_byte_data(self.LSM, self.ACC_Y_MSB), self.BUS.read_byte_data(self.LSM, self.ACC_Y_LSB))
        accz = self.twos_comp_combine(self.BUS.read_byte_data(self.LSM, self.ACC_Z_MSB), self.BUS.read_byte_data(self.LSM, self.ACC_Z_LSB))

        gyrox = self.twos_comp_combine(self.BUS.read_byte_data(self.GYRO, self.GYRO_X_MSB), self.BUS.read_byte_data(self.GYRO, self.GYRO_X_LSB))
        gyroy = self.twos_comp_combine(self.BUS.read_byte_data(self.GYRO, self.GYRO_Y_MSB), self.BUS.read_byte_data(self.GYRO, self.GYRO_Y_LSB))
        gyroz = self.twos_comp_combine(self.BUS.read_byte_data(self.GYRO, self.GYRO_Z_MSB), self.BUS.read_byte_data(self.GYRO, self.GYRO_Z_LSB))

        rate_gyrox = gyrox*self.GYRO_ADD
        rate_gyroy = gyroy*self.GYRO_ADD
        rate_gyroz = gyroz*self.GYRO_ADD

        if (not self.FIRST):
            self.FIRST = True
            self.gyroXangle = rate_gyrox*self.DT
            self.gyroYangle = rate_gyroy*self.DT
            self.gyroZangle = rate_gyroz*self.DT
        else:
            self.gyroXangle += rate_gyrox*self.DT
            self.gyroYangle += rate_gyroy*self.DT
            self.gyroZangle += rate_gyroz*self.DT

        roll = int(round(math.degrees(math.atan2(accx, accz))))
        pitch = int(round(math.degrees(math.atan2(accy, accz))))

        print "Przechylenie: ", int(round(roll,0)), " Pochylenie: ", int(round(pitch,0))

        self.FILTR_X = self.MULTIPLY*(roll)+(1-self.MULTIPLY)*self.gyroXangle
        self.FILTR_Y = self.MULTIPLY*(pitch)+(1-self.MULTIPLY)*self.gyroYangle

        print "Filtr przechylenie: ", int(round(self.FILTR_X,0)), " Filtr pochylenie: ", int(round(self.FILTR_Y,0))

        return str(roll)+';'+str(pitch)

def draw(x,y,count,isGrant,val,acceptCount):
	List_Of_Turtles = list()
	p = list()
	screen = turtle.getscreen()
	screen.setup( width = 2000, height = 2000, startx = None, starty = None) 
	for i in range(count):
	    screen.tracer(10)
	    List_Of_Turtles.append(turtle.Turtle())
	    #print "here",i,acceptCount,count
	    List_Of_Turtles[i].color("red")
	    if i >= acceptCount-1 and acceptCount!= 0:
		List_Of_Turtles[i].shape("square")
	    	List_Of_Turtles[i].color("green")
	    List_Of_Turtles[i].speed(1)
	    List_Of_Turtles[i].width(4)
	    angle = math.atan((y[i]-x[i])/300.0)
	    p.append(math.sqrt(90000+(x[i]-y[i])*(x[i]-y[i])))
	    if val == 200:
		Position_Set(List_Of_Turtles[i],val-275,200-x[i])
	    else:
		Position_Set(List_Of_Turtles[i],val,200-x[i])
	    if isGrant == 1:
		    List_Of_Turtles[i].right(180-math.degrees(angle)+90)
	    else:
		    List_Of_Turtles[i].right(math.degrees(angle)+90)
	    screen.update()
    
	for i in xrange(100):
		j=0
		for t in List_Of_Turtles:
			t.down()
			if (i*(i+1))/2 < p[j]:
				t.forward(i)
			j=j+1
		screen.update()

 def gazePix(self,anglex,angley):
     """Converts gaze angle to monitor pixel"""
     alphax = math.degrees(math.atan(self.ledx/self.monitordistance))
     pixx = self.deg2pix(anglex-alphax)
     alphay = math.degrees(math.atan(self.ledy/self.monitordistance))
     pixy = self.deg2pix(angley-alphay)
     return pixx,pixy

    def move_angle(self, ang, angvel=None, blocking=True):
        ''' move to angle (radians)
        '''
        if angvel == None:
            angvel = self.settings['max_speed']

        if angvel > self.settings['max_speed']:
            print 'lib_robotis.move_angle: angvel too high - %.2f deg/s' % (math.degrees(angvel))
            print 'lib_robotis.ignoring move command.'
            return

        if ang > self.settings['max_ang'] or ang < self.settings['min_ang']:
            print 'lib_robotis.move_angle: angle out of range- ', math.degrees(ang)
            print 'lib_robotis.ignoring move command.'
            return
        
        self.set_angvel(angvel)

        if self.settings['flipped']:
            ang = ang * -1.0
        enc_tics = int(round( ang / self.settings['rad_per_enc'] ))
        enc_tics += self.settings['home_encoder']
        self.move_to_encoder( enc_tics )

        if blocking == True:
            while(self.is_moving()):
                continue

def generate_wind(nodes, cells, prevailing_wind_bearing=270, average_wind_speed=7.0):
    import noise
    import math
    import numpy as np
    for n in nodes:
        wind_bearing = (int(round(180.0 * noise.pnoise3(n.x, n.y, n.z))) + 360) % 360  # in degrees
        wind_strength = (8 * noise.pnoise3(n.x, n.y, n.z)) + 16.0  # kmph

        # move the generated random wind bearing towards the prevailing wind a little
        wind_vector = np.add(np.array([wind_strength * math.cos(math.radians(90 - wind_bearing)),
                                       wind_strength * math.cos(math.radians(wind_bearing))]),
                             np.array([average_wind_speed * math.cos(math.radians(90 - prevailing_wind_bearing)),
                                       average_wind_speed * math.cos(math.radians(prevailing_wind_bearing))]))

        wind_strength = math.sqrt(math.pow(wind_vector[0], 2) + math.pow(wind_vector[1], 2)) #  sqrt((x2-x1)^2 + (y2-y1)^2) = vector magnitude
        wind_bearing = int(round(math.degrees(math.atan(wind_vector[1]/wind_vector[0])))) # tan-1((y2-y1)/x2-x1)) = vector bearing

        n.set_feature(['wind', 'bearing'], wind_bearing)
        n.set_feature(['wind', 'strength'], wind_strength)
        n.set_feature(['wind', 'vector', 'x'], wind_vector[0])
        n.set_feature(['wind', 'vector', 'y'], wind_vector[1])
        n.wind = Wind(wind_bearing, wind_strength, wind_vector)
    for c in cells:
        wind_vector = np.sum(np.array([n.wind.vector for n in c.nodes]), axis=0)
        wind_strength = math.sqrt(math.pow(wind_vector[0], 2) + math.pow(wind_vector[1], 2))  #  sqrt((x2-x1)^2 + (y2-y1)^2) = vector magnitude
        wind_bearing = int(round(math.degrees(math.atan(wind_vector[1]/wind_vector[0]))))  # tan-1((y2-y1)/x2-x1)) = vector bearing

        c.wind = Wind(wind_bearing, wind_strength, wind_vector)

	def compute(self):
		self.e.compute(datetime.datetime.utcnow())
		self.long  = math.degrees(float(self.e.sublong))
		self.lat   = math.degrees(float(self.e.sublat))
		self.x     = (self.long * 128/45) + 512
		self.y     = (self.lat * 128/45) + 256 + self.yoffset
		self.label = pyglet.text.Label(self.e.name, x=7,y=0, anchor_y="center", color=(255,255,255,255))

	def init_line(self):
		self.lines, self.vline_list, current_line = [], [], []
		for x in xrange(-total,total):
			temp = datetime.datetime.utcnow() + datetime.timedelta(seconds=interval*x)
			self.e.compute(temp)
			x = (math.degrees(float(self.e.sublong)) * 128/45) + 512
			y = (math.degrees(float(self.e.sublat)) * 128/45) + 256 + self.yoffset
			if len(current_line) > 1:
				# TO AVOID LINE FROM LEFT TO RIGHT
				temp_x, temp_y = current_line[-2], current_line[-1]
				if temp_x - x > 600:
					# From right edge to left edge
					current_line.extend((x+1024,y))
					self.lines.append(current_line)
					current_line = []
					current_line.extend((temp_x-1024,temp_y))
				elif temp_x - x < -600:
					# From left edge to right edge
					current_line.extend((x-1024,y))
					self.lines.append(current_line)
					current_line = []
					current_line.extend((temp_x+1024,temp_y))
			current_line.extend((x,y))
		self.lines.append(current_line)
		for x in self.lines:
			self.vline_list.append(pyglet.graphics.vertex_list(len(x)/2, ("v2f", x)))

    def connectionRound(self, first, last, pen, close):
        angle_1 = radians(degrees(self.prevAngle)+90)
        angle_2 = radians(degrees(self.currentAngle)+90)

        tempFirst = first - self.pointClass(cos(angle_1), sin(angle_1)) * self.miterLimit
        tempLast = last + self.pointClass(cos(angle_2), sin(angle_2)) * self.miterLimit

        newPoint = interSect((first, tempFirst), (last, tempLast))
        if newPoint is None:
            pen.lineTo(last)
            return
        #print "(%s, %s)," % (newPoint.x, newPoint.y)
        distance1 = newPoint.distance(first)
        distance2 = newPoint.distance(last)
        #print distance1, distance2
        if roundFloat(distance1) > self.miterLimit + self.contrast:
            distance1 = self.miterLimit + tempFirst.distance(tempLast) * .7
        if roundFloat(distance2) > self.miterLimit + self.contrast:
            distance2 = self.miterLimit + tempFirst.distance(tempLast) * .7

        distance1 *= self.magicCurve
        distance2 *= self.magicCurve

        bcp1 = first - self.pointClass(cos(angle_1), sin(angle_1)) * distance1
        bcp2 = last + self.pointClass(cos(angle_2), sin(angle_2)) * distance2
        pen.curveTo(bcp1, bcp2, last)

 def initialize(self):
     pos_std = self.parameter("pos_std", default=0.05)
     if isinstance(pos_std, dict):
         self._pos_std_dev = pos_std
     else:
         self._pos_std_dev = {'x': float(pos_std), 'y': float(pos_std), 'z': float(pos_std)}
     rot_std = self.parameter("rot_std", default=radians(5))
     if isinstance(rot_std, dict):
         self._rot_std_dev = rot_std
     else:
         self._rot_std_dev = {'roll': float(rot_std), 'pitch': float(rot_std), 'yaw': float(rot_std)}
     self._2D = bool(self.parameter("_2D", default=False))
     if self._2D:
         logger.info("Noise modifier standard deviations: x:%.4f, y:%.4f, yaw:%.3f deg",
                     self._pos_std_dev.get('x', 0),
                     self._pos_std_dev.get('y', 0),
                     degrees(self._rot_std_dev.get('yaw', 0)))
     else:
         logger.info("Noise modifier standard deviations: x:%.4f, y:%.4f, z:%.4f, "
                     "roll:%.3f deg, pitch:%.3f deg, yaw:%.3f deg",
                     self._pos_std_dev.get('x', 0),
                     self._pos_std_dev.get('y', 0),
                     self._pos_std_dev.get('z', 0),
                     degrees(self._rot_std_dev.get('roll', 0)),
                     degrees(self._rot_std_dev.get('pitch', 0)),
                     degrees(self._rot_std_dev.get('yaw', 0)))

def get_location(tle, now=None, lat=None, lng=None):
    """Compute the current location of the ISS"""
    now = now or datetime.datetime.utcnow()
    lat = lat or 37.7701
    lng = lng or -122.4664

    satellite = ephem.readtle(str(tle[0]), str(tle[1]), str(tle[2]))

    # Compute for current location
    observer = ephem.Observer()
    observer.lat = lat
    observer.lon = lng
    observer.elevation = 0
    observer.date = now
    satellite.compute(observer)
    lon = degrees(satellite.sublong)
    lat = degrees(satellite.sublat)

    # Return the relevant timestamp and data
    data = {'position': {'latitude': lat,
                         'longitude': lon},
            'visible': float(repr(satellite.alt)) > 0 and float(repr(satellite.alt)) < math.pi,
            'altitude': satellite.alt,
            'azimuth': satellite.az,
            'range': satellite.range,
            'velocity': satellite.range_velocity,
            'name': satellite.name}
    return data

def nearby(lat,lng,radius):
    radius = float(radius) / 1000
    lat = float(lat)
    lng = float(lng)

    maxLat = lat + degrees(radius/EARTH_RADIUS)
    minLat = lat - degrees(radius/EARTH_RADIUS)

    maxLng = lng + degrees(radius/EARTH_RADIUS/cos(radians(lat)))
    minLng = lng - degrees(radius/EARTH_RADIUS/cos(radians(lat)))

    lat = radians(lat)
    lng = radians(lng)


    rawSet = Spot.objects.raw('Select id, name, lat, lng, \
            acos(sin(%s)*sin(radians(lat)) + cos(%s)*cos(radians(lat))*cos(radians(lng)-%s)) * %s As D \
        From ( \
            Select id, name, lat, lng \
            From spots_spot \
            Where lat Between %s And %s \
              And lng Between %s And %s \
        ) As FirstCut \
        Where acos(sin(%s)*sin(radians(lat)) + cos(%s)*cos(radians(lat))*cos(radians(lng)-%s)) * %s < %s \
        Order by D', [lat,lat,lng,EARTH_RADIUS,minLat,maxLat,minLng,maxLng,lat,lat,lng,EARTH_RADIUS,radius])
    return rawSet

  def pitch_roll(self, px, pz):
    """works out the pitch (rx) and roll (rz) to apply to an object
    on the surface of the map at this point

    * returns a tuple (pitch, roll) in degrees

    Arguments:
      *px*
        x location
      *pz*
        z location
    """
    px -= self.unif[0]
    pz -= self.unif[2]
    halfw = self.width/2.0
    halfd = self.depth/2.0
    dx = self.width/self.ix
    dz = self.depth/self.iy
    x0 = int(math.floor((halfw + px)/dx + 0.5))
    if x0 < 0: x0 = 0
    if x0 > self.ix-1: x0 = self.ix-1
    z0 = int(math.floor((halfd + pz)/dz + 0.5))
    if z0 < 0: z0 = 0
    if z0 > self.iy-1: z0 = self.iy-1
    normp = array(self.buf[0].array_buffer[z0*self.ix + x0,3:6])
    # slight simplification to working out cross products as dirctn always 0,0,1
    #sidev = cross(normp, dirctn)
    sidev = array([normp[1], -normp[0], 0.0])
    sidev = sidev / sqrt(sidev.dot(sidev))
    #forwd = cross(sidev, normp)
    forwd = array([-normp[2]*normp[0], -normp[2]*normp[1],
                  normp[0]*normp[0] + normp[1]*normp[1]])
    forwd = forwd / sqrt(forwd.dot(forwd))
    return (degrees(asin(-forwd[1])), degrees(atan2(sidev[1], normp[1])))

def plot_planets(ax, plot, date, hill_sphere=False):
#     # only add the planets that would actually fall in this plot
#     plot_polygon = Polygon.Polygon(((plot[0],plot[2]),
#                                    (plot[0],plot[3]),
#                                    (plot[1],plot[3]),
#                                    (plot[0],plot[3]),
#                                    (plot[0],plot[2])))
# #    print plot_polygon
    mass = {"Sun":1.989*10**30, "Mars":639*10**21, "Jupiter":1.898*10**27, "Saturn":568.3*10**24, "Uranus":86.81*10**24, "Neptune":102.4*10**24}  # kg
    for planet in [ephem.Mars(), ephem.Jupiter(), ephem.Saturn(), ephem.Uranus(), ephem.Neptune()]:
        planet.compute(ephem.date(date))
        pos = (math.degrees(planet.ra), math.degrees(planet.dec))
#        if plot_polygon.isInside(math.degrees(planet.ra), math.degrees(planet.dec)):
        ax.scatter(pos[0], pos[1],
                 marker='o',
                 s=30,
                 facecolor='#E47833',
                 edgecolor='#E47833')
        ax.annotate(planet.name, (pos[0]+.4, pos[1]+0.15)) #(pos[0]+.9, pos[1]+0.5))  # offset to make it readable

        if hill_sphere:
            print planet.name, planet.sun_distance, mass[planet.name],
            hs_radius = (planet.sun_distance*ephem.meters_per_au)*((mass[planet.name]/3*mass['Sun'])**(1/3.))
            angular_size = planet.earth_distance*hs_radius  # FIXME
            print 'Hill sphere', hs_radius, hs_radius/ephem.meters_per_au, angular_size
            ax.add_patch(plt.Circle(pos, radius=angular_size, fill=False))

    return ax

def followCentre(data,desired_trajectory):
	global alpha

	a = getRange(data,130)
	b = getRange(data,179.9)
	swing = math.radians(50)
	print "center distances: ", a, b
	alpha = -math.atan((a*math.cos(swing)-b)/(a*math.sin(swing)))
	print "Alpha left",math.degrees(alpha)
	curr_dist1 = b*math.cos(alpha)
	future_dist1 = curr_dist1-car_length*math.sin(alpha)



	a = getRange(data,50)
	b = getRange(data,0)
	swing = math.radians(50)
	alpha = math.atan((a*math.cos(swing)-b)/(a*math.sin(swing)))
	print "Alpha right",math.degrees(alpha)
	curr_dist2 = b*math.cos(alpha)

	future_dist2 = curr_dist2+car_length*math.sin(alpha)

	desired_trajectory = (future_dist1 + future_dist2)/2

	print "dist 1 : ",future_dist1
	print "dist 2 : ",future_dist2
	# print "dist : ",future_dist
	error = future_dist1 - future_dist2
	print "Error : ",error
	return error, curr_dist2-curr_dist1

def toKepler(u, which = 'Pueyo', mass = 1, referenceTime = None):
    """
    """
    if which == 'Pueyo':
        res = np.zeros(6)
        res[1] = u[1]
        res[5] = u[5]
        
        res[0] = semimajoraxis(math.exp(u[0]), starMass = mass)
        res[2] = math.degrees(math.acos(u[2]))
        res[3] = np.mod((u[3]-u[4])*0.5,360)
        res[4] = np.mod((u[3]+u[4])*0.5,360)
        return res
    elif which == 'alternative':
        res = np.zeros(6)
        res[1] = u[1]
        res[5] = u[5]
        
        res[0] = semimajoraxis(math.exp(u[0]), starMass = mass)
        res[2] = math.degrees(math.acos(u[2]))
        res[3] = u[3]
        res[4] = u[4]
        return res        
    elif which == 'Chauvin':
        stat = StatisticsMCMC()
        res = stat.xFROMu(u,referenceTime,mass)    
        return res
    
    return None

def select_destination(origin='',
                       angle='',
                       radius='',
                       access_type='personal',
                       config_path='config/'):
    """
    Given a distance and polar angle, calculate the geocode of a destination point from the origin.
    """
    if origin == '':
        raise Exception('origin cannot be blank.')
    if angle == '':
        raise Exception('angle cannot be blank.')
    if radius == '':
        raise Exception('radius cannot be blank.')

    if isinstance(origin, str):
        origin_geocode = geocode_address(origin, access_type, config_path)
    elif isinstance(origin, list) and len(origin) == 2:
        origin_geocode = origin
    else:
        raise Exception('origin should be a list [lat, lng] or a string address.')

    # Find the location on a sphere a distance 'radius' along a bearing 'angle' from origin
    # This uses haversines rather than simple Pythagorean distance in Euclidean space
    #   because spheres are more complicated than planes.
    r = 3963.1676  # Radius of the Earth in miles
    bearing = radians(angle)  # Bearing in radians converted from angle in degrees
    lat1 = radians(origin_geocode[0])
    lng1 = radians(origin_geocode[1])
    lat2 = asin(sin(lat1) * cos(radius / r) + cos(lat1) * sin(radius / r) * cos(bearing))
    lng2 = lng1 + atan2(sin(bearing) * sin(radius / r) * cos(lat1), cos(radius / r) - sin(lat1) * sin(lat2))
    lat2 = degrees(lat2)
    lng2 = degrees(lng2)
    return [lat2, lng2]