def gps_to_kmeters(long1, lat1, long2, lat2):
        """Converts gps longitudes and latitudes distance to meters

        Keyword arguments:
        long1 -- longitude for pos 1
        lat1 -- latitute for pos 1
        long2 -- longitude for pos 2
        lat2 -- latitude for pos 2
        """

        #approx radius of earth
        R = 6373.0

        long1 = radians(long1)
        lat1 = radians(lat1)
        long2 = radians(long2)
        lat2 = radians(lat2)

        dist_long = long2 - long1
        dist_lat = lat2 - lat1

        a = sin(dist_lat / 2)**2 + cos(lat1) * cos(lat2) * sin(dist_long)**2
        c = atan2(sqrt(a), sqrt(1 - a))

        dist = R * c

        return dist

def calc_dist(plyr_angl,plyr_vel):
    gravity = 9.81                #gravity 9.81 m/s^2
    vel_init_x = plyr_vel * math.cos(math.radians(plyr_angl))    #creates Vix
    vel_init_y = plyr_vel * math.sin(math.radians(plyr_angl))    #creates Viy
    time_in_air = 2 * (vel_init_y / gravity)              #solves for time
    distance_x = vel_init_x * time_in_air                #solves for distance
    return (distance_x, time_in_air)                    #returns horizontal distance, time in air

    def __create_board(self):
        """
        Creates a Canvas widget where Hexagon fields are marked on it
        """
        m = self.size[0]
        n = self.size[1]
        edge_length = 24
        # calculates the size of the Hexagon
        y_top, x_right = self.__Hex_size(edge_length)
        canvas_width = x_right * n + x_right * m / 2 + 50
        canvas_height = y_top * m + 100

        self.w = Canvas(self.master, width=canvas_width, height=canvas_height)
        self.w.configure(background=self.bg)
        self.w.pack()
        # creates Hexagon Grid
        for j in range(m):
            for i in range(n):
                x = 40 + x_right * i + x_right / 2 * j
                y = 50 + y_top * j
                k = 0
                for angle in range(0, 360, 60):
                    y += math.cos(math.radians(angle)) * edge_length
                    x += math.sin(math.radians(angle)) * edge_length
                    self.point_coordinates[j][i][k] = x
                    self.point_coordinates[j][i][k + 1] = y
                    k += 2
                # draws Hexagon to the canvas widget
                self.w.create_polygon(list(
                    self.point_coordinates[j][i]),
                    outline=self.tile_outline, fill=self.tile, width=3)

    def sumVectors(self, vectors):
        """ sum all vectors (including targetvector)"""
        endObstacleVector = (0,0)

        ##generate endvector of obstacles
        #sum obstaclevectors
        for vector in vectors:
            vectorX = math.sin(math.radians(vector[1])) * vector[0] # x-position
            vectorY = math.cos(math.radians(vector[1])) * vector[0] # y-position
            endObstacleVector = (endObstacleVector[0]+vectorX,endObstacleVector[1]+vectorY)
        #mean obstaclevectors
        if len(vectors) > 0:
            endObstacleVector = (endObstacleVector[0]/len(vectors), endObstacleVector[1]/len(vectors))

        #add targetvector
        targetVector = self.target
        if targetVector != 0 and targetVector != None:
            vectorX = math.sin(math.radians(targetVector[1])) * targetVector[0] # x-position
            vectorY = math.cos(math.radians(targetVector[1])) * targetVector[0] # y-position
            endVector = (endObstacleVector[0]+vectorX,endObstacleVector[1]+vectorY)
            #endVector = (endVector[0]/2, endVector[1]/2)
        else:
            endVector = endObstacleVector


        return endVector

    def update(self):
        if self.turning_right:
            self.rot -= 5
        if self.turning_left:
            self.rot += 5

        a = [0.0,0.0]
        if self.boost_endtime > rabbyt.get_time():
            f = 3*(self.boost_endtime - rabbyt.get_time())/self.boost_length
            a[0] += cos(radians(self.boost_rot))*f
            a[1] += sin(radians(self.boost_rot))*f
            self.create_boost_particle()

        if self.accelerating:
            a[0] += cos(radians(self.rot))*.9
            a[1] += sin(radians(self.rot))*.9
            self.create_dust_particle(self.dust_r)
            self.create_dust_particle(self.dust_l)

        ff = .9 # Friction Factor

        self.velocity[0] *= ff
        self.velocity[1] *= ff

        self.velocity[0] += a[0]
        self.velocity[1] += a[1]

        self.x += self.velocity[0]
        self.y += self.velocity[1]

    def set_3d(self):
        """ Configure OpenGL to draw in 3d.

        """
        width, height = self.get_size()

        gl.glEnable(gl.GL_DEPTH_TEST)

        gl.glViewport(0, 0, width, height)
        gl.glMatrixMode(gl.GL_PROJECTION)
        gl.glLoadIdentity()
        gl.gluPerspective(65.0, width / float(height), 0.1, DIST)
        gl.glMatrixMode(gl.GL_MODELVIEW)
        gl.glLoadIdentity()

        x, y = self.rotation
        gl.glRotatef(x, 0, 1, 0)
        gl.glRotatef(-y, math.cos(math.radians(x)), 0, math.sin(math.radians(x)))
        x, y, z = self.position
        gl.glTranslatef(-x, -y, -z)

        gl.glEnable(gl.GL_LIGHTING)
        gl.glLightModelfv(gl.GL_LIGHT_MODEL_AMBIENT, GLfloat4(0.05,0.05,0.05,1.0))
        gl.glEnable(gl.GL_COLOR_MATERIAL)
        gl.glColorMaterial(gl.GL_FRONT, gl.GL_AMBIENT_AND_DIFFUSE)
        #gl.glLightfv(gl.GL_LIGHT1,gl.GL_SPOT_DIRECTION, GLfloat3(0,0,-1))
        gl.glLightfv(gl.GL_LIGHT1, gl.GL_AMBIENT, GLfloat4(0.5,0.5,0.5,1.0))
        gl.glLightfv(gl.GL_LIGHT1, gl.GL_DIFFUSE, GLfloat4(1.0,1.0,1.0,1.0))
        gl.glLightfv(gl.GL_LIGHT1, gl.GL_POSITION, GLfloat4(0.35,1.0,0.65,0.0))
        #gl.glLightfv(gl.GL_LIGHT0,gl.GL_SPECULAR, GLfloat4(1,1,1,1))
        gl.glDisable(gl.GL_LIGHT0)
        gl.glEnable(gl.GL_LIGHT1)

    def calculate_initial_compass_bearing(self, pointA, pointB):
        """
        Calculates direction between two points.
        Code based on compassbearing.py module
        https://gist.github.com/jeromer/2005586

        pointA: latitude/longitude for first point (decimal degrees)
        pointB: latitude/longitude for second point (decimal degrees)
    
        Return: direction heading in degrees (0-360 degrees, with 90 = North)
        """

        if (type(pointA) != tuple) or (type(pointB) != tuple):
            raise TypeError("Only tuples are supported as arguments")
    
        lat1 = math.radians(pointA[0])
        lat2 = math.radians(pointB[0])
    
        diffLong = math.radians(pointB[1] - pointA[1])
    
        # Direction angle (-180 to +180 degrees):
        # θ = atan2(sin(Δlong).cos(lat2),cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong))

        x = math.sin(diffLong) * math.cos(lat2)
        y = math.cos(lat1) * math.sin(lat2) - (math.sin(lat1) * math.cos(lat2) * math.cos(diffLong))
    
        initial_bearing = math.atan2(x, y)
    
        # Direction calculation requires to normalize direction angle (0 - 360)
        initial_bearing = math.degrees(initial_bearing)
        compass_bearing = (initial_bearing + 360) % 360
    
        return compass_bearing

