This function returns code for Theta_i j as defined in the equation labeled Theta. in terms
	of detunings. It recieves indexes i,j starting from 1.
432953
def Theta(i,j,theta,omega_rescaled,omega_min,
			detunings,detuningsij,combinations,detuning_indices,
			Lij,i_d,I_nd,Nnd,
			states=None,verbose=1,other_the=None):
	""""""

	if i==j:
		return ''
	elif j>i:
		raise ValueError,'i should never be less than j.'

	Ne=len(omega_rescaled[0]); Nl=len(theta[0])
	if other_the!=None:
		the=other_the[:]
	else:
		the=[theta[j-1][l]-theta[i-1][l] for l in range(Nl)]

	###########################################################################################
	#This part is about finding detunings that reduce to -omega_ij
	if the==[0 for l in range(Nl)]:
		if omega_rescaled[i-1][j-1]==0:
			return ''
		#We need a combination of detunings that yields -omega_i,j.
		#According to equation labeled "detuning-exception1"
		#-omega_ij= delta^l_ik - delta^l_jk
		#so we seek a k that is lower than i and j
		#and an l such that l is in L_ik and in L_jk

		ii=min(i,j)
		band1=False
		band2=False
		for k in range(1,ii):
			for l in range(1,Nl+1):
				if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
					s   ='I found that -omega_'+str(i)+','+str(j)
					s+='= delta^'+str(l)+'_'+str(i)+','+str(k)
					s+= '-delta^'+str(l)+'_'+str(j)+','+str(k)
					s+='\n and also l='+str(l)
					s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
					s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'

					#print s
					band1=True
					break
			if band1:
				break

		if not band1:
			#If no such a k exists then we follow te equation labeled "detuning-exception2"
			#-omega_ij= delta^l_kj - delta^l_ki
			#to look for a k greater than i and j
			#and an l such that l is in L_ik and in L_jk
			for k in range(ii+1,Ne+1):
				for l in range(1,Nl+1):
					if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
						s   ='I found that -omega_'+str(i)+','+str(j)
						s+='= delta^'+str(l)+'_'+str(k)+','+str(j)
						s+='-delta^'+str(l)+'_'+str(k)+','+str(i)
						s+='\n and also l='+str(l)
						s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
						s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'
						#print s
						band2=True
						break

				if band2:
					break

		if band1:
			The=''
			#We need to find which detunings are delta^l_i,k and delta^l_j,k.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(i) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(j) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		elif band2:
			The=''
			#We need to find which detunings are delta^l_k,j and delta^l_k,i.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(j):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(i):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		else:
			if verbose>1: print 'WARNING: Optical frequencies will be used instead for -omega_',i,j,'=',omega_rescaled[i-1][j-1],'\n'
			return format_double(-omega_rescaled[i-1][j-1])
	###########################################################################################
	#This part is about finding detunings that reduce to (theta_j -theta_i -omega_ij)

	#We establish to which omegaij the detunings should reduce to
	omegaij=[0 for k in range(Nnd)]
	omegaij[I_nd(i)-1]=1;omegaij[I_nd(j)-1]=-1

	#We search for a combination of detunings that reduces to omegaij
	#That is a linear combination of detunings with coefficients the
	#that reduces to omegaij. Here
	detuning_failure=True
	for comb in combinations:
		suma=[sum([the[l]*detunings[l][comb[l]][ii] for l in range(Nl)])  for ii in range(Nnd)]
		if suma==omegaij:
			detuning_failure=False
			break

	#We stop if no combination reduces to the needed expression.
	if detuning_failure:

		if verbose>1: print 'We will see if it is possible to express Theta_'+str(i)+','+str(j)+'=theta_'+str(j)+'-theta_'+str(i)+'-omega_'+str(i)+','+str(j)
		if verbose>1: print 'in terms of other indices a,b such that omega_ab=omega_ij and transition i -> j is allowed by Lij.'

		band3=False
		for a in range(Ne):
			for b in range(a):
				if omega_rescaled[a][b]==omega_rescaled[i-1][j-1] and Lij[a][b]!=[]:
					band3=True
					break
			if band3: break
		if verbose>1: print band3,a,b,i,j

		a=a+1; b=b+1
		if band3 and (i!=a or j!=b):
			if verbose>1: print omega_rescaled[i-1][j-1],omega_rescaled[a-1][b-1]
			if verbose>1: print the
			if verbose>1: print 'This was possible for omega_'+str(a)+','+str(b)
			return Theta(a,b,theta,omega_rescaled,omega_min,
					detunings,detuningsij,combinations,detuning_indices
					,Lij,other_the=the,verbose=verbose,states=states)
		else:
			#verbose=2
			#print 111
			#qi=states[i-1].quantum_numbers
			#qj=states[j-1].quantum_numbers
			#print verbose
			if verbose>0: print 'WARNING: It was impossible to express laser frequencies',the
			if verbose>0: print 'and atomic transition -omega_',i,j,'=',-omega_rescaled[i-1][j-1]
			if verbose>0: print 'in terms of detunings from the transitions given by Lij'
			#for i in range(Ne):
			#	print i+1,Lij[i]
			#print
			if verbose>0: print 'I Will use optical frequencies instead.'
			The=''
			#We give the optical frequencies
			for l in range(Nl):
				a=the[l]
				if a==1:
					The+='+'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				elif a==-1:
					The+='-('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'
				elif a==0:
					The+=''
				elif a>0:
					The+='+'+str(a)+'*'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				else:
					The+=    str(a)+'*('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'

			#We substract omega_ij
			The+=format_double(-omega_rescaled[i-1][j-1])
			if verbose>1: print The
			if verbose>1: print
			return The

	#For each optical frequency in the, we write the corresponding detuning.
	#This way of assigining a global index ll to the detunings ammounts to
	#   ll=   number_of_previous_detunings
	#       + number_of_detuning_ordered_by_row_and_from_left_to_right_column
	The=''
	acum=0
	####################################################################
	#print 'comb',comb,'the',the,'i,j',i,j
	for l in range(Nl):
		if the[l]==1:
			The+='+detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==-1:
			The+='-detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==0:
			pass
		elif the[l]>0:
			The+='+'+str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		else:
			The+= str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		acum+=detuning_indices[l]

	if The[:1]=='+':
		The=The[1:]

	####################################################################
	#print 'The',The
	return The
