def _gen_texture_coords(self, tex_coords=None, mapping_type=TextureMappingConstants.SPHERICAL):
     
     if len(tex_coords) == 0 :
         tex_coords = []
         fast = False
         
         if mapping_type == TextureMappingConstants.SPHERICAL:
             
             if fast:
                 for norm in self._normals:
                     u = norm.x / 2 + 0.5
                     v = norm.z / 2 + 0.5
                     tex_coords.append([u,v])
             else:
                 for norm in self._normals:
                     u = math.asin(norm.x)/math.pi + 0.5
                     v = math.asin(norm.z)/math.pi + 0.5
                     tex_coords.append([u,v])
     
     else:
         # Expand the texture coordinates array to match the format of the transformed vertex array
         new_tex_coords = []
         
         for vi in range(len(self._vert_index)):
             
             v_ind = self._vert_index[vi]
             new_tex_coords.append(tex_coords[v_ind])
         
         tex_coords = new_tex_coords
                 
     self._tex_coords = tex_coords

    def get_pos_command(self):
        ts  = ovr.getTrackingState(self.session, ovr.getTimeInSeconds(), True)
        if ts.StatusFlags & (ovr.Status_OrientationTracked | ovr.Status_PositionTracked):
            pose = ts.HeadPose
        
        # Get queternions
        q0 = pose.ThePose.Orientation.w;
        q1 = pose.ThePose.Orientation.x;
        q2 = pose.ThePose.Orientation.y;
        q3 = pose.ThePose.Orientation.z;

        # Calculate Euler angles
        x = int(1500 + 637*math.asin(2*(q0*q1 - q3*q2)))    #pitch
        y = int(1450 - 637*math.atan2(2*(q0*q2 - q1*q3),1-2*(q2*q2+q1*q1))) #yaw
        z = int(1500 + 637*math.asin(2*(q0*q3 + q1*q2)))    #roll
        
        pos_command = ""
        # Anti-shaking
        if abs(self.xtemp-x)>1 or abs(self.ytemp-y)>1 or abs(self.ztemp-z)>1:
            pos_command = (str(x) + "," + str(y) + ","+str(z)+'\n')
            sys.stdout.flush()
            self.xtemp = x
            self.ytemp = y
            self.ztemp = z
        return pos_command

def day_length(doy, yr_days, latitude):
    """ Daylength in hours

    Eqns come from Leuning A4, A5 and A6, pg. 1196
    
    Reference:
    ----------
    Leuning et al (1995) Plant, Cell and Environment, 18, 1183-1200.
    
    Parameters:
    -----------
    doy : int
        day of year, 1=jan 1
    yr_days : int
        number of days in a year, 365 or 366
    latitude : float
        latitude [degrees]

    Returns:
    --------
    dayl : float
        daylength [hrs]

    """
    deg2rad = pi / 180.0
    latr = latitude * deg2rad
    sindec = -sin(23.5 * deg2rad) * cos(2.0 * pi * (doy + 10.0) / yr_days)
    a = sin(latr) * sindec
    b = cos(latr) * cos(asin(sindec))
    dayl = 12.0 * (1.0 + (2.0 / pi) * asin(a / b))
    
    return dayl

 def getAngle(self,cc):
         r = ( (cc[2] - cc[0]) ** 2 + (cc[3] - cc[1] ) ** 2 ) ** 0.5
         print('current R is:%8.4f' % r)
         print (cc)
         x = cc[2] - cc[0]
         y = cc[3] - cc[1]
         print(x,y)
         sina = y / r
         cosa = x / r
         starta = float()
         if y > 0 and x > 0:
             starta=math.asin ( sina ) ## 0 ~ PI/2
             print ("q1: %f" % starta)
         else:
             if  y > 0 and x <=0:
                 starta = math.acos(cosa) ## PI/2 ~ PI
                 print ("q2: %f" % starta)
             else:
                 if y <=0 and x <=0:
                     starta = math.pi*2 - math.acos(cosa)
                     print("q3: %f" % starta)
                 else:
                     starta = math.pi * 2 + math.asin(sina)
                     print("q4: %f" % starta)
         print("The angle of the hand is: %f" % starta)
         print('The coords of the hand are (after calculated from the Angle: is %s' % cc)
         return (starta)

  def _FromXYZ(coord):
    x = coord[0] * 4 - 2
    y = coord[1] * 4 - 2
    z = coord[2] * 4 - 2

    assert -2 <= x <= 2
    assert -2 <= y <= 2
    assert -2 <= z <= 2

    r = GeoApi._GetDistance((0, 0, 0), (x, y, z))
    if not r:
      return (0, 0, -GeoApi._EARTH_RADIUS)

    phi = math.asin(z / r)
    r_xy = r * math.cos(phi)
    if not r_xy:
      gamma = 0
    else:
      gamma = math.asin(y / r_xy)
    if x < 0:
      if y > 0:
        gamma = math.pi - gamma
      else:
        gamma = -math.pi - gamma

    elevation =  (r - 1) * GeoApi._EARTH_RADIUS

    phi = phi * 180.0 / math.pi
    gamma = gamma * 180.0 / math.pi

    return (phi, gamma, elevation)