def distance_to_coords_formula(latitude, longitude, bearing1, bearing2):
    """ math formula for calculating the lat, lng based on
        distance and bearing """

    # one degree of latitude is approximately 10^7 / 90 = 111,111 meters.
    # http://stackoverflow.com/questions/2187657/calculate-second-point-
    # knowing-the-starting-point-and-distance
    # one degree of latitude is approximately 10^7 / 90 = 111,111 meters
    # http://stackoverflow.com/questions/13836416/geohash-and-max-distance
    distance = 118  # meters

    east_displacement_a = distance * sin(radians(bearing1)) / 111111
    north_displacement_a = distance * cos(radians(bearing1)) / 111111

    east_displacement_b = distance * sin(radians(bearing2)) / 111111
    north_displacement_b = distance * cos(radians(bearing2)) / 111111

    # calculate the total displacement for N, S respectively
    waypoint_latitude_a = latitude + north_displacement_a
    waypoint_longitude_a = longitude + east_displacement_a

    waypoint_latitude_b = latitude + north_displacement_b
    waypoint_longitude_b = longitude + east_displacement_b

    return [(waypoint_latitude_a, waypoint_longitude_a),
            (waypoint_latitude_b, waypoint_longitude_b)]

 def distance(a, b):
     """Calculates distance between two latitude-longitude coordinates."""
     R = 3963  # radius of Earth (miles)
     lat1, lon1 = math.radians(a[0]), math.radians(a[1])
     lat2, lon2 = math.radians(b[0]), math.radians(b[1])
     return math.acos( math.sin(lat1)*math.sin(lat2) +
         math.cos(lat1)*math.cos(lat2)*math.cos(lon1-lon2) ) * R

