Have you written the code for the function 
? Is this function designed to work with non-scalar θ I.e. if you call it with:
? Is this function designed to work with non-scalar θ I.e. if you call it with:S(10, linspace(0, 2*pi, 10)) %f = 10, theta is a vector of 10 elements
do you get an error or do you get a vector of 10 elements as output?
The answer to this will change the way to code that summation.