def rpy(R):
    """Converts a rotation matrix to a roll,pitch,yaw angle triple.
    The result is given in radians."""
    sign = lambda x: 1 if x > 0 else (-1 if x < 0 else 0)

    m = matrix(R)
    b = -math.asin(m[2][0]) # m(2,0)=-sb
    cb = math.cos(b)
    if abs(cb) > 1e-7:
        ca = m[0][0]/cb   #m(0,0)=ca*cb
        ca = min(1.0,max(ca,-1.0))
        if sign(m[1][0]) == sign(cb): #m(1,0)=sa*cb
            a = math.acos(ca);
        else:
            a = 2*math.pi - math.acos(ca)

        cc = m[2][2] / cb  #m(2,2)=cb*cc
        cc = min(1.0,max(cc,-1.0))
        if sign(m[2][1]) == sign(cb): #m(2,1)=cb*sc
            c = math.acos(cc)
        else:
            c = math.pi*2 - math.acos(cc)
    else: 
        #b is close to 90 degrees, i.e. cb=0
        #this reduces the degrees of freedom, so we can set c=0
        c = 0
        #m(0,1)=-sa
        a = -math.asin(m[0][1]);
        if sign(math.cos(a)) != sign(m[1][1]): #m(1,1)=ca
            a = math.pi - a;
    return c,b,a

def phase3():
	print "Phase 3"
	targetAP = math.acos(.5 * 400 / 300)
	thetaAP = state.joint[ha.RAP].pos
	# change in x position of CoM during bend
	baseAR = math.asin(-88.4/600)
	deltaAR = math.asin(-88.4 / 400) - baseAR
	elapsed = 0
	#print "AP start: ", thetaAP
	#print "AP target: ", targetAP
	while (elapsed < 10):
		s.get(state, wait=False, last=False)
		t1 = state.time
		thetaAP = filterTargetPos(targetAP, elapsed, freq)
		thetaAR = filterTargetPos(deltaAR, elapsed, freq) + baseAR
		# print elapsed, "\t", thetaAR
		# shift z
		ref.ref[ha.RHP] = -thetaAP
		ref.ref[ha.RAP] = -thetaAP
		ref.ref[ha.RKN] = 2 * thetaAP
		# shift x
		ref.ref[ha.RHR] = -thetaAR
		ref.ref[ha.RAR] = thetaAR
		ref.ref[ha.LHR] = -thetaAR
		ref.ref[ha.LAR] = thetaAR
		r.put(ref)
		elapsed += TIME_STEP
		s.get(state, wait=False, last=False)
		t2 = state.time
		delay(TIME_STEP - (t2 - t1))