def create_circle(cx, cy, r, a1, a2):
    points = []
    for a in range(a1, a2 + 1):
        x = cx + r * cos(radians(a))
        y = cy + r * sin(radians(a))
        points.append((x, y))
    return points

def dayLength(date, latitude):
    # Formulas from: http://www.gandraxa.com/length_of_day.xml
    # I don't really understand them, and the numbers don't quite match
    # what I get from calendar sites. Perhaps this is measuring the exact
    # moment the center of the sun goes down, as opposed to civil twilight?
    # But I think it's close enough to get the idea across.

    # Tilt of earth's axis relative to its orbital plane ("obliquity of ecliptic")
    axis = math.radians(23.439)

    # Date of winter solstice in this year. Not quite right, but good
    # enough for our purposes.
    solstice = date.replace(month=12, day=21)

    # If a year is a full circle, this is the angle between the solstice
    # and this date, in radians. May be negative if we haven't reached the
    # solstice yet.
    dateAngle = (date - solstice).days * 2 * math.pi / 365.25

    latitude = math.radians(latitude)
    m = 1 - math.tan(latitude) * math.tan(axis * math.cos(dateAngle))

    # If m is less than zero, the sun never rises; if greater than two, it never sets.
    m = min(2, max(0, m))
    return math.acos(1 - m) / math.pi

    def set_3d(self):
        """ Configure OpenGL to draw in 3d.

		"""
        width, height = self.get_size()
        glEnable(GL_DEPTH_TEST)
        glViewport(0, 0, width, height)
        glMatrixMode(GL_PROJECTION)
        glLoadIdentity()
        # gluPerspective(45.0f,(GLfloat)width/(GLfloat)height,near distance,far distance);
        # gluPerspective(65.0, width / float(height), 0.1, 60.0)
        gluPerspective(65.0, width / float(height), 0.1, 6000.0)
        glMatrixMode(GL_MODELVIEW)
        glLoadIdentity()

        # glDepthMask(GL_FALSE)
        # drawSkybox()
        # glDepthMask(GL_TRUE)

        x, y = self.rotation
        # http://wiki.delphigl.com/index.php/glRotate
        # procedure glRotatef(angle: TGLfloat; x: TGLfloat; y: TGLfloat; z: TGLfloat);
        glRotatef(x, 0, 1, 0)
        glRotatef(-y, math.cos(math.radians(x)), 0, math.sin(math.radians(x)))
        x, y, z = self.position
        glTranslatef(-x, -y, -z)

    def precisionToDim(self):
        """Convert precision from Wikibase to GeoData's dim and return the latter.

        dim is calculated if the Coordinate doesn't have a dimension, and precision is set.
        When neither dim nor precision are set, ValueError is thrown.

        Carrying on from the earlier derivation of precision, since
        precision = math.degrees(dim/(radius*math.cos(math.radians(self.lat)))), we get
            dim = math.radians(precision)*radius*math.cos(math.radians(self.lat))
        But this is not valid, since it returns a float value for dim which is an integer.
        We must round it off to the nearest integer.

        Therefore::
            dim = int(round(math.radians(precision)*radius*math.cos(math.radians(self.lat))))

        @rtype: int or None
        """
        if self._dim is None and self._precision is None:
            raise ValueError('No values set for dim or precision')
        if self._dim is None and self._precision is not None:
            radius = 6378137
            self._dim = int(
                round(
                    math.radians(self._precision) * radius * math.cos(math.radians(self.lat))
                )
            )
        return self._dim

