/* Auxiliary function */
g1(n, q, phi, chi) := [
Local(s);
s := q^2-n^2;
N(Cos(n*chi) * If(s=0,
1/2, /* Special case of s=0: avoid division by 0 */
Sin(Sqrt(s)*phi)/Sin(2*Sqrt(s)*phi) /* now s != 0 */
/* note that Sqrt(s) may be imaginary here */
)
);
];
/* Main function */
g(q, phi, chi) := [
Local(M, n);
M := 16; /* Exp(-M) will be the precision */
/* Use N() to always force evaluation */
N(1/2*Sin(q*phi)/Sin(2*q*phi)) +
/* Estimate the necessary number of terms in the series */
Sum(n, 1, N(1+Sqrt(q^2+M^2/phi^2)), g1(n, q, phi, chi)) ;
];
/* Parameters */
q:=3.5;
phi:=2;
/* Make a function for plotting: it must have only one argument */
f(x) := g(q, phi, x);
/* Plot from 0 to 2*Pi with 80 points */
GnuPlot(0, N(2*Pi), 80, f(x));
|