164616
def getGradient(self,j):
        
        i = int(self.indicator['term'][j])
        r = int(self.indicator['row'][j])
        c = int(self.indicator['col'][j])
        rv = -np.kron(self.Fstar()[i][:,[r]],self.Astar()[i][[c],:])
        return rv
212898
def dot(self, vector, theta=None):
        
        if theta is not None:
            return (self.magnitude() * vector.magnitude() *
                    math.degrees(math.cos(theta)))
        return (reduce(lambda x, y: x + y,
                [x * vector.vector[i] for i, x in self.to_list()()]))
546782
function(jd) {
			// Jean Meeus algorithm, "Astronomical Algorithms", 1991
			var z = Math.floor(jd + 0.5);
			var a = Math.floor((z - 1867216.25) / 36524.25);
			a = z + 1 + a - Math.floor(a / 4);
			var b = a + 1524;
			var c = Math.floor((b - 122.1) / 365.25);
			var d = Math.floor(365.25 * c);
			var e = Math.floor((b - d) / 30.6001);
			var day = b - d - Math.floor(e * 30.6001);
			var month = e - (e > 13.5 ? 13 : 1);
			var year = c - (month > 2.5 ? 4716 : 4715);
			if (year <= 0) { year--; } // No year zero
			return this.newDate(year, month, day);
		}
471348
function _get_zs($i, $j = -1) {
        if ($j == -1) $j = $i;
        return substr($this->_txt, $this->_zmap[$i]['start'], $this->_zmap[$j]['end'] - $this->_zmap[$i]['start']);
    }
