def find_perp_error(self):
        global path_list
        xall = path_list[0]
        zall = path_list[1]
        kmin = 1000.0 #arbitrarily high value

        self.delta = 1
        self.eps = -tan(pi/2 - self.thetac)
        self.zeta = -self.xc - self.zc*tan(pi/2 - self.thetac)

        a1 = self.alpha 
        b1 = self.beta
        c1 = self.gamma
        a2 = self.delta
        b2 = self.eps
        c2 = self.zeta

        xdm = (b2*c1 - b1*c2)/ fabs(a1*b2 - a2*b1)
        zdm = (-a2*c1 + a1*c2)/ fabs(a1*b2 - a2*b1)
               
       #Perpendicular error message
        self.perppt.x = xdm
        self.perppt.z = zdm
        self.perppt.y = 0.0

       #Perpendicular error 
        xde = self.lateralpt.x
        zde = self.lateralpt.z
        
        kerr = sqrt((xde-xdm)**2 + (zde-zdm)**2)

        return float(kerr)

def s3r_ab(x1, y1, theta1, x2, y2, Db, Da, t1) :
    theta_seg1 = math.atan2(y2 - y1, x2 - x1)
    
    delta_theta1 = theta1 - theta_seg1
    delta_theta1 = modulo_angle(delta_theta1)
        
    if delta_theta1 > 0.0 :
        theta1_a = theta1 - pi/2.0
    else :
        theta1_a = theta1 + pi/2.0
    theta1_a = modulo_angle(theta1_a)
        
    [s_bx, s_by] = signe_delta_xy(theta1)
    #~ print(theta1)
    print("s_bx: {0}, s_by: {1}".format(s_bx, s_by))
    
    [s_ax, s_ay] = signe_delta_xy(theta1_a)
    
    bx = s_bx * ((Db/t1) / sqrt(1.0 + pow(tan(theta1), 2)))
    by = s_by * sqrt(pow(Db/t1, 2) - pow(bx, 2))
    
    ax = s_ax * (((6.0*Da)/pow(t1, 2)) / sqrt(1.0 + pow(tan(theta1_a), 2)))
    ay = s_ay * sqrt(pow((6.0*Da)/pow(t1, 2), 2) - pow(ax, 2))
    
    return ([ax, ay, bx, by])

def planetstransit(date):
    """returns SHA and meridian passage for the navigational planets
    """
    v = ephem.Venus()
    mars = ephem.Mars()
    j = ephem.Jupiter()
    sat = ephem.Saturn()
    obs = ephem.Observer()

    obs.date = date
    v.compute(date)
    vsha = nadeg(2 * math.pi - ephem.degrees(v.g_ra).norm)
    vtrans = time(obs.next_transit(v))
    hpvenus = "%0.1f" % ((math.tan(6371 / (v.earth_distance * 149597870.7))) * 60 * 180 / math.pi)

    obs.date = date
    mars.compute(date)
    marssha = nadeg(2 * math.pi - ephem.degrees(mars.g_ra).norm)
    marstrans = time(obs.next_transit(mars))
    hpmars = "%0.1f" % ((math.tan(6371 / (mars.earth_distance * 149597870.7))) * 60 * 180 / math.pi)

    obs.date = date
    j.compute(date)
    jsha = nadeg(2 * math.pi - ephem.degrees(j.g_ra).norm)
    jtrans = time(obs.next_transit(j))

    obs.date = date
    sat.compute(date)
    satsha = nadeg(2 * math.pi - ephem.degrees(sat.g_ra).norm)
    sattrans = time(obs.next_transit(sat))

    return [vsha, vtrans, marssha, marstrans, jsha, jtrans, satsha, sattrans, hpmars, hpvenus]

