Calculate the two-theta array
970075
def twotheta(self) -> ErrorValue:
        """"""
        row, column = np.ogrid[0:self.shape[0], 0:self.shape[1]]
        rho = (((self.header.beamcentery - row) * self.header.pixelsizey) ** 2 +
               ((self.header.beamcenterx - column) * self.header.pixelsizex) ** 2) ** 0.5
        assert isinstance(self.header.pixelsizex, ErrorValue)
        assert isinstance(self.header.pixelsizey, ErrorValue)
        assert isinstance(rho, ErrorValue)
        return (rho / self.header.distance).arctan()
680423
function(x, y) {
        this.rho = Math.sqrt(x * x + y * y);
        this.theta = Math.atan2(y, x);
        if(this.theta < 0) this.theta += Math.PI * 2;
    }
543612
def N_to_Ntriangles(N):
    

    theta = np.array([np.pi/2*(k-0.5)/N for k in range(1, N+1)])
    phi = np.array([[np.pi*(l-0.5)/Mk for l in range(1, Mk+1)] for Mk in np.array(1 + 1.3*N*np.sin(theta), dtype=int)])
    Ntri = 2*np.array([len(p) for p in phi]).sum()

    return Ntri
762531
def gamma(phi1,phi2,theta1,theta2):
    
    
    if phi1 == phi2 and theta1 == theta2:
        gamma = 0
    else:
        gamma = atan( sin(theta2)*sin(phi2-phi1) / \
                      (cos(theta1)*sin(theta2)*cos(phi1-phi2) - \
                       sin(theta1)*cos(theta2)) )

    dummy_arg = (cos(gamma)*cos(theta1)*sin(theta2)*cos(phi1-phi2) + \
                 sin(gamma)*sin(theta2)*sin(phi2-phi1) - \
                 cos(gamma)*sin(theta1)*cos(theta2))
        
    if dummy_arg >= 0:
        return gamma
    else:
        return pi + gamma
825631
def rotation_matrix_from_axis_and_angle(u, theta):
    ''''''
    
    x = u[0]
    y = u[1]
    z = u[2]
    s = np.sin(theta)
    c = np.cos(theta)
    
    return np.array([[c + x**2 * (1 - c), x * y * (1 - c) - z * s, x * z * (1 - c) + y * s],
                     [y * x * (1 - c) + z * s, c + y**2 * (1 - c), y * z * (1 - c) - x * s ],
                     [z * x * (1 - c) - y * s, z * y * (1 - c) + x * s, c + z**2 * (1 - c) ]])
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
382105
def to_polar(x, y, theta_units="radians"):
    
    assert theta_units in [
        "radians",
        "degrees",
    ], "kwarg theta_units must specified in radians or degrees"

    theta = atan2(y, x)
    r = sqrt(x ** 2 + y ** 2)

    if theta_units == "degrees":
        theta = to_degrees(theta)

    return r, theta
76014
def lincon(self, x, theta=0.01):
        """"""
        if x[0] < 0:
            return np.NaN
        return theta * x[1] + x[0]
175389
def _tm(self, theta, phi, psi, dx, dy, dz):
        
        matrix = self.get_matrix(theta, phi, psi, dx, dy, dz)
        coord = matrix.dot(self.coord2)
        dist = coord - self.coord1

        d_i2 = (dist * dist).sum(axis=0)

        tm = -(1 / (1 + (d_i2 / self.d02)))

        return tm
949032
def _calc_theta(self):
        """"""
        if self.decaying_prior:
            n_sampled = np.clip(self.alpha_ + self.beta_, 1, np.inf)
            prior_weight = 1/n_sampled
            alpha = self.alpha_ + prior_weight * self.alpha_0
            beta = self.beta_ + prior_weight * self.beta_0
        else:
            alpha = self.alpha_ + self.alpha_0
            beta = self.beta_ + self.beta_0

        # Mean of Beta-distributed rv
        self.theta_ = alpha / (alpha + beta)

        # NEW: calculate theta assuming weak prior
        if self.store_wp:
            alpha = self.alpha_ + self._wp_weight * self.alpha_0
            beta = self.beta_ + self._wp_weight * self.beta_0
            self.theta_wp_ = alpha / (alpha + beta)