def CountAngle(p1, p2, p3):
    l1 = GetLength(p1, p2)
    l2 = GetLength(p2, p3)
    
    if l1 * l2 < EPS:
        return 0

    v1 = [p2[0] - p1[0], p2[1] - p1[1]]
    v2 = [p2[0] - p3[0], p2[1] - p3[1]]

    sin_arg = (v1[0] * v2[1] - v2[0] * v1[1]) / (l1 * l2)
    if abs(abs(sin_arg) - 1) < EPS:
        x = Sign(sin_arg)
        aSin = math.asin(x)
    else:
        aSin = math.asin(sin_arg)

    cos_arg = (v1[0] * v2[0] + v1[1] * v2[1]) / (l1 * l2)
    if abs(abs(cos_arg) - 1) < EPS:
        x = Sign(cos_arg)
        aCos = math.acos(x)
    else :
        aCos = math.acos(cos_arg)
    
    if sin_arg >= 0 and cos_arg >= 0: # first quarter
        return (aSin + aCos) / 2.0
    elif sin_arg >= 0 and cos_arg < 0: # second quarter
        return ((PI - aSin) + aCos) / 2.0
    elif sin_arg < 0 and cos_arg < 0: # third quarter
        return ((PI - aSin) + (2 * PI - aCos)) / 2.0
    else: # fourth quarter
        return 2 * PI + (aSin - aCos) / 2.0

    def _anglesToVirtualAngles(self, pos, wavelength):
        """
        Return dictionary of all virtual angles in radians from VliegPosition
        object win radians and wavelength in Angstroms. The virtual angles are:
        Bin, Bout, azimuth and 2theta.
        """

        # Create transformation matrices
        [ALPHA, DELTA, GAMMA, OMEGA, CHI, PHI] = createVliegMatrices(
            pos.alpha, pos.delta, pos.gamma, pos.omega, pos.chi, pos.phi)
        [SIGMA, TAU] = createVliegsSurfaceTransformationMatrices(
            self._getSigma() * TORAD, self._getTau() * TORAD)

        S = TAU * SIGMA
        y_vector = matrix([[0], [1], [0]])

        # Calculate Bin from equation 15:
        surfacenormal_alpha = OMEGA * CHI * PHI * S * matrix([[0], [0], [1]])
        incoming_alpha = ALPHA.I * y_vector
        minusSinBetaIn = dot3(surfacenormal_alpha, incoming_alpha)
        Bin = asin(bound(-minusSinBetaIn))

        # Calculate Bout from equation 16:
        #  surfacenormal_alpha has just ben calculated
        outgoing_alpha = DELTA * GAMMA * y_vector
        sinBetaOut = dot3(surfacenormal_alpha, outgoing_alpha)
        Bout = asin(bound(sinBetaOut))

        # Calculate 2theta from equation 25:

        cosTwoTheta = dot3(ALPHA * DELTA * GAMMA * y_vector, y_vector)
        twotheta = acos(bound(cosTwoTheta))
        psi = self._anglesToPsi(pos, wavelength)

        return {'Bin': Bin, 'Bout': Bout, 'azimuth': psi, '2theta': twotheta}

 def _get_turbulent_cf(self,Rex):
     tol = 1.0e-8
     r = 0.88
     Te = 222.0
     Mach = self.fc.Mach
     m = (self.gamma - 1.0)*Mach*Mach/2
     TawTe = 1.0 + r*m
     F = self.TwTaw*TawTe
     Tw = F*Te
     A = (r*m/F)**0.5
     B = (1.0+r*m-F)/F
     denom = (4*A*A + B**2)**0.5
     alpha = (2.0*A*A-B)/denom
     beta = B / denom
     if Mach>0.1:
         Fc = r*m/((asin(alpha)+asin(beta))**2)
     else:
         Fc = ((1.0+F**0.5)/2)**2
     xden = 1.0 + 122.0**(-5/Tw)/Te
     xnum = 1.0 + 122.0**(-5/Tw)/Tw
     Ftheta = xnum/xden * (1/F)**0.5
     Fx = Ftheta / Fc
     RexBar = Fx * Rex
     Cfb = 0.074/(RexBar**0.2)
     itr = 0
     err = tol+1.0
     while err>tol and itr<100:
         Cf0 = Cfb
         xnum = 0.242 - Cfb**0.5*log10(RexBar*Cfb)
         xden = 0.121 + Cfb**0.5 / log(10.0)
         Cfb = Cfb*(1.0 + xnum/xden)
         err = abs(Cf0-Cfb)
         itr += 1
     return Cfb/Fc

 def check_fisher(self, path):
     # Будем делить пополам
     L = []
     for line in open(path, "r"):
         line = line.split(" ")
         L.append(int(line[1]))
     eff = 0
     noeff = 0
     for i in range(59):
         eff += L[i]
         noeff += L[i + 59]
     eff = eff * 1500 / max(eff, noeff)
     noeff = noeff * 1500 / max(eff, noeff)
     sum_l = eff + noeff
     eff = noeff / float(eff + noeff)
     print "eff:", eff
     for name in sorted(self.S.keys()):
         eff_s = 0
         noeff_s = 0
         for i in range(59):
             eff_s += self.S[name][i]
             noeff_s += self.S[name][i + 59]
         eff_s = eff_s * 1500 / max(eff_s, noeff_s)
         noeff_s = noeff_s * 1500 / max(eff_s, noeff_s)
         sum_s = eff_s + noeff_s
         eff_s = noeff_s / float(eff_s + noeff_s)
         print "eff_s:", eff
         phi1 = asin(sqrt(eff / 100.0))
         phi2 = asin(sqrt(eff_s / 100.0))
         phi_emp = abs(phi1 - phi2) * sqrt(sum_l * sum_s / float(sum_s + sum_l))
         print name, phi_emp
         print "____________________________"