def apply_2d_transforms(context, box):
    # "Transforms apply to block-level and atomic inline-level elements,
    #  but do not apply to elements which may be split into
    #  multiple inline-level boxes."
    # http://www.w3.org/TR/css3-2d-transforms/#introduction
    if box.style.transform and not isinstance(box, boxes.InlineBox):
        border_width = box.border_width()
        border_height = box.border_height()
        origin_x, origin_y = box.style.transform_origin
        origin_x = percentage(origin_x, border_width)
        origin_y = percentage(origin_y, border_height)
        origin_x += box.border_box_x()
        origin_y += box.border_box_y()

        context.translate(origin_x, origin_y)
        for name, args in box.style.transform:
            if name == 'scale':
                context.scale(*args)
            elif name == 'rotate':
                context.rotate(args)
            elif name == 'translate':
                translate_x, translate_y = args
                context.translate(
                    percentage(translate_x, border_width),
                    percentage(translate_y, border_height),
                )
            else:
                if name == 'skewx':
                    args = (1, 0, math.tan(args), 1, 0, 0)
                elif name == 'skewy':
                    args = (1, math.tan(args), 0, 1, 0, 0)
                else:
                    assert name == 'matrix'
                context.transform(cairo.Matrix(*args))
        context.translate(-origin_x, -origin_y)

    def recalcCameraSphere(self):
        nearPlaneDist = base.camLens.getNear()
        hFov = base.camLens.getHfov()
        vFov = base.camLens.getVfov()
        hOff = nearPlaneDist * math.tan(deg2Rad(hFov / 2.0))
        vOff = nearPlaneDist * math.tan(deg2Rad(vFov / 2.0))
        camPnts = [Point3(hOff, nearPlaneDist, vOff),
         Point3(-hOff, nearPlaneDist, vOff),
         Point3(hOff, nearPlaneDist, -vOff),
         Point3(-hOff, nearPlaneDist, -vOff),
         Point3(0.0, 0.0, 0.0)]
        avgPnt = Point3(0.0, 0.0, 0.0)
        for camPnt in camPnts:
            avgPnt = avgPnt + camPnt

        avgPnt = avgPnt / len(camPnts)
        sphereRadius = 0.0
        for camPnt in camPnts:
            dist = Vec3(camPnt - avgPnt).length()
            if dist > sphereRadius:
                sphereRadius = dist

        avgPnt = Point3(avgPnt)
        self.ccSphereNodePath.setPos(avgPnt)
        self.ccSphereNodePath2.setPos(avgPnt)
        self.ccSphere.setRadius(sphereRadius)