405450
def get_G(self, k, add_noise=True):
        
        G = np.ones((self.N - 1, self.N - 1))
        if (add_noise):
            noise = pi * 0.1 * np.random.rand(
                (self.N - 1) * (self.N - 1)).reshape((self.N - 1, self.N - 1))
        other_indices = np.delete(range(self.N), k)
        for idx, i in enumerate(other_indices):
            for jdx, j in enumerate(other_indices):
                if (add_noise and
                        i != j):  # do not add noise on diagonal elements.
                    thetak_ij = self.get_inner_angle(k,
                                                     (i, j)) + noise[idx, jdx]
                else:
                    thetak_ij = self.get_inner_angle(k, (i, j))
                G[idx, jdx] = cos(thetak_ij)
                G[jdx, idx] = cos(thetak_ij)
        return G
89639
def function(self, x, y, Rs, theta_Rs, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        e = min(abs(1. - q), 0.99)
        xt1 = (cos_phi*x_shift+sin_phi*y_shift)*np.sqrt(1 - e)
        xt2 = (-sin_phi*x_shift+cos_phi*y_shift)*np.sqrt(1 + e)
        R_ = np.sqrt(xt1**2 + xt2**2)
        rho0_input = self.nfw._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs)
        if Rs < 0.0000001:
            Rs = 0.0000001
        f_ = self.nfw.nfwPot(R_, Rs, rho0_input)
        return f_
1344565
private function _theta(PointInterface $a, PointInterface $b)
    {
        list($x1, $y1) = $a->getCoordinates();
        list($x2, $y2) = $b->getCoordinates();
        return atan2($y2-$y1, $x2 - $x1) + 2 * M_PI;
    }
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)]
    ])
1226088
def _g(self, z):
        
        return np.exp(np.multiply(-self.theta, z)) - 1
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
646314
def bbox(self):
        """"""
        # a(t) = radians(self.theta + self.delta*t)
        #      = (2*pi/360)*(self.theta + self.delta*t)
        # x'=0: ~~~~~~~~~
        # -rx*cos(phi)*sin(a(t)) = ry*sin(phi)*cos(a(t))
        # -(rx/ry)*cot(phi)*tan(a(t)) = 1
        # a(t) = arctan(-(ry/rx)tan(phi)) + pi*k === atan_x
        # y'=0: ~~~~~~~~~~
        # rx*sin(phi)*sin(a(t)) = ry*cos(phi)*cos(a(t))
        # (rx/ry)*tan(phi)*tan(a(t)) = 1
        # a(t) = arctan((ry/rx)*cot(phi))
        # atanres = arctan((ry/rx)*cot(phi)) === atan_y
        # ~~~~~~~~
        # (2*pi/360)*(self.theta + self.delta*t) = atanres + pi*k
        # Therfore, for both x' and y', we have...
        # t = ((atan_{x/y} + pi*k)*(360/(2*pi)) - self.theta)/self.delta
        # for all k s.t. 0 < t < 1
        from math import atan, tan

        if cos(self.phi) == 0:
            atan_x = pi/2
            atan_y = 0
        elif sin(self.phi) == 0:
            atan_x = 0
            atan_y = pi/2
        else:
            rx, ry = self.radius.real, self.radius.imag
            atan_x = atan(-(ry/rx)*tan(self.phi))
            atan_y = atan((ry/rx)/tan(self.phi))

        def angle_inv(ang, k):  # inverse of angle from Arc.derivative()
            return ((ang + pi*k)*(360/(2*pi)) - self.theta)/self.delta

        xtrema = [self.start.real, self.end.real]
        ytrema = [self.start.imag, self.end.imag]

        for k in range(-4, 5):
            tx = angle_inv(atan_x, k)
            ty = angle_inv(atan_y, k)
            if 0 <= tx <= 1:
                xtrema.append(self.point(tx).real)
            if 0 <= ty <= 1:
                ytrema.append(self.point(ty).imag)
        xmin = max(xtrema)
        return min(xtrema), max(xtrema), min(ytrema), max(ytrema)