def bird_blave(timestamp, place, tilt=0, azimuth=180, cloudCover=0.0):
    #SUN Position
    o = ephem.Observer()
    o.date = timestamp #'2000/12/21 %s:00:00' % (hour - self.tz)
    latitude, longitude = place
    o.lat = radians(latitude)
    o.lon = radians(longitude)
    az = ephem.Sun(o).az
    alt = ephem.Sun(o).alt

    #Irradiance
    day = dayOfYear(timestamp)
    record = {}
    record['utc_datetime'] = timestamp
    Z = pi/2-alt
    aaz = radians(azimuth+180)
    slope = radians(tilt)

    #incidence angle
    theta = arccos(cos(Z)*cos(slope) + \
            sin(slope)*sin(Z)*cos(az - pi - aaz))
    ETR = apparentExtraterrestrialFlux(day)
    #pressure?
    t, Bh, Gh = bird(theta, 1010.0)
    Dh = diffuseHorizontal(alt, Bh, day)

    record['DNI (W/m^2)'] = Bh #8 Direct normal irradiance
    record['GHI (W/m^2)'] = Gh #5 Global horizontal irradiance
    record['DHI (W/m^2)'] = Dh #11 Diffuse horizontal irradiance
    record['ETR (W/m^2)'] = ETR
    return record

	def getDec(self,j):
		j=float(j)
		m = 357.0+(0.9856*j) 
		c = (1.914*sin(radians(m))) + (0.02*sin(radians(2.0*m)))
		l = 280.0 + c + (0.9856*j)
		sinDec= 0.3978*sin(radians(l))
		return degrees(asin(sinDec))

 def _convert_to_cartesian(self, location):
     lat = radians(location.lat)
     lon = radians(location.long)
     x = cos(lat) * cos(lon)
     y = cos(lat) * sin(lon)
     z = sin(lat)
     return (x, y, z)

    def __init__(self,figura_1, figura_2, figura_1_punto=(0,0), figura_2_punto=(0,0), angulo_minimo=None,angulo_maximo=None, fisica=None, con_colision=True):
        """ Inicializa la constante
        :param figura_1: Una de las figuras a conectar por la constante.
        :param figura_2: La otra figura a conectar por la constante.
        :param figura_1_punto: Punto de rotación de figura_1
        :param figura_2_punto: Punto de rotación de figura_2
        :param angulo_minimo: Angulo minimo de rotacion para figura_2 con respecto a figura_1_punto
        :param angulo_maximo: Angulo maximo de rotacion para figura_2 con respecto a figura_1_punto
        :param fisica: Referencia al motor de física.
        :param con_colision: Indica si se permite colisión entre las dos figuras.
        """
        if not fisica:
            fisica = pilas.escena_actual().fisica

        if not isinstance(figura_1, Figura) or not isinstance(figura_2, Figura):
            raise Exception("Las dos figuras tienen que ser objetos de la clase Figura.")

        constante = box2d.b2RevoluteJointDef()
        constante.Initialize(bodyA=figura_1._cuerpo, bodyB=figura_2._cuerpo,anchor=(0,0))
        constante.localAnchorA = convertir_a_metros(figura_1_punto[0]), convertir_a_metros(figura_1_punto[1])
        constante.localAnchorB = convertir_a_metros(figura_2_punto[0]), convertir_a_metros(figura_2_punto[1])       
        if angulo_minimo != None or angulo_maximo != None:
            constante.enableLimit = True
            constante.lowerAngle = math.radians(angulo_minimo)
            constante.upperAngle = math.radians(angulo_maximo)
        constante.collideConnected = con_colision
        self.constante = fisica.mundo.CreateJoint(constante)