def uproj_tmerc(x, y):
    easting = x - x0
    northing = y
    # Meridional Arc
    M = northing / k0
    mu = M / (a * (1.0 - e * e / 4.0 - 3.0 * e * e * e * e / 64.0 -
        5.0 * e * e * e * e * e * e / 256.0))

    # Footprint latitude
    fp = mu + J1 * math.sin(2.0 * mu) + J2 * math.sin(4.0 * mu) + \
            J3 * math.sin(6.0 * mu) + J4 * math.sin(8.0 * mu)
    C1 = ep2 * math.cos(fp) * math.cos(fp)
    T1 = math.tan(fp) * math.tan(fp)
    # Radius of curvature in meridian plane
    R1 = a * (1.0 - e * e) / \
            math.pow(1.0 - e * e * math.sin(fp) * math.sin(fp), 1.5)
    # Radius of curvature perpendicular to meridian plane
    N1 = a / math.sqrt(1.0 - e * e * math.sin(fp) * math.sin(fp))
    D = easting / (N1 * k0)

    Q1 = N1 * math.tan(fp) / R1
    Q2 = D * D / 2.0
    Q3 = (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C1 * C1 - 9.0 * ep2) * \
            D * D * D * D / 24.0
    Q4 = (61.0 + 90.0 * T1 + 298.0 * C1 + 45.0 * T1 * T1 - 3.0 * C1 * C1 -
            252.0 * ep2) * D * D * D * D * D * D / 720.0
    Q5 = D
    Q6 = (1.0 + 2.0 * T1 + C1) * D * D * D / 6.0
    Q7 = (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C1 * C1 + 8.0 * ep2 +
            24.0 * T1 * T1) * D * D * D * D * D / 120.0

    lat = fp - Q1 * (Q2 - Q3 + Q4)
    lon = math.radians(long0) + (Q5 - Q6 + Q7) / math.cos(fp)

    return ( math.degrees(lat), math.degrees(lon) )

	def mapcoords(self, center_lat, center_lon, zoomlevel, width, height):
		center_lon_rad = math.radians(center_lon);
		center_lat_rad = math.radians(center_lat);
		pos_lon_rad = math.radians(self.lon);
		pos_lat_rad = math.radians(self.lat);
		n = pow(2.0, zoomlevel);

		centerTileX = ((center_lon + 180) / 360) * n;
		centerTileY = (1 - (math.log(math.tan(center_lat_rad) + 1.0/math.cos(center_lat_rad)) / math.pi)) * n / 2.0;
		
		#print("centerTileX = {0}".format(centerTileX))
		#print("centerTileY = {0}".format(centerTileY))
		
		posTileX = ((self.lon + 180) / 360) * n;
		posTileY = (1 - (math.log(math.tan(pos_lat_rad) + 1.0/math.cos(pos_lat_rad)) / math.pi)) * n / 2.0;
		
		#print("posTileX = {0}".format(posTileX))
		#print("posTileY = {0}".format(posTileY))
		
		diffTileX = posTileX - centerTileX
		diffTileY = posTileY - centerTileY
		
		#print("diffTileX = {0}".format(diffTileX))
		#print("diffTileY = {0}".format(diffTileY))
		
		diffPixlesX = 256 * diffTileX
		diffPixlesY = 256 * diffTileY
		
		#print("diffPixlesX = {0}".format(diffPixlesX))
		#print("diffPixlesY = {0}".format(diffPixlesY))
		
		xPos = (width / 2) + diffPixlesX
		yPos = (height / 2) + diffPixlesY
		
		return [int(xPos), int(yPos)]

	def getZD(self, md, lat, decl, umd):
		'''Calculates Regiomontan zenith distance '''

		zd = 0.0
		if md == 90.0:
			zd = 90.0-math.degrees(math.atan(math.sin(math.fabs(math.radians(lat))))*math.tan(math.radians(decl)))
		elif md < 90.0:
			A = math.degrees(math.atan(math.cos(math.radians(lat))*math.tan(math.radians(md))))
			B = math.degrees(math.atan(math.tan(math.fabs(math.radians(lat)))*math.cos(math.radians(md))))

			C = 0.0
			if (decl < 0 and lat < 0) or (decl >= 0 and lat >= 0):
				if umd:
					C = B-math.fabs(decl)
				else:
					C = B+math.fabs(decl)
			elif (decl < 0 and lat > 0) or (decl > 0 and lat < 0):
				if umd:
					C = B+math.fabs(decl)
				else:
					C = B-math.fabs(decl)

			F = math.degrees(math.atan(math.sin(math.fabs(math.radians(lat)))*math.sin(math.radians(md))*math.tan(math.radians(C)))) #C and F can be negative
			zd = A+F

		return zd

	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 tex_triangle(tex, sommet, mesure, tex_mesure):
    # Coefficient tel que hauteur = gamma * base
    gamma = math.tan(math.radians(mesure[0]))*math.tan(math.radians(mesure[1]))/(math.tan(math.radians(mesure[0]))+math.tan(math.radians(mesure[1])))
    # Coordonnée des sommets pour que la base ou la hauteur mesure (max-min) cm
    min = -4
    max = 4
    if gamma < 1:
        coordonnee = [min, min, max, min, round((max-min)*gamma/math.tan(math.radians(mesure[0]))+min, 4),round((max-min)*gamma+min,4)]
    else:
        coordonnee = [min, min, round((max-min)/gamma+min, 4), min, round((max-min)/math.tan(math.radians(mesure[0])) +min, 4) ,max]
    # Angle des symboles des angles
    angle = [0, mesure[0], 180 - mesure[1], 180, 180 + mesure[0], 180 + mesure[0] + mesure[2]]
    ## Construction
    tex.append("\\begin{center}")
    tex.append("\\psset{unit=0.5cm}")
    tex.append("\\begin{pspicture}(%s,%s)(%s,%s)" %(coordonnee[0]-1, coordonnee[1]-1, coordonnee[2]+1, coordonnee[5]+1))
    # Triangle
    tex.append("\\pstTriangle(%s,%s){%s}(%s,%s){%s}(%s,%s){%s}" %(coordonnee[0], coordonnee[1], sommet[0], coordonnee[2], coordonnee[3], sommet[1], coordonnee[4], coordonnee[5], sommet[2]))
    # Symbole et légende de chaque angle
    for i in range(len(mesure)):
        if mesure[i] == 90:
           tex.append("\\pstRightAngle{%s}{%s}{%s}" % ( sommet[(i+1)%3], sommet[i], sommet[(i-1)%3]))
        elif mesure[0] == mesure[1] and i < 2:
           tex.append("\\pstMarkAngle{%s}{%s}{%s}{}" % ( sommet[(i+1)%3], sommet[i], sommet[(i-1)%3]))
           tex.append("\\uput{0.25}[%s]{%s}(%s,%s){\\psline(0,0)(0.5,0)}" % ((angle[2*i]+angle[2*i+1])/2, (angle[2*i]+angle[2*i+1])/2, coordonnee[2*i], coordonnee[2*i+1]))
           if i == 0:
               tex.append("\\pstMarkAngle{%s}{%s}{%s}{$%s$}" % ( sommet[(i+1)%3], sommet[i], sommet[(i-1)%3], tex_mesure[i]))
        else:
           tex.append("\\pstMarkAngle{%s}{%s}{%s}{$%s$}" % ( sommet[(i+1)%3], sommet[i], sommet[(i-1)%3], tex_mesure[i]))
    tex.append("\\end{pspicture}")
    tex.append("\\end{center}")

    def _convert_transform_to_tikz(self, transform):
        """Convert a SVG transform attribute to a list of TikZ transformations"""
        #return ""
        if not transform:
            return []

        options = []
        for cmd, params in transform:
            if cmd == 'translate':
                x, y = params
                options.append("shift={(%s,%s)}" % (round(x, 5) or '0', round(y, 5) or '0'))
                # There is bug somewere.
                # shift=(400,0) is not equal to xshift=400

            elif cmd == 'rotate':
                if params[1] or params[2]:
                    options.append("rotate around={%s:(%s,%s)}" % params)
                else:
                    options.append("rotate=%s" % round(params[0], 5))
            elif cmd == 'matrix':
                options.append("cm={{%s,%s,%s,%s,(%s,%s)}}" % tuple(map(lambda x: round(x, 5), params)))
            elif cmd == 'skewX':
                options.append("xslant=%.3f" % math.tan(params[0] * math.pi / 180))
            elif cmd == 'skewY':
                options.append("yslant=%.3f" % math.tan(params[0] * math.pi / 180))
            elif cmd == 'scale':
                if params[0] == params[1]:
                    options.append("scale=%.3f" % params[0])
                else:
                    options.append("xscale=%.3f,yscale=%.3f" % params)

        return options

