      SUBROUTINE ABSRP(DIN,DOUT)
C
C ROUTINE TO PERFORM ANALYTICAL ABSORPTION FUNCTION OF
C   KATAYAMA (1986) ACTA CRYST A42, 19-23.
C
      INCLUDE 'ABSORB.INC'
      DIMENSION DIN(3),DOUT(3),PNMP(11),PNMS(11)
C
      ABSFCT=1.0
C NO ABSORPTION TO BE CALCULATED
      IF(NABSRN.EQ.0) RETURN
C ABSORPTION TO BE CALCULATED - GET ANGLES
      PHIP=ATAN2(-DIN(2),-DIN(1))
      PHIS=ATAN2(DOUT(2),DOUT(1))
      NP=0
      DO 120 N=1,NABSRN
C EVALUATE ASSOCIATED LEGENDRE FUNCTIONS
      CALL ASSLEG(PNMP,N,2*N+1,-DIN(3),IER)
      CALL ASSLEG(PNMS,N,2*N+1,DOUT(3),IER)
C SUM COS(M*PHI) TYPE TERMS
      DO 100 M=0,N
      NP=NP+1
      DABSCF(NP)=PNMP(N+M+1)*COS(M*PHIP)+PNMS(N+M+1)*COS(M*PHIS)
      ABSFCT=ABSFCT+ABSCF(NP)*DABSCF(NP)
100   CONTINUE
C SUM SIN(M*PHI) TYPE TERMS
      DO 110 M=1,N
      NP=NP+1
      DABSCF(NP)=PNMP(N+M+1)*SIN(M*PHIP)+PNMS(N+M+1)*SIN(M*PHIS)
      ABSFCT=ABSFCT+ABSCF(NP)*DABSCF(NP)
110   CONTINUE
120   CONTINUE
      RETURN
      END
      SUBROUTINE ASSLEG(PLM,L,N,X,IER)
C
C Evaluates the Associated Legendre functions Plm(x) for m=-1,1 using
C  relationships in Abramowitz and Stegun, Chapter 8.  Use is made of the
C  recurrence relation (8.5.3) to compute functions of order m=0 for
C  degree l=0 to L, and relation (8.5.1) for the positive orders m.
C  Relation (8.2.5) provides values for the functions of negative order.
C
C  ARGUMENTS : (* indicates an overwritten argument)
C
C  PLM*   Associated Legendre function Plm(x) contained in PLM(l+m+1)
C  L      Degree of function, l
C  N      Dimension of array PLM, equal to 2*l+1
C  X      Real argument x, mod(x).le.1.0
C  IER*   Error flag, if 0 - no error, if 1 - invalid argument.
C
C  Routine obtained from TCM Library, Cavendish Laboratory.
C    Author unknown.
C
      DIMENSION PLM(N)
      DATA EMACH/1.0E-5/
C
      DO 100 I=1,N
        PLM(I)=0.0
100   CONTINUE
      IF(ABS(ABS(X)-1.0) .LT. EMACH) GO TO 110
      IF(ABS(X)-1.0)120,110,180
110   PLM(L+1)=X**L
      RETURN
120   PLM(1)=1.0
      IF(L.EQ.0)RETURN
      PLM(2)=X
      IF(L.EQ.1)GO TO 140
      DO 130 I=2,L
        PLM(I+1)=((2*I-1)*X*PLM(I)-(I-1)*PLM(I-1))/I
130   CONTINUE
140   FAC=1.0/SQRT(1.0-X*X)
      DO 160 M=1,L
        LM=L+M
        PLM(LM+1)=FAC*(-(L-M+1)*X*PLM(LM)+(LM-1)*PLM(L))
        IF(M.EQ.L)GOTO 160
        NN=M+1
        DO 150 I=NN,L
          II=L-I+NN
          PLM(II)=FAC*(-(II-M)*X*PLM(II)+(II+M-2)*PLM(II-1))
150     CONTINUE
160   CONTINUE
      FAC=1.0
      DO 170 M=1,L
        FAC=FAC/((L+M)*(L-M+1))
        LM=L+1+M
        LN=L+1-M
        PLM(LN)=FAC*PLM(LM)
170   CONTINUE
      IER=0
180   IER=1
      RETURN
      END
