C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
      SUBROUTINE FNNMAT
      IMPLICIT REAL*8 (A-H,O-Z)
C      PARAMETER(MAXLP=19894)
      INCLUDE 'prmfit.f'
      CHARACTER*4 NCCAR
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/SCL/NATOM,NKAT,IPR,ISO,E0,DE,NNEE,NELEC,ENEND,AM
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT,NBAND,NB1,NB2
      COMMON/LCM/COMP(MAXEIG,MAXKPT,MAXCMP),NCOMP,
     $           NKCOMP(MAXCMP),NCCAR(MAXCMP)
      COMMON/LATVC/LATOM(3,MAXLP),ALM(MAXLP),MP(MAXLP),NLAT,NSHEL,ICD
	  COMMON/FITB/B(MAXKPT,MAXEIG,2),FF(MAXKPT,MAXKPT),JPMAT(MAXKPT)
     &   ,NFSHEL
	  DIMENSION F(MAXKPT,MAXKPT),C(MAXKPT)
	  DIMENSION FT(MAXKPT),BS(MAXKPT),BB(MAXKPT)
	  DIMENSION IPP(MAXKPT)
	  TPAI=DATAN(1.0D0)*8.0D0
	  WRITE(6,*) NKPINT,NSHEL,ISO
      NFSHEL=NKPINT
	  IF(ISO.EQ.2) NFSHEL=NKPINT/2
	  WRITE(6,*) ' NFSHEL=',NFSHEL,'NKPINT=',NKPINT
	  IP=1
	  JP=0
	  DO 12 J=1,NSHEL
C
      IF(IPR.GE.5) write(6,*) ' J=',J,IP,MP(J),' JP=',JP
      K=0
	  DO 11 KK=1,NKPINT
	  IF(IUDM(KK).NE.1) GO TO 11
	  K=K+1
C	  WRITE(6,*) (KXM(IK,KK),IK=1,3)

	  FT(K)=0.0D0
	  DO 13 I=IP,MP(J) 
	  RK=LATOM(1,I)*KXM(1,KK)+LATOM(2,I)*KXM(2,KK)
     & +LATOM(3,I)*KXM(3,KK)
	  RK=(RK*TPAI)/(ICD*ICM(KK))
	  FT(K)=FT(K)+COS(RK)
C	  WRITE(6,*) K,I,COS(RK)
   13 CONTINUE
      FT(K)=FT(K)/(MP(J)-IP+1)
C	  WRITE(6,*) ' FT(K)',K,J,FT(K)
   11 CONTINUE
C      WRITE(6,*) ' J=',J
C      WRITE(6,602) (FT(K),K=1,NFSHEL)
      IF(JP.EQ.0) GO TO 20
	  DO 15 JJ=1,JP
	     DO 16 K=1,NFSHEL
	        IF(DABS(FT(K)-F(K,JJ)).GT.1.0D-5) GO TO 15
   16    CONTINUE
C         WRITE(6,*) ' COIN ',JJ,J,JP,NFSHEL
      GO TO 19
   15 CONTINUE
   20 JP=JP+1
      JPMAT(JP)=J
      IF(IPR.GE.4) WRITE(6,*) ' BASE ',JP,J,MP(J),NFSHEL
      DO 17 K=1,NFSHEL
	     F(K,JP)=FT(K)
		 FF(K,JP)=FT(K)
   17 CONTINUE
      IF(JP.EQ.NFSHEL) GO TO 18   
   19 IP=MP(J)+1
   12 CONTINUE
	     WRITE(6,*) ' THE NUMBER OF BASES ARE TOO SMALL'
		 STOP ' STOP in FNNMAT'
   18 CONTINUE
C      DO 14 K=1,NFSHEL
C	  WRITE(6,*) (KXM(IK,K),IK=1,3),ICM(K)
C	  WRITE(6,602) (F(K,J),J=1,NFSHEL)
C  602 FORMAT(10F8.4)
C	  WRITE(9,602) (F(K,J),J=1,NFSHEL)	  
C   14 CONTINUE
      CALL MLUDIS(F,MAXKPT,NFSHEL,IPP,IER)
	  NSPIN=1
	  IF(ISO.EQ.2) NSPIN=2
	  DO 26 ISPIN=1,NSPIN
	  IB=0
	  DO 24 JB=NB1,NB2
	  IB=IB+1
	  I=0
	  DO 21 II=1,NKPINT
	     IF(IUDM(II).NE.ISPIN) GO TO 21
		 I=I+1
	     BB(I)=EIGM(JB,II)
		 BS(I)=BB(I)
   21 CONTINUE
       WRITE(6,*) ' INPUT ENERGY JB=',JB,' ISPIN=',ISPIN
       WRITE(6,610) (BB(K),K=1,NFSHEL)
	  CALL MLULEQ(F,MAXKPT,NFSHEL,BB,IPP,IERR)
      WRITE(6,*) ' COEFFICIENT JB=',JB	  
      WRITE(6,610) (BB(I),I=1,NFSHEL)
  610 FORMAT(10F8.5)
	  DO 22 K=1,NFSHEL
	     C(K)=0.0D0
		 DO 23 J=1,NFSHEL
		    C(K)=C(K)+BB(J)*FF(K,J)
   23    CONTINUE
         IF(DABS(C(K)-BS(K)).GT.1.0D-13) THEN
		    WRITE(6,611) K,C(K),BS(K),C(K)-BS(K)
			STOP ' STOP IN FNNMAT COEFF. IS INCORRECT'
  611       FORMAT(' STOP IN FNNMAT K=',I5,3E12.4)
         END IF
   22 CONTINUE
      DO 25 I=1,NFSHEL
	     B(I,IB,ISPIN)=BB(I)
   25 CONTINUE
   24 CONTINUE
   26 CONTINUE
	  RETURN
      END

	  
