Writing a double sum in Python

Question

I am new to StackOverflow, and I am extremely new to Python.

My problem is this... I am needing to write a double-sum, as follows:

The motivation is that this is the angular correction to the gravitational potential used for the geoid.

1. I am having difficulty writing the sums. And please, before you say "Go to such-and-such a resource," or get impatient with me, this is the first time I have ever done coding/programming/whatever this is.

2. Is this a good place to use a "for" loop?

3. I have data for the two indices (n,m) and for the coefficients c_{nm} and s_{nm} in a .txt file. Each of those items is a column. When I say usecols, do I number them 0 through 3, or 1 through 4?

---

(the equation above)

\begin{equation}
V(r, \phi, \lambda) = \sum_{n=2}^{360}\left(\frac{a}{r}\right)^{n}\sum_{m=0}^{n}\left[c_{nm}*\cos{(m\lambda)} + s_{nm}*\sin{(m\lambda)}\right]*\sqrt{\frac{(n-m)!}{(n+m)!}(2n + 1)(2 - \delta_{m0})}P_{nm}(\sin{\lambda})
\end{equation}

Answer 1

(2) Yes, a "for" loop is fine. As @jpmc26 notes, a generator expression is a good alternative to a "for" loop. IMO, you'll want to use numpy if efficiency is important to you.

(3) As @askewchan notes, "usecols" refers to an argument of genfromtxt; as specified in that documentation, column indexes start at 0, so you'll want to use 0 to 3.

A naive implementation might be okay since the larger factorial is the denominator, but I wouldn't be surprised if you run into numerical issues. Here's something to get you started. Note that you'll need to define P() and a. I don't understand how "0 through 3" relates to c and s since their indexes range much further. I'm going to assume that each (and delta) has its own file of values.

import math
import numpy
c = numpy.getfromtxt("the_c_file.txt")
s = numpy.getfromtxt("the_s_file.txt")
delta = numpy.getfromtxt("the_delta_file.txt")
def V(r, phi, lam):
  ret = 0
  for n in xrange(2, 361):
    for m in xrange(0, n + 1):
      inner = c[n,m]*math.cos(m*lam) + s[n,m]*math.sin(m*lam)
      inner *= math.sqrt(math.factorial(n-m)/math.factorial(n+m)*(2*n+1)*(2-delta[m,0]))
      inner *= P(n, m, math.sin(lam))
    ret += math.pow(a/r, n) * inner
  return ret

Make sure to write unittests to check the math. Note that "lambda" is a reserved word.