def MatrixLog6(T):#Takes a T matrix SE(3) and returns a 6vector Stheta
    '''
Example Input: 
T = [[1,0,0,0], [0,0,-1,0], [0,1,0,3], [0,0,0,1]]
Output:
[1.5707963267948966, 0.0, 0.0, 0.0, 2.3561944901923448, 2.3561944901923457]
    '''
    R,p = TransToRp(T)
    Rtrace = R[0][0]+R[1][1]+R[2][2]
    if(R==np.eye(3)).all():
        w=0
        v=Normalise(p)
        th=Magnitude(p)
    
    elif(Rtrace == -1):
        th = pi
        w = MatrixLog3(R)
        G = (1/th)*np.eye(3) - 0.5*np.asarray(VecToso3(w)) + ((1/th)-((1/(tan(th/2.0)))/2.0))*(matmult(VecToso3(w),VecToso3(w)))
        v = np.dot(G,p)

    else:
        th = acos((Rtrace-1)/2.0)
        w = so3ToVec((1/(2*np.sin(th)))*(np.subtract(R, RotInv(R))))
        G = (1/th)*np.eye(3) - 0.5*np.asarray(VecToso3(w)) + ((1/th)-((1/(tan(th/2.0)))/2.0))*(matmult(VecToso3(w),VecToso3(w)))
        v = np.dot(G,p)     

    return ([w[0]*th,w[1]*th,w[2]*th,v[0]*th,v[1]*th,v[2]*th])

