comment for the pl_leg function says that it returns 0 if n exceeds the...

comment for the pl_leg function says that it returns 0 if n exceeds the maximum Legendre order, but I don't think this was actually implemented, so I added it. I also changed the Fourier function to by default calculate A for the more common case of n > 0, and only calculate if n=0, instead of the reverse. Refs #8

examples/integration_example/integration_example.usr

@@ -535,41 +535,47 @@ c-----------------------
       return
       end
 c-----------------------------------------------------------------------
+! calculates Legendre polynomials using a recurrence relationship. If
+! n > the maximum Legendre order, the function returns 0.0.
       function pl_leg(x,n)
-!======================================
-!    calculates Legendre polynomials Pn(x)
-!    using the recurrence relation
-!    if n > 100 the function retuns 0.0
-!======================================
+      parameter (nl_max=100)
       real*8 pl,pl_leg
       real*8 x
       real*8 pln(0:n)
       integer n, k
+
       pln(0) = 1.0
       pln(1) = x
+
       if (n.le.1) then
         pl = pln(n)
-          else
-            do k=1,n-1
-              pln(k+1)=
-     &  ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
-            end do
-            pl = pln(n)
+      else if (n.le.nl_max) then
+        do k=1,n-1
+          pln(k+1) = ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
+        end do
+        pl = pln(n)
+      else
+        pl = 0.0
       end if
-      pl_leg=pl*sqrt(dble(2*n+1)/2.0)
-      return 
+
+      pl_leg = pl*sqrt(dble(2*n+1)/2.0)
+
+      return
       end
-!======================================
+c-----------------------------------------------------------------------
+! calculates Fourier polynomials Fn(x)
       function fl_four(x,n)
       real*8 fl_four,pi
       real*8 x, A
       integer n, k
+
       pi=4.0*atan(1.0)
-      A=1.0/sqrt(2*pi)
-      if (n.gt.0) A=1.0/sqrt(pi)  
-      fl_four=A*cos(n*x)    
-c       fl_four=1.0/sqrt(pi)
-      return 
+      A=1.0/sqrt(pi)
+
+      if (n.eq.0) A=1.0/sqrt(2*pi)
+      fl_four=A*cos(n*x)
+
+      return
       end
 c-----------------------------------------------------------------------
       subroutine heat_balance(fflux)