docs: 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):
  	"""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."""
  
  	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