def k_eq(R, H, beta_deg):
    k=0.
    beta = math.radians(beta_deg)
    tan_35 = math.pow(math.tan(beta), .35)
    tan_23 = math.pow(math.tan(beta), .23)
    k=R/H*1.15/tan_35*math.pow(H/(2.*R), .9/tan_23)
    return k

def ground_offset(height, pitch, roll, yaw):
    '''
    find the offset on the ground in meters of the center of view of the plane
    given height above the ground in meters, and pitch/roll/yaw in degrees.

    The yaw is from grid north. Positive yaw is clockwise
    The roll is from horiznotal. Positive roll is down on the right
    The pitch is from horiznotal. Positive pitch is up in the front

    return result is a tuple, with meters east and north of GPS position

    This is only correct for small values of pitch/roll
    '''

    # x/y offsets assuming the plane is pointing north
    xoffset = -height * math.tan(math.radians(roll))
    yoffset = height * math.tan(math.radians(pitch))

    # convert to polar coordinates
    distance = math.hypot(xoffset, yoffset)
    angle    = math.atan2(yoffset, xoffset)

    # add in yaw
    angle -= math.radians(yaw)

    # back to rectangular coordinates
    x = distance * math.cos(angle)
    y = distance * math.sin(angle)

    return (x, y)

