docs: def write_equations_code(path, name, laser, omega, gamma, r, Lij, states=None,
                           excluded_mu=[], verbose=1):
      r"""Write code for the equations."""
      Ne = len(omega[0])
      Nl = len(laser)
      N_excluded_mu = len(excluded_mu)
  
      if states is None: states = range(1, Ne+1)
  
      omega_rescaled = omega[:]
  
      # We determine whether it is possible to eliminate explicit time-dependance
      theta = find_phase_transformation(Ne, Nl, r, Lij)
  
      # We construct the correspondence i <-> I between degenerate and
      # non-degenerate indices.
      i_d, I_nd, Nnd = calculate_iI_correspondence(omega)
      # We find the non-degenerate detunings
      detunings, detuningsij = laser_detunings(Lij, Nl, i_d, I_nd, Nnd)
      aux = find_omega_min(omega_rescaled, Nl, detuningsij, i_d, I_nd)
      omega_min, omega_min_indices = aux
  
      detuning_indices = [len(detunings[i]) for i in range(Nl)]
      Nd = sum([len(detunings[l]) for l in range(Nl)])
      combinations = detuning_combinations(detuning_indices)
  
      ##########################################################################
  
      # We add the code to caculate the detunings for each laser.
      code0 = ''
      code0 += '	!We calculate the detunings.\n'
      # We find the minimal frequency corresponding to each laser.
      code0 += '	!The list of detunings has the following meaning:\n'
      conta = 0
      for ll in range(Nl):
          for kk in range(len(detuningsij[ll])):
              conta += 1
              ii, jj = detuningsij[ll][kk]
  
              if states is None:
                  code0 += '	!detuning('+str(conta)+')= delta^'+str(ll+1)+'_'
                  code0 += str(i_d(ii+1))+','+str(i_d(jj+1))+'\n'
              else:
                  state_i = str(states[i_d(ii+1)-1])[5:]
                  state_j = str(states[i_d(jj+1)-1])[5:]
                  code0 += '	!detuning('+str(conta)+')= delta^'+str(ll+1)+'_'
                  code0 += state_i+','+state_j+'\n'
  
      det_index = 1
      for l in range(Nl):
          # omega0 = omega_min[l]
          i_min, j_min = omega_min_indices[l]
          for ii, jj in detuningsij[l]:
              code0 += '	detuning('+str(det_index)+')='
              # code0+=format_double(omega0-omega_rescaled[ii][jj])+'+(detuning_knob('+str(l+1)+'))\n'
              code0 += 'detuning_knob('+str(l+1)+') '
              code0 += '-('+format_double(
                  omega_rescaled[i_d(ii+1)-1][i_d(i_min+1)-1])+')'
              code0 += '-('+format_double(
                  omega_rescaled[i_d(j_min+1)-1][i_d(jj+1)-1])+')\n'
              det_index += 1
      code0 += '\n'
  
      ###########################################################################
      # We add here the code for the equations
      ###########################################################################
      code = ''
      ###########################################################################
      # We need to check that the resulting matrix A doesn't have any
      # row or column made of zeroes. If there is such a row or column
      # it means that for that particular i(mu),j(mu) or i(nu),j(nu)
      # neither the hamiltonian part of the correspondin equation nor the
      # phase transformation part nor the decay part were nonzero
      # (at least for the given light polarizations).
      # In this cases it is necessary to re-calculate all of the problem
      # without the corresponding mu components of A.
      # Also the right hand side will be checked for non zero elements.
      # Now we start the check for these exceptions.
  
      row_check = [False for mu in range(Ne**2-1-N_excluded_mu)]
      col_check = [False for nu in range(Ne**2-1-N_excluded_mu)]
      rhs_check = [False for nu in range(Ne**2-1-N_excluded_mu)]
      ##########################################################################
      # We give the code to calculate the independent vector.
      code += '	!We calculate the independent vector.\n'
      for i in range(2, Ne+1-N_excluded_mu):
          for s in [1, -1]:
              nu = Mu(i, 1, s, Ne, excluded_mu)
              rhs_check[nu-1] = True
  
              for l in Lij[i-1][0]:
                  dp = dot_product(laser[l-1], 1, r, i, 1)
                  dp = s*part(dp, -s)
                  if dp != 0:
                      code += '	B('+str(nu)+',1)=B('+str(nu)+',1) +E0('+str(l)
                      code += ')*('+format_double(dp)+')\n'
  
      code += '\n'
      code += '	B=B/2.0d0\n\n'  # +str(1/sqrt(2.0))+'d0\n\n'
  
      ###########################################################################
      # We give the code to calculate the equations for populations.
      code += '	!We calculate the equations for populations.\n'
      for i in range(2, Ne+1):
          mu = Mu(i, i, 1, Ne, excluded_mu)
          for k in range(1, Ne+1):
              if k < i:
                  for l in Lij[k-1][i-1]:
                      dp1 = dot_product(laser[l-1], -1, r, k, i)
                      dp2 = dot_product(laser[l-1], 1, r, i, k)
  
                      real_coef = part(dp1-dp2, -1)
                      imag_coef = part(dp1+dp2, 1)
  
                      if real_coef != 0:
                          nu = Mu(i, k, 1, Ne, excluded_mu)
                          code += '	A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                          code += str(nu)+')'
                          code += '+E0('+str(l)+')*('+format_double(real_coef)
                          code += ')\n'
                          row_check[mu-1] = True; col_check[nu-1] = True
  
                      if imag_coef != 0:
                          nu = Mu(i, k, -1, Ne, excluded_mu)
                          code += '	A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                          code += str(nu)+')'
                          code += '+E0('+str(l)+')*('+format_double(imag_coef)
                          code += ')\n'
                          row_check[mu-1] = True; col_check[nu-1] = True
  
              if k > i:
                  for l in Lij[k-1][i-1]:
                      dp1 = dot_product(laser[l-1], -1, r, i, k)
                      dp2 = dot_product(laser[l-1], 1, r, k, i)
  
                      real_coef = -part(dp1-dp2, -1)
                      imag_coef = -part(dp1+dp2, 1)
  
                      if real_coef != 0:
                          nu = Mu(k, i, 1, Ne, excluded_mu)
                          code += '	A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                          code += str(nu)+')'
                          code += '+E0('+str(l)+')*('+format_double(real_coef)
                          code += ')\n'
                          row_check[mu-1] = True; col_check[nu-1] = True
  
                      if imag_coef != 0:
                          nu = Mu(k, i, -1, Ne, excluded_mu)
                          code += '	A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                          code += str(nu)+')'
                          code += '+E0('+str(l)+')*('+format_double(imag_coef)
                          code += ')\n'
                          row_check[mu-1] = True; col_check[nu-1] = True
  
      code += '\n'
      ###########################################################################
      # We give the code to calculate the equations for coherences
      # given in equations with label "stationary-coherences"
  
      code += '	!The code to calculate the equations for coherences.\n'
  
      for i in range(2, Ne+1):
          for j in range(1, i):
              for s in [1, -1]:
                  mu = Mu(i, j, s, Ne, excluded_mu)
                  for k in range(1, Ne+1):
                      for l in Lij[k-1][j-1]:
                          if k < i:
                              if k < j:
                                  dp = s*dot_product(laser[l-1], -1, r, k, j)
                                  dp1 = part(dp, -s)
                                  dp2 = part(s*dp, +s)
                                  nu = Mu(i, k, 1, Ne, excluded_mu)
                                  if dp1 != 0:
                                      code += '	A('+str(mu)+','+str(nu)+')=A('
                                      code += str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp1)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                                  nu = Mu(i, k, -1, Ne, excluded_mu)
                                  if dp2 != 0:
                                      code += '	A('+str(mu)+','+str(nu)+')=A('
                                      code += str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp2)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                              elif k > j:
                                  dp = s*dot_product(laser[l-1], +1, r, k, j)
                                  dp1 = part(dp, -s)
                                  dp2 = part(s*dp, s)
                                  nu = Mu(i, k, 1, Ne, excluded_mu)
                                  if dp1 != 0:
                                      code += '	A('+str(mu)+','+str(nu)+')=A('
                                      code += str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp1)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                                  nu = Mu(i, k, -1, Ne, excluded_mu)
                                  if dp2 != 0:
                                      code += '	A('+str(mu)+','+str(nu)+')=A('
                                      code += str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp2)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                          elif k > i:
                              dp = s*dot_product(laser[l-1], 1, r, k, j)
                              dp1 = part(dp, -s)
                              dp2 = part(-s*dp, s)
                              nu = Mu(k, i, 1, Ne, excluded_mu)
                              if dp1 != 0:
                                  code += '	A('+str(mu)+','+str(nu)
                                  code += ')=A('+str(mu)+','+str(nu)+')'
                                  code += '+E0('+str(l)+')*('
                                  code += format_double(dp1)+')\n'
                                  row_check[mu-1] = True; col_check[nu-1] = True
                              nu = Mu(k, i, -1, Ne, excluded_mu)
                              if dp2 != 0:
                                  code += '	A('+str(mu)+','+str(nu)+')=A('
                                  code += str(mu)+','+str(nu)+')'
                                  code += '+E0('+str(l)+')*('+format_double(dp2)
                                  code += ')\n'
                                  row_check[mu-1] = True; col_check[nu-1] = True
                      for l in Lij[i-1][k-1]:
                          if k > j:
                              if k > i:
                                  dp = -s*dot_product(laser[l-1], -1, r, i, k)
                                  dp1 = part(dp, -s)
                                  dp2 = part(s*dp, s)
                                  nu = Mu(k, j, 1, Ne, excluded_mu)
                                  if dp1 != 0:
                                      code += '	A('+str(mu)+','+str(nu)
                                      code += ')=A('+str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp1)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                                  nu = Mu(k, j, -1, Ne, excluded_mu)
                                  if dp2 != 0:
                                      code += '	A('+str(mu)+','+str(nu)
                                      code += ')=A('+str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp2)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                              elif k < i:
                                  dp = -s*dot_product(laser[l-1], 1, r, i, k)
                                  dp1 = part(dp, -s)
                                  dp2 = part(s*dp, s)
                                  nu = Mu(k, j, 1, Ne, excluded_mu)
                                  if dp1 != 0:
                                      code += '	A('+str(mu)+','+str(nu)
                                      code += ')=A('+str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp1)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                                  nu = Mu(k, j, -1, Ne, excluded_mu)
                                  if dp2 != 0:
                                      code += '	A('+str(mu)+','+str(nu)
                                      code += ')=A('+str(mu)+','+str(nu)+')'
                                      code += '+E0('+str(l)+')*('
                                      code += format_double(dp2)+')\n'
                                      row_check[mu-1] = True
                                      col_check[nu-1] = True
                          elif k < j:
                              dp = -s*dot_product(laser[l-1], 1, r, i, k)
                              dp1 = part(dp, -s)
                              dp2 = part(-s*dp, s)
                              nu = Mu(j, k, 1, Ne, excluded_mu)
                              if dp1 != 0:
                                  code += '	A('+str(mu)+','+str(nu)+')=A('
                                  code += str(mu)+','+str(nu)+')'
                                  code += '+E0('+str(l)+')*('
                                  code += format_double(dp1)+')\n'
                                  row_check[mu-1] = True; col_check[nu-1] = True
                              nu = Mu(j, k, -1, Ne, excluded_mu)
                              if dp2 != 0:
                                  code += '	A('+str(mu)+','+str(nu)+')=A('
                                  code += str(mu)+','+str(nu)+')'
                                  code += '+E0('+str(l)+')*('+format_double(dp2)
                                  code += ')\n'
                                  row_check[mu-1] = True; col_check[nu-1] = True
                  for l in Lij[i-1][j-1]:
                      dp = s*part(dot_product(laser[l-1], 1, r, i, j), -s)
                      nu = Mu(i, i, 1, Ne, excluded_mu)
                      if dp != 0:
                          code += '	A('+str(mu)+','+str(nu)+')=A('
                          code += str(mu)+','+str(nu)+')'
                          code += '+E0('+str(l)+')*('+format_double(dp)+')\n'
                          row_check[mu-1] = True; col_check[nu-1] = True
                          nu = Mu(j, j, 1, Ne, excluded_mu)
                          if nu == 0:
                              for n in range(1, Ne):
                                  code += '	A('+str(mu)+','+str(n)+')=A('
                                  code += str(mu)+','+str(n)+')'
                                  code += '+E0('+str(l)+')*('+format_double(dp)
                                  code += ')\n'
                                  row_check[mu-1] = True; col_check[n-1] = True
                          else:
                              code += '	A('+str(mu)+','+str(nu)+')=A('
                              code += str(mu)+','+str(nu)+')'
                              code += '+E0('+str(l)+')*('+format_double(-dp)
                              code += ')\n'
                              row_check[mu-1] = True; col_check[nu-1] = True
  
      code += '\n'
      code += '	A=A/2.0d0\n\n'  # +str(1/sqrt(2.0))+'d0\n\n'
      ###########################################################################
      # We add the terms associated with the phase transformation.
      code += '	!We calculate the terms associated with the phase '
      code += 'transformation.\n'
  
      # for i in range(2,10):
      for i in range(2, Ne+1):
          for j in range(1, i):
              extra = Theta(i, j, theta, omega_rescaled, omega_min, detunings,
                            detuningsij, combinations, detuning_indices, Lij,
                            i_d, I_nd, Nnd, verbose=verbose, states=states)
  
              if extra != '':
                  for s in [1, -1]:
                      mu = Mu(i, j, s, Ne, excluded_mu)
                      nu = Mu(i, j, -s, Ne, excluded_mu)
  
                      code += '    A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                      code += str(nu)+')'
                      if s == 1:
                          code += '-('+str(extra)+')\n'
                      elif s == -1:
                          code += '+('+str(extra)+')\n'
  
                      row_check[mu-1] = True; col_check[nu-1] = True
  
      code += '\n'
      ###########################################################################
      # We add the terms associated with spontaneous decay.
      code += '	!We calculate the terms associated with spontaneous decay.\n'
      # First for populations.
      for i in range(2, Ne+1):
          mu = Mu(i, i, 1, Ne, excluded_mu)
          for k in range(1, Ne+1):
              gams = 0
              if k < i:
                  gams += gamma[i-1][k-1]
              elif k > i:
                  nu = Mu(k, k, 1, Ne, excluded_mu)
                  ga = gamma[i-1][k-1]
                  if ga != 0:
                      code += '    A('+str(mu)+','+str(nu)+')=A('+str(mu)+','
                      code += str(nu)+')'
                      code += '-('+format_double(ga)+')\n'
                      row_check[mu-1] = True; col_check[nu-1] = True
              if gams != 0:
                  code += '    A('+str(mu)+','+str(mu)+')=A('+str(mu)+','
                  code += str(mu)+')'
                  code += '-('+format_double(gams)+')\n'
                  row_check[mu-1] = True; col_check[mu-1] = True
  
      # And now for coherences
      for i in range(1, Ne+1):
          for j in range(1, i):
              gams = gamma[i-1][j-1]/2
              if gams != 0:
                  for a in range(i+1, Ne+1):
                      mu = Mu(a, i, 1, Ne, excluded_mu)
                      code += '    A('+str(mu)+','+str(mu)+')=A('+str(mu)+','
                      code += str(mu)+')'
                      code += '-('+format_double(gams)+')\n'
                      row_check[mu-1] = True; col_check[mu-1] = True
                      mu = Mu(a, i, -1, Ne, excluded_mu)
                      code += '    A('+str(mu)+','+str(mu)+')=A('+str(mu)+','
                      code += str(mu)+')'
                      code += '-('+format_double(gams)+')\n'
                      row_check[mu-1] = True; col_check[mu-1] = True
  
                  for b in range(1, i):
                      mu = Mu(i, b, 1, Ne, excluded_mu)
                      code += '    A('+str(mu)+','+str(mu)+')=A('+str(mu)+','
                      code += str(mu)+')'
                      code += '-('+format_double(gams)+')\n'
                      row_check[mu-1] = True; col_check[mu-1] = True
                      mu = Mu(i, b, -1, Ne, excluded_mu)
                      code += '    A('+str(mu)+','+str(mu)+')=A('+str(mu)+','
                      code += str(mu)+')'
                      code += '-('+format_double(gams)+')\n'
                      row_check[mu-1] = True; col_check[mu-1] = True
  
      code += '\n'
      code = code0+code
      Nd = sum([len(detunings[l]) for l in range(Nl)])
      res = [code, Nd, row_check, col_check, rhs_check, Ne, N_excluded_mu,
             states, omega_min, detuningsij, omega_rescaled]
      return res