def vert_sep():
    """This function collects all of the coordinates and elevations necessary to calculate the vertical separation between a drone being flown under
    an established approach path, runs all of the calculations, and prints the separation in feet. All coordinates are entered by the user in degrees
    minutes seconds format with no special characters. Degrees minutes and seconds are separated with spaces. All elevations and altitudes are entered
    in feet. The approach path glide slope angle is entered in degrees. The example coordinates, elevations, and altitudes are
    for Gainesville Runway 5 ILS and RNAV approaches."""
    drone_northing_dms = raw_input("Enter the northing of the location that the drone will fly in DMS format. ex. 34 14 21.25 :")
    drone_northing_dd = dms_to_dd(drone_northing_dms)
    drone_easting_dms = raw_input("Enter the easting of the location that the drone will fly in DMS format. ex. 83 52 1.75 :")
    drone_easting_dd = dms_to_dd(drone_easting_dms)
    iaf_northing_dms = raw_input("Enter the northing of the initial approach fix in DMS format. ex. 34 12 12.11 :")
    iaf_northing_dd = dms_to_dd(iaf_northing_dms)
    iaf_easting_dms = raw_input("Enter the easting of the initial approach fix in DMS format. ex. 83 54 22.74 :")
    iaf_easting_dd = dms_to_dd(iaf_easting_dms)
    
    long_degree_dist_ft = (math.cos(math.radians(float(drone_northing_dd))) * 69.172) * 5280
    easting_distance_ft = (abs(float(iaf_easting_dd)) - float(drone_easting_dd)) * float(long_degree_dist_ft)
    northing_distance_ft = (abs(float(iaf_northing_dd) - float(drone_northing_dd))) * 364320
    distance_ft = ((float(easting_distance_ft) ** 2.0) + (float(northing_distance_ft) ** 2.0)) ** (1/2.0)
    tdze = raw_input("Enter the elevation of the touch down zone in feet. ex 1276 :")
    iaf_alt = raw_input("Enter the altitude above mean sea level of the initial approach fix in feet. ex 3100 :")
    gsa = raw_input("Enter the glide slope angle in degrees. ex 3.00 :")
    study_area_elevation = raw_input("Enter the elevation of the location that the drone will fly in feet. ex. 1162 :")
    drone_alt = raw_input("Enter the altitude above ground level that the drone will fly in feet. ex 270 :")
    adjactnt_side = float(iaf_alt) - float(tdze)
    alt_above_tdze = (math.tan(math.radians(float(gsa)))) * float(distance_ft)
    elevation_diff = float(tdze) - float(study_area_elevation)
    approach_path_alt = float(alt_above_tdze) + float(elevation_diff)
    vert_sep = float(approach_path_alt) - float(drone_alt)
    print ("There will be {} feet of vertical separation between the drone and the approach path.").format(round(vert_sep, 1))

def distance(origin, destination):
    """ 
    Calculates both distance and bearing
    """
    lat1, lon1 = origin
    lat2, lon2 = destination
    if lat1>1000:
        (lat1,lon1)=dm2dd(lat1,lon1)
        (lat2,lon2)=dm2dd(lat2,lon2)
        print('converted to from ddmm to dd.ddd')
    radius = 6371 # km
    

    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c
    
    def calcBearing(lat1, lon1, lat2, lon2):
       dLon = lon2 - lon1
       y = math.sin(dLon) * math.cos(lat2)
       x = math.cos(lat1) * math.sin(lat2) \
           - math.sin(lat1) * math.cos(lat2) * math.cos(dLon)
       return math.atan2(y, x)
       
    bear= math.degrees(calcBearing(lat1, lon1, lat2, lon2))  
    return d,bear

def xy(lng, lat, truncate=False):
    """Convert longitude and latitude to web mercator x, y

    Parameters
    ----------
    lng, lat : float
        Longitude and latitude in decimal degrees.
    truncate : bool, optional
        Whether to truncate or clip inputs to web mercator limits.

    Returns
    -------
    x, y : float
        y will be inf at the North Pole (lat >= 90) and -inf at the
        South Pole (lat <= -90).
    """
    if truncate:
        lng, lat = truncate_lnglat(lng, lat)
    x = 6378137.0 * math.radians(lng)
    if lat <= -90:
        y = float('-inf')
    elif lat >= 90:
        y = float('inf')
    else:
        y = 6378137.0 * math.log(
            math.tan((math.pi * 0.25) + (0.5 * math.radians(lat))))
    return x, y