def determine_if_in_wake(xt, yt, xw, yw, k, r0, alpha):  # According to Jensen Model only
    # Eq. of centreline is Y = tan (d) (X - Xt) + Yt
    # Distance from point to line
    alpha = deg2rad(alpha + 180)
    distance_to_centre = abs(- tan(alpha) * xw + yw + tan(alpha) * xt - yt) / sqrt(1.0 + tan(alpha) ** 2.0)
        # print distance_to_centre
    # Coordinates of the intersection between closest path from turbine in wake to centreline.
    X_int = (xw + tan(alpha) * yw + tan(alpha) * (tan(alpha) * xt - yt)) / (tan(alpha) ** 2.0 + 1.0)
    Y_int = (- tan(alpha) * (- xw - tan(alpha) * yw) - tan(alpha) * xt + yt) / (tan(alpha) ** 2.0 + 1.0)
    # Distance from intersection point to turbine
    distance_to_turbine = sqrt((X_int - xt) ** 2.0+(Y_int - yt) ** 2.0)
    # Radius of wake at that distance
    radius = wake_radius(r0, k, distance_to_turbine)
    # print radius
    if (xw - xt) * cos(alpha) + (yw - yt) * sin(alpha) <= 0.0:
        if abs(radius) >= abs(distance_to_centre):
            if abs(radius) >= abs(distance_to_centre) + r0:
                fraction = 1.0
                value = True
                return fraction
            elif abs(radius) < abs(distance_to_centre) + r0:
                fraction = area.AreaReal(r0, radius, distance_to_centre).area()
                value = True
                return fraction
        elif abs(radius) < abs(distance_to_centre):
            if abs(radius) <= abs(distance_to_centre) - r0:
                fraction = 0.0
                value = False
                return fraction
            elif abs(radius) > abs(distance_to_centre) - r0:
                fraction = area.AreaReal(r0, radius, distance_to_centre).area()
                value = True
                return fraction
    else:
        return 0.0

	def get_trafo(self):
		x0 = self.origin_x.get_point_value()
		y0 = self.origin_y.get_point_value()
		sx = self.scale_x.get_value() / 100.0
		sy = self.scale_y.get_value() / 100.0
		shx = self.shear_x.get_value()
		shy = self.shear_y.get_value()

		if shx + shy > 85:
			if shx == self.transforms[3]: shy = 85 - shx
			else: shx = 85 - shy

		shx = math.pi * shx / 180.0
		shy = math.pi * shy / 180.0

		angle = math.pi * self.rotate.get_value() / 180.0

		trafo = [sx, 0.0, 0.0, sy, x0, y0]
		if angle:
			trafo2 = [math.cos(angle), math.sin(angle),
					- math.sin(angle), math.cos(angle), 0.0, 0.0]
			trafo = libgeom.multiply_trafo(trafo, trafo2)
		if shx or shy:
			trafo2 = [1.0, math.tan(shy), math.tan(shx), 1.0, 0.0, 0.0]
			trafo = libgeom.multiply_trafo(trafo, trafo2)

		self.transforms = [sx, sy, shx, shy, angle]
		return trafo, [sx, sy, shx, shy, angle]

	def setupKinectParams(self, version=1):
		self.kinect_hfov = 1.0140363  # 58.1 (degrees)
		self.kinect_vfov = 0.813323431  # 46.6 (degrees)
		self.kinect_height = 240
		self.kinect_width = 320

		if version == 2:  # Kinect v2
			self.kinect_hfov = 1.23220245  # 70.6 (degrees)
			self.kinect_vfov = 1.04719755  # 60 (degrees)
			self.kinect_height = 424
			self.kinect_width = 512

		projectionScale = [math.tan((self.kinect_hfov * 0.5)) / (self.kinect_width * 0.5),\
						   math.tan((self.kinect_vfov * 0.5)) / (self.kinect_height * 0.5)]


		imgCenter = [float(self.kinect_width) * 0.5, float(self.kinect_height * 0.5)]

		xIdxImg = np.tile(np.arange(self.kinect_width), self.kinect_height).reshape((self.kinect_height, self.kinect_width))
		xIdxImg = xIdxImg.astype(np.float16)
		yIdxImg = np.tile(np.arange(self.kinect_height), self.kinect_width).reshape((self.kinect_width, self.kinect_height))
		yIdxImg = yIdxImg.T
		yIdxImg = yIdxImg.astype(np.float16)

		self.xIdxImg = ((xIdxImg - imgCenter[0]) * projectionScale[0])
		self.yIdxImg = ((yIdxImg - imgCenter[1]) * projectionScale[1])

  def Refrac(self):    # Originally in CORRECT.PAS
    """Calculate the correction for atmospheric refraction for the given coordinates.
                   
       Returns a tuple of dRA and dDEC, which are OFFSETS from the current position, in arcseconds.
    """
    ObjRa = self.RaA / 54000.0       # Translate to hours from arcsec}
    ObjDec = self.DecA / 3600.0      # Translate to degrees from arcsec}

    z = DegToRad(90 - self.Alt)                      # Zenith distance in radians}
    if z <= 0:
      z = 1e-6
    h = DegToRad((self.Time.LST - ObjRa) * 15)     # Hour angle in radians}
    dummy = trunc(h / (2 * pi))
    h -= dummy * 2 * pi
    obs = DegToRad(prefs.ObsLat)                    # Observatory Lat in radians}
    R = -1
    NewR = 0
    twiggles = 0
    # Calculate the value R in arc seconds}
    while ((NewR - R) >= 1e-8) and (twiggles < 20):
      twiggles += 1
      R = NewR
      Tanof = tan(z - DegToRad(R / 3600))
      NewR = (R1 * Tanof) + (R2 * Tanof * Tanof * Tanof)

    if twiggles > 20:     # If we're very close to the horizon, the iterative solver can fail...
      logger.error('correct.CalcPosition.Refrac: Too many Twiggles in refraction code!')

    # Calculate dRA and dDEC in arcsec}
    R = NewR
    CurlR = R * 17 * (prefs.Press * 30 / 1015.92) / (460 + ((prefs.Temp * 9 / 5) + 32))   # convert Temp and Press to F and "Hg for correction}
    dRA = CurlR * sin(h) * cosec(z) * cos(obs) * sec(DegToRad(ObjDec))
    dDEC = CurlR * ((sin(obs) * cosec(z) * sec(DegToRad(ObjDec))) - (tan(DegToRad(ObjDec)) * cot(z)))
    return dRA, dDEC

