docs: def calculate_reduced_matrix_elements_0(fine_states):
  	'''This function calculates the reduced matrix elments <N,L,J||T^1(r)||N',L',J'> given a list of fine states.'''
  	#We calculate the reduced matrix elements starting from the list of fine_states
  
  	#The factor composed of physical quantities.
  	factor=sqrt(3*c**3*me**2*e**2/(16*Pi*epsilon0*hbar**3))
  	#We read the database to obtain the Einstein A coefficients in Hz.
  	einsteinA=get_einstein_A_matrix(fine_states)
  	#We read the database to obtain the transition frequencies in Hz.
  	omega_fine=calculate_omega_matrix(fine_states)
  
  	reduced_matrix_elements=[[0.0 for jj in range(len(fine_states))] for ii in range(len(fine_states))]
  
  	for ii in range(len(fine_states)):
  		i=fine_states[ii]
  		for jj in range(ii):
  			j=fine_states[jj]
  			t=Transition(i,j)
  			einsteinAij=einsteinA[ii][jj]
  			omega0=omega_fine[ii][jj]
  
  			#The formula is valid only for i =/= j so that omega0 =/= 0.
  			#Because fine states are asumed to be ordered by their energies we can asume that
  			# i decays in j.
  			Ji=i.j; Jj=j.j
  
  			rij=sqrt((2.0*Ji+1)/(2*Jj+1))*sqrt(einsteinAij/omega0**3)
  			rij=factor*rij
  
  			reduced_matrix_elements[ii][jj]=rij
  			#We add the matrix elements on the other side of the diagonal.
  			reduced_matrix_elements[jj][ii]=rij*(-1)**(Jj-Ji)
  
  	return reduced_matrix_elements