How to use Sum(...)?

rowlesmr

I'm writing up a preferred orientation correction, and I need to do a sum.

The techref (p. 16) says I can use:

Sum(returns summation_eqn, initializer, conditional_test, increment_eqn)

Taking guidance from the Multiplicities_Sum macro, I did this:

(1/N)*Sum( PO_P(r,alpha,Delta,(j+(1/2))*Pi/N), j=0,j<=N-1,j=j+1)

but I get an error: Uninitialized_Variable in equation: j

If I do "prm j 0", Topas just hangs
If I do "prm !j 0", Topas catches an exception
If I do "...Pi/N), prm j=0,j<=...", Topas catches an exception

I can get around it by explicitly expanding the sum, but that isn't very elegant.

Is there a way to do a Sum using any in-built functions?

alancoelho

Sum can only use the internal iterator Mi as in the Multiplicities_Sum macro; ie.

Sum( sum_eqn, Mi = 0, Mi < M, Mi = Mi+1) / M

where sum_eqn can be an equation that is a function of H, K, L, M, D_spacing, 2Th. You will find that explicitly expanding the sum will in fact be quicker in this case as the simplification routines will have a chance to see the whole summation. Note, TOPAS will load 1000,000 equations and simplify them in little time; don't be afraid to use large expressions.

rowlesmr

Thanks for that Alan.

I'll just hard-code in N=32.

.

Further to my looking at preferred orientation corrections, I was looking at implementing Legendre polynomials to put in with the spherical harmonics correction, and I'm getting answers that don't match up.

I tried implementing the recursive relationship for Jacobi polynomials, and then the same relationship for Legendre polynomials, but for n>=3 and 4, respectively, the results diverge from what they should be.

This is what I did:

macro & JP_0 { 1 }
macro & JP_1 { 0.5*a - 0.5*b + (1+ 0.5*a + 0.5*b)*x }
macro & JP_n { ((2*n+a+b-1)*((2*n+a+b)*(2*n+a+b-2)*x+a^2-b^2)*Jacobi_polynomial(n-1,a,b,x) - 2*(n+a-1)*(n+b-1)*(2*n+a+b)*Jacobi_polynomial(n-2,a,b,x))/(2*n*(n+a+b)*(2*n+a+b-2)) }

fn Jacobi_polynomial(n,a,b,x) = If(n==0, JP_0, If(n==1, JP_1, JP_n));


macro & LP_0 { 1 }
macro & LP_1 { x }
macro & LP_n { (x*Legendre_polynomial_fn(n-1,x)*(2*n-1) - Legendre_polynomial_fn(n-2,x)*(n-1))/n }

fn Legendre_polynomial_fn(n,x) = If(n==0, LP_0, If(n==1, LP_1, LP_n));


prm = Jacobi_polynomial(0,0,0,0.1); :  1.00000 '      1
prm = Jacobi_polynomial(1,0,0,0.1); :  0.10000 '     0.1
prm = Jacobi_polynomial(2,0,0,0.1); : -0.48500 '-0.4850000000
prm = Jacobi_polynomial(3,0,0,0.1); : -0.08283 '-0.1475000000
prm = Jacobi_polynomial(4,0,0,0.1); :  0.66538 ' 0.3379375000
prm = Jacobi_polynomial(5,0,0,0.1); :  0.13262 ' 0.1788287500
prm = Jacobi_polynomial(6,0,0,0.1); : -1.56397 '-0.2488293125
prm = Jacobi_polynomial(7,0,0,0.1); : -0.31738 '-0.1994929438
prm = Jacobi_polynomial(8,0,0,0.1); :  5.23586 ' 0.1803207215

prm = Legendre_polynomial_fn(0,0.1); :  1.00000 '      1
prm = Legendre_polynomial_fn(1,0.1); :  0.10000 '     0.1
prm = Legendre_polynomial_fn(2,0.1); : -0.48500 '-0.4850000000
prm = Legendre_polynomial_fn(3,0.1); : -0.14750 '-0.1475000000
prm = Legendre_polynomial_fn(4,0.1); :  0.67588 ' 0.3379375000
prm = Legendre_polynomial_fn(5,0.1); :  0.39943 ' 0.1788287500
prm = Legendre_polynomial_fn(6,0.1); : -1.47000 '-0.2488293125
prm = Legendre_polynomial_fn(7,0.1); : -1.43586 '-0.1994929438
prm = Legendre_polynomial_fn(8,0.1); :  4.06811 ' 0.1803207215

The commented numbers after the equations are the correct result according to Maple. I got the recursive relationship from the Maple help files.

I've tested the same code (literally copy/paste the macro contents) in Julia, and I get the same results as Maple.

alancoelho

There could be a bug here in regards to recursive functions; will get back to you on this hopefully in the next week.