def p_tet_sp(angle_a, angle_ab, angle_ad, angle_cd):
    """ Returns an array with angles of interest (angle_b, angle_bc, angle_ac). Summer 2016. """

    if (0 > angle_a) or (angle_a > 178) or (0 > angle_ab) or (angle_ab > 178) or (0 > angle_ad) or (angle_ad > 178) or \
            (0 > angle_cd) or (angle_cd > 178):
        return exit('Does not exist 1')  # Angles of the same triangle should be at least equal to 1

    if angle_ad >= angle_ab:
        return exit('Does not exist 2')  # Angle_ab = Angle_ad + Angle_bd

    if 0 < angle_ab - angle_ad < 180:
        angle_bd = angle_ab - angle_a
    else:
        return exit('Does not exist 3')

    angle_b = math.atan(
        (math.tan(math.radians(angle_a)) * math.tan(math.radians(angle_ad))) / math.tan(math.radians(angle_bd)))

    angle_bc = math.asin(math.cos(math.radians(angle_ad)) * math.sqrt(
        ((math.tan(math.radians(angle_ad)) ** 2) + (math.sin(math.radians(angle_cd)) ** 2))))

    angle_ac = math.asin(math.cos(math.radians(angle_bd)) * math.sqrt(
        ((math.tan(math.radians(angle_bd)) ** 2) + (math.sin(math.radians(angle_cd)) ** 2))))

    if 180 <= angle_a + angle_b:
        return exit('Does not exist 4')  # Angle_a + Angle_b + Angle_c = 180

    return [angle_b, angle_bc, angle_ac]

 def angulos(self):
     angles = []
     angles.append(2 * asin(self.b / (2 * self.dist_centro_vertice)))
     angles.append(pi)
     angles.append(pi + 2 * asin(self.b / (2 * self.dist_centro_vertice)))
     angles.append(0)
     return angles

 def __draw_cell(self, x, y, is_good, percent):  # horrible graphic code,don't touch it unless you are very sure
     arc_coord = x - GroundView.CELL_RADIUS, y - GroundView.CELL_RADIUS, \
                 x + GroundView.CELL_RADIUS, y + GroundView.CELL_RADIUS
     delta_h = 2 * GroundView.CELL_RADIUS * abs(percent) - GroundView.CELL_RADIUS
     delta_x = sqrt(abs(GroundView.CELL_RADIUS * GroundView.CELL_RADIUS - delta_h * delta_h))
     self.__canvas.create_oval(arc_coord, fill='white')
     if percent >= 0.1:
         if percent >= 0.99:
             self.__canvas.create_oval(arc_coord, fill='green', outline='black' if is_good else 'red')
         else:
             r = (int(510 - 510 * percent) if percent >= 0.5 else 255) * 16 * 16 * 16 * 16
             g = int(255 if percent >= 0.5 else percent * 2 * 255) * 16 * 16
             color = '#%06x' % (r + g)
             tri_coord = x - delta_x, y - delta_h, \
                         x + delta_x, y - delta_h, \
                         x, y
             self.__canvas.create_arc(arc_coord, start=degrees(asin(2 * percent - 1)),
                                      extent=-180 - 2 * degrees(asin(2 * percent - 1)), fill=color,
                                      outline=color)
             self.__canvas.create_polygon(tri_coord, fill=color if percent >= 0.5 else 'white',
                                          outline=color if percent >= 0.5 else 'white',
                                          width=1 if percent > 0.5 else 2)
     elif percent <= -0.1:
         if percent <= - 0.99:
             self.__canvas.create_oval(arc_coord, fill='grey', outline='black' if is_good else 'red')
         else:
             tri_coord = x - delta_x, y + delta_h, \
                         x + delta_x, y + delta_h, \
                         x, y
             self.__canvas.create_arc(arc_coord, start=-degrees(asin(2 * -percent - 1)),
                                      extent=180 + 2 * degrees(asin(2 * -percent - 1)), fill='grey', outline='grey')
             self.__canvas.create_polygon(tri_coord, fill='grey' if -percent >= 0.5 else 'white',
                                          outline='grey' if -percent >= 0.5 else 'white',
                                          width=1 if abs(percent) > 0.5 else 3)
     self.__canvas.create_oval(arc_coord, outline='black' if is_good else 'red')