def main():
    WIDTH, HEIGHT = (300, 300)
    parser = GPXParser("steig.gpx")

    fo = open("test.svg", "w")
    gp = open("gnuplot.dat", "w")
    surface = cairo.SVGSurface (fo, WIDTH, HEIGHT)
    ctx = cairo.Context (surface)
    ctx.scale (WIDTH/1.0, HEIGHT/1.0)

    for lat, lon, ele, time in parser.tracks["Nibelungensteig on GPSies.com"]:
        ctx.line_to(lon, -180/pi*log(tan(pi/4 + lat*(pi/180)/2.0)))
        gp.write("%f %f %f\n"%(lat, lon, ele))

    ctx.set_source_rgb(0.3, 0.2, 0.5) # Solid color
    ctx.set_line_width(0.02)
    ctx.stroke()
    ctx.close_path()

    f = open("start.txt", "r")
    f.readline()

    for line in f:
        if line == '\r\n':
            break
        lat, lon = map(float, line.split('\t')[:2])
        ctx.move_to(lon, -180/pi*log(tan(pi/4 + lat*(pi/180)/2.0)))
        ctx.arc(lon, -180/pi*log(tan(pi/4 + lat*(pi/180)/2.0)), 0.0005, 0.0, 2.0*pi)
        ctx.set_source_rgb(1.0, 0.0, 0.0)
        ctx.fill_preserve()
        ctx.close_path()

    surface.finish()

 def ApparentPlace(self):
   """Calculate annual aberration (I think :-)
   
      Returns dRA and dDEC corrections as a tuple, in arcseconds.
      
      # This is taken from Astronomical Formulae for Calculators, Jean Meeus,
      # 3rd Ed. 1985.  P:71-73.
   """
   Ra = DegToRad((self.RaA / 54000.0) * 15)     # Convert to degrees, and: radians
   Dec = DegToRad(self.DecA / 3600.0)         # Convert to degrees, and: radians
   if abs((Dec / 3600) + 90) < 1e-6:
     Dec = -89.999999 * 3600
   if abs((Dec / 3600) - 90) < 1e-6:
     Dec = 89.999999
   T = (self.Time.JD - 2415020) / 36525
   dPhi, dEpsi = self.Nutation(T)
   L = 279.69668 + (36000.76892 * T) + (0.0003025 * T * T)         # Sun's mean longitude
   M = 358.47583 + (35999.04975 * T) - (0.000150 * T * T) - (0.0000033 * T * T * T)
   L = DegToRad(Reduce(L))                  # Reduce to 0-360, convert to radians
   M = DegToRad(Reduce(M))
   Epsi = 23.452294 - (0.0130125 * T) - (0.00000164 * T * T) + (0.000000503 * T * T * T)
   Epsi = DegToRad(Epsi)                    # Convert to radians
   C = (((1.919460 - (0.004789 * T) - (0.000014 * T * T)) * sin(M)) +
        ((0.020094 - (0.000100 * T)) * sin(2 * M)) + (0.000293 * sin(3 * M)))
   Sun = L + DegToRad(C)                    # Sun's true longitude, in radians
   dRa1 = ((cos(Epsi) + (sin(Epsi) * sin(Ra) * tan(Dec))) * dPhi) - (cos(Ra) * tan(Dec) * dEpsi)
   dDec1 = (sin(Epsi) * cos(Ra) * dPhi) + (sin(Ra) * dEpsi)
   dRa2 = -20.49 * (((cos(Ra) * cos(Sun) * cos(Epsi)) + (sin(Ra) * sin(Sun))) / cos(Dec))
   dDec2 = -20.49 * ((cos(Sun) * cos(Epsi) * (tan(Epsi) * cos(Dec)) - (sin(Ra) * sin(Dec)))
                     + (cos(Ra) * sin(Dec) * sin(Sun)))
   dRA = dRa1 + dRa2                          # In arcsecs
   dDEC = dDec1 + dDec2                       # Also in arcsecs
   self.RaA += dRA
   self.DecA += dDEC