530875
def create_theta(self):
        
        import itertools
        from pylocus.basics_angles import from_0_to_pi
        theta = np.empty((self.M, ))
        corners = np.empty((self.M, 3))
        k = 0
        indices = np.arange(self.N)
        for triangle in itertools.combinations(indices, 3):
            for counter, idx in enumerate(triangle):
                corner = idx
                other = np.delete(triangle, counter)
                corners[k, :] = [corner, other[0], other[1]]
                theta[k] = self.get_inner_angle(corner, other)
                theta[k] = from_0_to_pi(theta[k])
                if DEBUG:
                    print(self.abs_angles[corner, other[0]],
                          self.abs_angles[corner, other[1]])
                    print('theta', corners[k, :], theta[k])
                k = k + 1
            inner_angle_sum = theta[k - 1] + theta[k - 2] + theta[k - 3]
            assert abs(inner_angle_sum - pi) < 1e-10, \
                'inner angle sum: {} {} {}'.format(
                    triangle, inner_angle_sum, (theta[k - 1], theta[k - 2], theta[k - 3]))
        self.theta = theta
        self.corners = corners
        return theta, corners
946904
function(i, a, d) {

      this.index = i;
      this.angle = a;
      this.distance = d;
   
      this.compare = function(p) {
         if (this.angle<p.angle)
            return -1;
         else if (this.angle>p.angle)
            return 1;
         else {
            if (this.distance<p.distance)
               return -1;
            else if (this.distance>p.distance)
               return 1;
         }
         return 0;
      }
   }
1217530
private final int m() {
        int n = 0;
        int i = 0;
        for (; cons(i); ++i) {
            if (i > j)
                return n;
        }
        i++;
        while(i <= j) {
            for (; !cons(i); ++i) {
                if (i > j)
                    return n;
            }
            i++;
            n++;

            if (i > j)
                return n;
            for (; cons(i); ++i) {
                if (i > j)
                    return n;
            }
            i++;
        }
        return n;
    }
449720
def unit_poly_verts(theta):
    
    x0, y0, r = [0.5] * 3
    verts = [(r*np.cos(t) + x0, r*np.sin(t) + y0) for t in theta]
    return verts
346562
def theta(v):
    
    v = np.atleast_2d(v)
    dir_z = v[:, 2]
    return theta_separg(dir_z)
600535
def Ry_matrix(theta):
    """"""
    return np.array([
        [np.cos(theta), 0, np.sin(theta)],
        [0, 1, 0],
        [-np.sin(theta), 0, np.cos(theta)]
    ])
285202
def ij_jlk_to_ilk(A, B):
    
    return A.dot(B.reshape(B.shape[0], -1)).reshape(A.shape[0], B.shape[1], B.shape[2])
275521
def sbesselh1(x, N):
    ""
    
    jn = sbesselj(x, N)
    yn = sbessely(x, N)

    return jn + 1j * yn
936340
private static boolean inCone(int i, int j, int n, int[] verts, int[] indices) {
        int pi = (indices[i] & 0x0fffffff) * 4;
        int pj = (indices[j] & 0x0fffffff) * 4;
        int pi1 = (indices[next(i, n)] & 0x0fffffff) * 4;
        int pin1 = (indices[prev(i, n)] & 0x0fffffff) * 4;
        // If P[i] is a convex vertex [ i+1 left or on (i-1,i) ].
        if (leftOn(verts, pin1, pi, pi1)) {
            return left(verts, pi, pj, pin1) && left(verts, pj, pi, pi1);
        }
        // Assume (i-1,i,i+1) not collinear.
        // else P[i] is reflex.
        return !(leftOn(verts, pi, pj, pi1) && leftOn(verts, pj, pi, pin1));
    }
1046587
def _theta_E_convert(self, theta_E, w_c, w_t):
        
        if not w_t > w_c:
            w_t, w_c = w_c, w_t
        s_scale_1 = w_c
        s_scale_2 = w_t
        f_x_1, f_y_1 = self.nie.derivatives(1, 0, theta_E=1, e1=0, e2=0, s_scale=s_scale_1)
        f_x_2, f_y_2 = self.nie.derivatives(1, 0, theta_E=1, e1=0, e2=0, s_scale=s_scale_2)
        f_x = f_x_1 - f_x_2
        theta_E_convert = theta_E / f_x
        return theta_E_convert, w_c, w_t
1983
def funTransSectX(theta, rho):
    
    return np.matrix([[np.cos(theta), rho * np.sin(theta)], [-np.sin(theta) / rho, np.cos(theta)]], dtype=np.double)
393189
def get_index_nested(x, i):
    
    for ind in range(len(x)):
        if i == x[ind]:
            return ind
    return -1