def solar_intensity(lat = 40.713, lng = -74.006, date_time = datetime.datetime.now()):
    # day_of_year = datetime.datetime.strptime(date_time, '%m/%d/%Y').timetuple().tm_yday
    day_of_year = date_time.timetuple().tm_yday
    declination_angle = degrees(asin(sin(radians(23.45)) * sin(2 * pi/365 * (day_of_year - 81))))
    # print 'declination_angle = ' + str(declination_angle)
    # ~~~equation of time~~~ to account for orbital eccentricity and axial wobble
    e_o_t = 9.87 * sin(2.0 * 2 * pi / 365 * (day_of_year - 81)) - 7.53 * cos(2 * pi / 365 * (day_of_year - 81)) - 1.5 * sin(2 * pi / 365 * (day_of_year - 81))
    gmt_diff = floor(lng / 15) # hour difference from GMT (+5 for EST)
    lstm = 15.0 * gmt_diff # local standard time meridian (edge of time zone)
    tc = 4.0 * (lng - lstm) + e_o_t # time correction factor
    local_solar_time = date_time + datetime.timedelta(0, tc/60.0)
    lst_dec = local_solar_time.hour + local_solar_time.minute / 60.0 + local_solar_time.second / 3600.0
    hour_angle = 15.0 * (lst_dec - 12.0)
    # print 'hour angle past solar noon = ' + str(hour_angle)
    # with a flat panel, the elevation angle is the incident angle
    elevation_angle = degrees(asin(sin(radians(declination_angle)) * sin(radians(lat)) + cos(radians(declination_angle)) * cos(radians(lat)) * cos(radians(hour_angle))))
    # print elevation_angle
    rad_from_vert = pi / 2 - radians(elevation_angle) # zenith angle in radians for airmass formula
    if rad_from_vert < pi / 2:
        air_mass = 1 / (cos(rad_from_vert) + 0.50572 * (96.07995 - degrees(rad_from_vert)) ** (-1.6364))
    else:
        air_mass = None
    e_0 = 1367 * (1 + 0.033 * cos(2 * pi * (day_of_year - 3) / 365))
    if air_mass:
        intensity = cos(rad_from_vert) * e_0 * 0.7 ** (air_mass ** 0.678)
        return intensity
    else:
        return 0

    def _determineBinAndBoutInFourAndFiveCirclesModes(self, hklNorm):
        """(Bin, Bout) = _determineBinAndBoutInFourAndFiveCirclesModes()"""
        BinModes = ('4cBin', '5cgBin', '5caBin')
        BoutModes = ('4cBout', '5cgBout', '5caBout')
        BeqModes = ('4cBeq', '5cgBeq', '5caBeq')
        azimuthModes = ('4cAzimuth')
        fixedBusingAndLeviWmodes = ('4cFixedw')

        # Calculate RHS of equation 20
        # RHS (1/K)(S^-1*U*B*H)_3 where H/K = hklNorm
        UB = self._getUBMatrix()
        [SIGMA, TAU] = createVliegsSurfaceTransformationMatrices(
            self._getSigma() * TORAD, self._getTau() * TORAD)
        #S = SIGMA * TAU
        S = TAU * SIGMA
        RHS = (S.I * UB * hklNorm)[2, 0]

        if self._getMode().name in BinModes:
            Bin = self._getParameter('betain')
            check(Bin != None, "The parameter betain must be set for mode %s" %
                  self._getMode().name)
            Bin = Bin * TORAD
            sinBout = RHS - sin(Bin)
            check(fabs(sinBout) <= 1, "Could not compute Bout")
            Bout = asin(sinBout)

        elif self._getMode().name in BoutModes:
            Bout = self._getParameter('betaout')
            check(Bout != None, "The parameter Bout must be set for mode %s" %
                  self._getMode().name)
            Bout = Bout * TORAD
            sinBin = RHS - sin(Bout)
            check(fabs(sinBin) <= 1, "Could not compute Bin")
            Bin = asin(sinBin)

        elif self._getMode().name in BeqModes:
            sinBeq = RHS / 2
            check(fabs(sinBeq) <= 1, "Could not compute Bin=Bout")
            Bin = Bout = asin(sinBeq)

        elif self._getMode().name in azimuthModes:
            azimuth = self._getParameter('azimuth')
            check(azimuth != None, "The parameter azimuth must be set for "
                  "mode %s" % self._getMode().name)
            del azimuth
            # TODO: codeit
            raise NotImplementedError()

        elif self._getMode().name in fixedBusingAndLeviWmodes:
            bandlomega = self._getParameter('blw')
            check(bandlomega != None, "The parameter abandlomega must be set "
                  "for mode %s" % self._getMode().name)
            del bandlomega
            # TODO: codeit
            raise NotImplementedError()
        else:
            raise RuntimeError("AngleCalculator does not know how to handle "
                               "mode %s" % self._getMode().name)

        return (Bin, Bout)

    def monitorArm(self, task):
        def clamp1(x): return min(1, max(-1, x))

        direction = self.ralph.get_pos() - self.robotarm.get_pos()
        direction.z += 1
        direction.normalize()

        # camera starts facing along x
        dec = math.asin(direction.x)
        cosdec = math.cos(dec) 
        if cosdec > 1e-05:
          ra = math.asin(clamp1(direction.z / cosdec))
          ra2 = math.acos(clamp1(direction.y / cosdec))
        else: ra = ra2 = math.pi/2#float('nan')
        #print(cosdec, direction)
        #print(dec, ra, ra2, cosdec)

        if direction.z < 0: 
          if ra2 < math.pi/2: ra2 = 0
          else: ra2 = math.pi

        self.jointForearm.setH(-dec * 180/math.pi)
        self.jointBase.setP(ra2 * 180/math.pi) 

        self.dirText.setText(str(direction)) 
        self.anglesText.setText(str((dec, ra, ra2, cosdec)))

        a = self.jointForearm.getH() / 90.0 * 300 + 512
        b = self.jointBase.getP()    / 90.0 * 300 + 212
        #print(a, b)
        baseID = 9
        servos.setAngle({baseID:int(round(b)), (baseID+1):int(round(a))})
        return task.again

    def _determineBinAndBoutInZaxisModes(self, Hw3OverK):
        """(Bin, Bout) = _determineBinAndBoutInZaxisModes(HwOverK)"""
        BinModes = ('6czBin')
        BoutModes = ('6czBout')
        BeqModes = ('6czBeq')

        if self._getMode().name in BinModes:
            Bin = self._getParameter('betain')
            check(Bin != None, "The parameter betain must be set for mode %s" %
                  self._getMode().name)
            Bin = Bin * TORAD
            # Equation 32a:
            Bout = asin(Hw3OverK - sin(Bin))

        elif self._getMode().name in BoutModes:
            Bout = self._getParameter('betaout')
            check(Bout != None, "The parameter Bout must be set for mode %s" %
                  self._getMode().name)
            Bout = Bout * TORAD
            # Equation 32b:
            Bin = asin(Hw3OverK - sin(Bout))

        elif self._getMode().name in BeqModes:
            # Equation 32c:
            Bin = Bout = asin(Hw3OverK / 2)

        return (Bin, Bout)

