*     *****************************
      FUNCTION Y(L,N,COSTHE,COSPHI)
*     *****************************
      INTEGER L,N
      REAL COSTHE,COSPHI,Y
*
*     This subroutine evaluates spherical harmonic functions
*
*     Y   (theta,phi),   where m=ABS(N) and s is sign of N.
*      lms
*
      INTEGER M
      REAL SINTHE
*
      M=ABS(N)
      SINTHE=SQRT(1-COSTHE**2)
*
*     First operate with l only - branch according to it.
*     Odd values of l are invalid for powder neutron diffraction.
*
      GO TO (900,20,900,40,900,60,900,80),L
      GO TO 900
*                             l = 2
   20 IF(M .EQ. 0) THEN
         YL=(3*COSTHE**2-1)/2
      ELSE IF(M .EQ. 1) THEN
         YL=2*COSTHE*SINTHE
      ELSE
         YL=SINTHE**2
      END IF
      GO TO 100
*                             l = 4
   40 IF(M .EQ. 0) THEN
         YL=(COSTHE**2*(35*COSTHE**2-30)+3)/8
      ELSE IF(M. EQ. 1) THEN
         YL=0.94695*(7*COSTHE**2-3)*COSTHE*SINTHE
      ELSE IF(M .EQ. 2) THEN
         YL=7*(7*COSTHE**2-1)*SINTHE**2/9
      ELSE IF(M .EQ. 3) THEN
         YL=3.079201436*COSTHE*SINTHE**3
      ELSE
         YL=SINTHE**4
      END IF
      GO TO 100
*                             l = 6
   60 IF(M .EQ. 0) THEN
         YL=(((231*COSTHE**2-315)*COSTHE**2+105)*COSTHE**2-5)/16
      ELSE IF(M .EQ. 1) THEN
         GO TO 900
      ELSE IF(M .EQ. 2) THEN
         YL=0.64549*((33*COSTHE**2-18)*COSTHE**2+1)*SINTHE**2
      ELSE IF(M .EQ. 3) THEN
         YL=1.41685*(11*COSTHE**2-3)*COSTHE*SINTHE**3
      ELSE IF(M .EQ. 4) THEN
         YL=0.81675*(11*COSTHE**2-1)*SINTHE**4
      ELSE IF(M .EQ. 5) THEN
         GO TO 900
      ELSE
         YL=SINTHE**6
      END IF
*                             l = 8
      GO TO 100
   80 IF(M .EQ. 0) THEN
         YL=((((6435*COSTHE**2-12012)*COSTHE**2+6930)*COSTHE**2-1260)*
     *          COSTHE**2+35)/128
      ELSE
         GO TO 900
      END IF
*
*     Now m, that is, phi-dependence is added.
*
  100 IF(N .GT. 0) THEN
         Y=YL*COS(M*ACOS(COSPHI))
      ELSE IF(N .LT. 0) THEN
         Y=YL*SIN(M*ACOS(COSPHI))
      ELSE
         Y=YL
      END IF
*
      RETURN
*
*     Error: L is illegal.
*
  900 WRITE(6,1000) L,M
      STOP
 1000 FORMAT(' In Y: Invalid indices; L =',I3,' M =',I3,'.')
      END