def ApproximateTransitTime(t = float('nan'), P = float('nan'), RS = float('nan'), a = float('nan')):
	"""Calculates the approximate transit time, using equation 2.6 in TE

	IN:
	t: transit duration, in s
	P: period of the planet, in s
	RS: radius of the host star, in m
	a: semi-major axis of the orbit, in m

	OUT:
	[t, P, RS, a]: list of all of the above
	"""

	#TODO: Extend this formula according to equation 3.4 in TE

	if(math.isnan(t)):
		t = P/math.pi * math.sin(RS/a)
	elif(math.isnan(P)):
		P = t * math.pi * math.sin(RS/a)
	elif(math.isnan(RS)):
		RS = a * math.asin(t*math.pi/P)
	elif(math.isnan(a)):
		a = RS / math.asin(t*math.pi/P)
	else:
		print ('TransitTime: WARNING: Nothing to solve')

	if(math.isnan(t) or math.isnan(P) or math.isnan(RS) or math.isnan(a)):
		print ('TransitTime: ERROR: Could not solve')

	return [t, P, RS, a]

    def get_spherical_coords( self ):
        '''returns the spherical coordinates from a cartesian coordinates

        using this formula:

            - http://www.math.montana.edu/frankw/ccp/multiworld/multipleIVP/spherical/body.htm#converting

        :rtype: ( radius, zenith, azimuth )
        '''
        eye = self.target.camera.eye - self.target.camera.center

        radius = math.sqrt( pow(eye.x,2) + pow(eye.y,2) + pow(eye.z,2) )
        s = math.sqrt( pow(eye.x,2) + pow(eye.y,2) )
        if s == 0:
            s = 0.000000001
        r = radius
        if r == 0:
            r = 0.000000001

        angle_z = math.acos( eye.z / float(r) )
        if eye.x < 0:
            angle_x = math.pi - math.asin( eye.y / s )
        else:
            angle_x = math.asin( eye.y / s )

        radius = radius / self.target.camera.get_z_eye()
        return radius, angle_z, angle_x