C************************************************************************
C      BAND PLOT PROGRAM                                                *
C               MODIFIED BY A.YANASE
C                         1995/02/10
C     01:D1 FILE                                                        *
C     02:ENERGY OUTPUT FILE                                             *
C     03:FORMAT FILE                                                    *
C     40:40-BAN FILE                                                    *
C************************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
      CHARACTER*28 D
C
      CALL TSPPRP
      WRITE(6,*) ' END OF TSPPRP'
C
      CALL STRUCT(NLCOMP,NSPIN,JMARK)
      WRITE(6,*) ' END OF STRUCT'
       write(6,*) 'jmark=',jmark
C
      CALL ENERED(NLCOMP)
      WRITE(6,*) ' END OF ENERED'
C
      IPRR=IPR
      JD=0
      IF(ISO.EQ.3) JD=1
      CALL AXENER(JD,IPRR,NSPIN)
      WRITE(6,*) ' END OF AXENR' 
C
c      CALL TEST(7)
      CALL DRFSTR(0)
      CALL DEGREE
      CALL ORIGIN(50.0,40.0)
      write(6,*) ' END OF ORIGIN'
C
      CALL PLOTR
      write(6,*) ' END OF PLOTR'
C
      IF(NSPIN.NE.2) IUD=1
      IF(NSPIN.EQ.2) IUD=2
      CALL PLOTSC(IUD)
      IF(JMARK.EQ.1) CALL MARKPL(IUD)
      IF(NSPIN.EQ.3) THEN
         CALL SETLN2(2.0,1.0,0.0,1)
         CALL PLOTSC(2)
         IF(JMARK.EQ.1) CALL MARKPL(2)
      END IF
      write(6,*) ' END OF PLOTSC'
      CALL SETPEN(6)
      READ(3,*,END=10)EF
      GOTO 20
   10 WRITE(6,*)' FERMI LEVEL'
      READ(5,*,ERR=30)EF
   20 CALL SETPEN(5)
      CALL FERMIE(EF)
      write(6,*) ' END OF FERMIE'
   30 READ(3,'(A28)',err=40)D
      GOTO 50
   40 WRITE(6,*)' TITLE <28 1****+****2****+***'
      READ(5,'(A28)',END=60)D
   50 CALL SETPEN(4)
      CALL TITLE(D,28)
      WRITE(6,*) ' END OF TITLE'
   60 CALL DRFEND
      STOP
      END
      SUBROUTINE TITLE(D,K)
C***********************************************************************
C WRITE TITLE                                                          *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
      CHARACTER*28 D
      REAL*4 YM
      YM=YMM
      CALL BTYPE(2.0,YM+5.0,7.0,D,0.0,K,1,1,-1.0)
      RETURN
      END
      SUBROUTINE FERMIE(EF)
C***********************************************************************
C WRITE FERMI ENERGY                                                   *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
      REAL*4 XO,XE,YF,XMM,YMM
      YMM=YM
      XMM=XM
      YF=(EF-EO)*WY
C         write(6,*) naxm,'fermie','yf=',yf,ymm
      IF(YF.GT.YMM.OR.YF.LT.0.0) RETURN
      XO=0.0
      DO 10 N=1,NAXM
C         write(6,*) n,icon(n),naxm
      IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
         XE=WX*PS(2,N)
C          write(6,*) xo,yf,xe
         CALL LINE3(XO,YF,XE,YF,1.0,2.0,1)
         IF(N.NE.NAXM) XO=WX*PS(1,N+1)
      ENDIF
   10 CONTINUE
      CALL BTYPE(XMM+2.0,YF-2.5,5.0,'E$DF',0.0,4,1,1,-1.0)
      RETURN
      END
      SUBROUTINE TSPPRP
C***********************************************************************
C PREPERATION TO CALL TSPACE                                           *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
      DIMENSION JB(2,3)
      READ(1,*)
      READ(1,*) IL,NGEN,INV
      CALL TSPACE(IL)
      DO 1 I=1,NGEN
      READ(1,*) JA,((JB(J,K),J=1,2),K=1,3)
      CALL TSGENR(JA,JB)
    1 CONTINUE
      CALL TSPGRP(INV)
      IF(IPR.GE.3) CALL TSPGDS
      READ(1,104) A,B,C
      READ(1,104) CA,CB,CC
  104 FORMAT(3D23.16)
      CALL TSLATC(A,B,C,CA,CB,CC)
C 102 FORMAT(14I5)
C     READ(1,102) NATOM,NKAT
C     NK=NKAT
C     CALL ENRRED(2,NK)
C     PRINT ,'  END OF ENRRED'
C     IF(IPR.GE.4)CALL ENERDS(6)
      RETURN
      END
      SUBROUTINE STRUCT(NLCOMP,NSPIN,JMARK)
C***********************************************************************
C SELECTION OF POINTS TO BE PLOTTED                                    *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)                
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX)
     &   ,MK(3,NAXMAX),NAXM
      COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
     &      ,ICC1M(NAXMAX),ICC2M(NAXMAX)
      COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
      DIMENSION KK1(3),KK2(3),KK3(3),MR(3),KTV(3),KG(3)
      CHARACTER*2 CCR1,CCR2,CCR3
      CHARACTER*4 CHARA4
      CHARACTER*10 MAGNET
C
      READ(3,*) MAGNET
      CHARA4=MAGNET
      IF(CHARA4.EQ.'NONM') ISO=1
      IF(CHARA4.EQ.'MAGN') ISO=2
      IF(CHARA4.EQ.'SPIN') ISO=3
      READ(3,*) NLCOMP,NSPIN
      IF((ISO.EQ.1.OR.ISO.EQ.3).AND.NSPIN.NE.0) THEN
         WRITE(6,*) ' FOR PARAMAGNETIC BANDS NSPIN MUST BE 0',NSPIN
         STOP
      END IF
      IF(ISO.EQ.2.AND.(NSPIN.LE.0.OR.NSPIN.GT.3)) THEN
         WRITE(6,*) ' FOR FERROMAGNETIC BANDS'
         WRITE(6,*) ' NSPIN=1 ONLY FOR UP   SPIN BANDS'
         WRITE(6,*) ' NSPIN=2 ONLY FOR DOWN SPIN BANDS'
         WRITE(6,*) ' NSPIN=3 FOR BOTH SPIN BANDS'
         STOP
      END IF
      READ(3,*) IPR,JMARK
      READ(3,*) EO,EM,YM,XM
      WY=YM/(EM-EO)
C
      PS(1,1)=0.0
      NK=0
      READ(3,*) NAXM
      IF(NAXM.GT.NAXMAX) THEN
         WRITE(6,*) ' NAXMAX=',NAXMAX,' NAXM=',NAXM
         STOP
      END IF
      DO 99 NAX=1,NAXM
      READ(3,*) (KK1(J),J=1,3),ICC1,(KK2(J),J=1,3),ICC2
      IF(IPR.GE.2) WRITE(6,600) (KK1(J),J=1,3),ICC1,(KK2(J),J=1,3),ICC2
  600 FORMAT( '(',3I4,')/',I4,'---','(',3I4,')/',I4)
C
      DO 75 J=1,3
      IIX=J
      IF(KK1(J)*ICC2.NE.KK2(J)*ICC1) GOTO 81
   75 CONTINUE
      write(6,*)  'KK1 AND KK2 ARE THE SAME POINT'
      STOP
   81 IX(NAX)=IIX
      DO 71 J=1,3
      KK1M(J,NAX)=KK1(J)
      KK2M(J,NAX)=KK2(J)
   71 CONTINUE
      ICC1M(NAX)=ICC1
      ICC2M(NAX)=ICC2
      KTV(1)=KK2(1)*ICC1-KK1(1)*ICC2
      KTV(2)=KK2(2)*ICC1-KK1(2)*ICC2
      KTV(3)=KK2(3)*ICC1-KK1(3)*ICC2
      ICC=ICC1*ICC2
      CALL ZZZY37(KTV(1),KTV(2),KTV(3),ICC,WW)
      RT(NAX)=WW
      BLK=0.0
      ICON(NAX)=0
      IF(NAX.NE.1) THEN
        CALL EQUIKK(KK1,ICC1,KK3,ICC3,KG,IND)
        IF(IND.EQ.0) THEN
            ICON(NAX-1)=1
            BLK=0.3*ABS(RT(1))
        ENDIF
        PS(1,NAX)=PS(2,NAX-1)+BLK
      ENDIF
      PS(2,NAX)=RT(NAX)+PS(1,NAX)
C      IF(KTV(IIX).LT.0.) RT(NAX)=-RT(NAX)
      IF(IPR.GE.2) WRITE(6,601) PS(1,NAX),PS(2,NAX),RT(NAX)
  601 FORMAT(3F10.5)
      DO 34 J=1,3
      KK3(J)=KK2(J)
   34 CONTINUE
      ICC3=ICC2
      KX=KK1(1)
      KY=KK1(2)
      KZ=KK1(3)
      ICC=ICC1
      CALL KPNAME(KX,KY,KZ,ICC,CCR1,KG)
      READ(CCR1,'(A2)')MR(1)
      MK(1,NAX)=MR(1)
      KX=KK2(1)
      KY=KK2(2)
      KZ=KK2(3)
      ICC=ICC2
      CALL KPNAME(KX,KY,KZ,ICC,CCR3,KG)
      READ(CCR3,'(A2)')MR(1)
      MK(3,NAX)=MR(1)
      KX=KK1(1)*ICC2+KK2(1)*ICC1
      KY=KK1(2)*ICC2+KK2(2)*ICC1
      KZ=KK1(3)*ICC2+KK2(3)*ICC1
      ICC=ICC1*ICC2*2
      CALL KPNAME(KX,KY,KZ,ICC,CCR2,KG)
      READ(CCR2,'(A2)')MR(1)
      MK(2,NAX)=MR(1)
      IF(IPR.GE.2) WRITE(6,602) CCR1,CCR2,CCR3
  602 FORMAT(3(3X,A2))
   99 CONTINUE
      WX=XM/PS(2,NAXM)
      RETURN
      END
      SUBROUTINE PLOTR
C***********************************************************************
C PLOT FRAMES AND SCALES                                               *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
      REAL*4 X,Y,XO,XE,YM,E
      WRITE(6,*) ' NAXM=',NAXM
      YM=YMM
      XO=0.
      DO 10 N=1,NAXM
      IF(IPR.GE.3) WRITE(6,600) N,PS(1,N),PS(2,N),ICON(N)
  600 FORMAT(I4,2F10.6,I2)
      IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
         XE=WX*PS(2,N)
         CALL LINE2(XO,0.0,2)
         CALL LINE2(XE,0.0,1)
         CALL LINE2(XE,YM,2)
         CALL LINE2(XO,YM,1)
         XO=WX*PS(1,N+1)
      ENDIF
   10 CONTINUE
      DO 301 I=1,NAXM
      JM=1
      IF(ICON(I).EQ.1.OR.I.EQ.NAXM) JM=2
      DO 302 J=1,JM
      X=PS(J,I)*WX
      CALL LINE2(X,0.0,2)
      CALL LINE2(X,YM,1)
  302 CONTINUE
  301 CONTINUE
c      WRITE(6,*) ' 301'
      CALL LINE2(-7.0,0.0,2)
      CALL LINE2(-10.0,0.0,1)
      CALL LINE2(-10.0,YM,1)
      CALL LINE2(-7.0,YM,1)
      ESDV=0.1
      IF(EM-EO.LT.0.1) ESDV=0.01
      E0A=(EO-0.00001)/ESDV
      IF(E0A.GE.0.0) I0=INT(E0A)+1
      IF(E0A.LT.0.0) I0=INT(E0A)
      EMA=(EM+0.00001)/ESDV
      IF(EMA.GE.0.0) IM=INT(EMA)
      IF(EMA.LT.0.0) IM=INT(EMA)-1
c      WRITE(6,*) IM,I0,ESDV,E0A,EMA
      DO 100 I=1,IM-I0+1
      E=(I0+I-1)*ESDV
      Y=(E-EO)*WY
      IF(ABS(Y).LT.1.0) GO TO 101
      IF(ABS(Y-YM).LT.1.0) GO TO 101
      CALL LINE2(-10.0,Y,2)
      CALL LINE2(-7.0,Y,1)
c      WRITE(6,*) I,E,Y
  101 IF(ESDV.GE.0.1)
     &CALL FTYPE2(-13.0,Y-2.5,5.0,E,0.0,4,1,2,1,-1.0)
      IF(ESDV.LT.0.1)
     &CALL FTYPE2(-13.0,Y-2.5,5.0,E,0.0,4,2,2,1,-1.0)
  100 CONTINUE
C      write(6,*) ' btype'
      CALL BTYPE(-35.0,YM/2.0,5.0,'Energy (Ry)',90.0,11,3,1,-1.0)
C      write(6,*) ' btype'
      DO 110 N=1,NAXM
      MMM=MK(1,N)
      X=WX*PS(1,N)
c      write(6,*) ' ctype',x
      CALL CTYPE(X,-10.0,5.0,MMM)
      MMM=MK(2,N)
      X=WX*(PS(1,N)+PS(2,N))/2.0
c      write(6,*) ' ctype',x
      CALL CTYPE(X,-7.5,4.0,MMM)
      IF(ICON(N).EQ.1.OR.N.EQ.NAXM) THEN
      MMM=MK(3,N)
      X=WX*PS(2,N)
c      write(6,*) ' ctype',x
      CALL CTYPE(X,-10.0,5.0,MMM)
      ENDIF
  110 CONTINUE
      RETURN
      END
      SUBROUTINE CTYPE(X,Y,H,M)
C***********************************************************************
C WRITE NAMES OF POINTS AND AXES                                       *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*4 D
      REAL*4 X,Y,H
      WRITE(D(1:2),'(A2)')M
      K=2
      IF(D(2:2).EQ.' ')K=1
      IF(D(1:2).EQ.'GM') THEN
      D(1:3)='$GC'
      K=3
      ELSEIF(D(1:2).EQ.'XD') THEN
      D(1:1)='X'
      K=1
      ELSEIF(D(1:2).EQ.'SM') THEN
      D(1:3)='$GS'
      K=3
      ELSEIF(D(1:2).EQ.'LD') THEN
      D(1:3)='$GL'
      K=3
      ELSEIF(D(1:2).EQ.'DT') THEN
      D(1:3)='$GD'
      K=3
      else if(D(1:2).eq.'TP') then
      D(1:2)="T'"
      else if(D(1:2).eq.'SP') then
      D(1:2)="S'"
      ELSEIF(D(1:2).EQ.'MX') THEN
      D(1:1)='M'
      K=1
      ELSEIF(D(1:2).EQ.'LX') THEN
      D(1:1)='L'
      K=1
      ENDIF
      CALL BTYPE(X,Y,H,D(1:K),0.0,K,3,1,-1.0)
      RETURN
      END
      SUBROUTINE PLOTSC(IUD)
C***********************************************************************
C PLOT BAND                                                            *
C***********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (MAXFKP=500,MAXIRP=33)
      PARAMETER (N65=65)
      COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
     &      ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),KPM
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/SCL/WX,WY,EO,EM,XM,YMM,IPR,ISO
      REAL*4 YM,XC,XA,YA,XB,YB
      CHARACTER*4 DDD
      YM=YMM
      DO 110 K=1,KPM
C        write(6,*) iud,k,jud(k)
      IF(IUD.NE.JUD(K)) GO TO 110
      NP=NKPT(K)
      WRITE(DDD(1:2),'(A2)') MK(2,JRR(K))
      IF(IPR.GE.1) WRITE(6,600) K,NP,DDD,NLIN(K)
  600 FORMAT(2I5,1X,A2,I4)
      IF(NLIN(K).LE.0) GO TO 110
      DO 100 J=1,NLIN(K)
      IF(IPR.GE.4) WRITE(6,*) ' BNAD NO=',J
      XA=WX*PS(1,JRR(K))
      YA=(EAX(1,J,K)-EO)*WY
      IF(YA.GT.0.0.AND.YA.LT.YM) CALL LINE4(XA,YA,2)
C      IF(IPR.GE.2) WRITE(6,*) '2',XA,YA
      DO 180 I=2,NP
C
      XB=XA
      YB=YA
      XA=PS(1,JRR(K))+(RT(JRR(K))*(I-1))/(NP-1)
      XA=XA*WX
      YA=(EAX(I,J,K)-EO)*WY
C      write(6,*) xa,ya,eax(i,j,k)
      IF(YB.GT.YM) GO TO 291
      IF(YA.LT.YM) GO TO 295
      XC=XB+(XA-XB)*((YM-YB)/(YA-YB))
      CALL LINE4(XC,YM,1)
      GO TO 180
  291 IF(YA.GT.YM) GO TO 180
      XC=XB+(XA-XB)*((YB-YM)/(YB-YA))
      CALL LINE4(XC,YM,2)
      GO TO 299
  295 IF(YB.LT.0.0) GO TO 296
      IF(YA.GT.0.0) GO TO 299
      XC=XB+(XA-XB)*(YB/(YB-YA))
      CALL LINE4(XC,0.0,1)
      GO TO 180
  296 IF(YA.LT.0.0) GO TO 180
      XC=XB+(XA-XB)*(-YB/(YA-YB))
      CALL LINE4(XC,0.0,2)
  299 CALL LINE4(XA,YA,1)
C      WRITE(6,*) '1',XA,YA
  180 CONTINUE
  100 CONTINUE
  110 CONTINUE
      RETURN
      END
      SUBROUTINE ENERED(NLCOMP)
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*2 MRR,MRRM
      COMMON/SCL/WX,WY,EO,EM,XM,YM,IPR,ISO
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
      COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
      DIMENSION KX(3),kxx(3),EIG(MAXEIG),KBB(3,10),KG(3,10),itcr(12)
      REWIND 2 
C
      N=0
      ifrag=0
   40 READ(2,150,END=41) KK,(KX(J),J=1,3),IC,IUD,MRR,MRN,MWEI,NST,NEIG
      IF(NEIG.GT.MAXEIG) THEN
         WRITE(6,*) ' MAXEIG=',MAXEIG,' NEIG=',NEIG
         STOP
      END IF
      IF(NEIG.GT.0) READ(2,*) (EIG(J),J=1,NEIG)
  150 FORMAT(I3,2X,3I3,1X,2I3,2X,A2,I2,3I3)
C  151 FORMAT(8F8.4)
      IF(NEIG.EQ.0) GO TO 40
      IF(NLCOMP.NE.0) THEN
         NFACT=(NLCOMP-1)/2+1
         DO 50 K=1,NEIG*NFACT
   50       READ(2,*)
      END IF
      N=N+1
      if(ifrag.eq.1) go to 40
      IF(N.GT.MAXKPT) THEN
         ifrag=1
         go to 40
      end if
      CALL NEAREC(KX,IC,KBB,KG,NG)
      DO 10 J=1,3
         kxx(j)=kbb(j,1)
         KXM(J,N)=KBB(J,1)
   10 CONTINUE
      icx=ic
      jdob=0
      if(iso.eq.3) jdob=1
      call corres(kxx,icx,kx,ic,jdob,itcr,nrr,ind)
C      write(6,*) kxx,kx,ic,(itcr(j),j=1,nrr)
      if(ind.eq.0) then
         write(6,*) ' kx=',kx,ic,' kxx=',kxx,icx
         write(6,*) ' stop at aorres in ENERED'
         stop
      end if
      ICM(N)=IC
C       write(6,*) n,(kxm(j,n),j=1,3),icm(n),KX,IC
      IUDM(N)=IUD
      MRRM(N)=MRR
      MRNM(N)=itcr(MRN)
      MWEIM(N)=MWEI
      NSTM(N)=NST
      NEIGM(N)=NEIG
      DO 20 J=1,NEIG
        EIGM(J,N)=EIG(J)
   20 CONTINUE
      GO TO 40
   41 CONTINUE
      if(ifrag.eq.1) then
         WRITE(6,*) ' MAXKPT=',MAXKPT,' NKPINT=',N
         STOP
      END IF
      NKPINT=N
      DO 1 I=1,NKPINT
         NSIT(I)=0
    1 CONTINUE
      RETURN
      END
      SUBROUTINE MARKPL(IUD)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (NAXMAX=20)                
      CHARACTER*2 MRRM
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
      COMMON/SCL/WX,WY,EO,EM,XMM,YMM,IPR,ISO
      COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)

      REAL*4 XM,YM,XA,YA
      DIMENSION NMARK(12)
      DATA NMARK/1,2,3,4,5,6,7,8,9,11,12,13/
C      DATA NMARK/73, 74, 75, 76, 80, 63, 68, 94, 70, 69, 72, 71/
C
      XM=XMM
      YM=YMM
      DO 10 I=1,NKPINT
         IF(IUD.NE.IUDM(I)) GO TO 10
C       write(6,*) i,mrrm(i),mrnm(i),nsit(i)
         IF(NSIT(I).EQ.0) GO TO 10
         DO 11 N=1,NSIT(I)
            XA=(PS(1,MSIT(N,I))+RT(MSIT(N,I))*XIMS(N,I))*WX
            DO 12 KK=1,NEIGM(I)
               YA=(EIGM(KK,I)-EO)*WY
               IF(YA.LT.YM.AND.YA.GT.0.0) THEN
                  CALL MARK(XA,YA,4.0,NMARK(MRNM(I)))
               END IF
   12       CONTINUE
   11    CONTINUE
   10 CONTINUE
      RETURN
      END  
      SUBROUTINE AXENER(JDOUB,IPR,NSPIN)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)                
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (MAXFKP=500,MAXIRP=33)
      PARAMETER (N65=65)
      CHARACTER*2 MRRM
      COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
     &      ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
     &      ,ICC1M(NAXMAX),ICC2M(NAXMAX)
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
      DIMENSION KB(3),JTR(12),IPA(12),ND(12)
      DIMENSION NRECPT(3,20)
      JD=JDOUB
      NAXEN=0
      DO 11 I=1,NAXM
         II=I
         IF(IPR.GE.2) WRITE(6,601) II
  601    FORMAT(' RECFND',I3)
         CALL RECFND(II,NRECPT,NRPT,IPR)
         IF(IPR.GE.2) WRITE(6,602) II,NRPT
  602    FORMAT(' KPFIEN',I3,' NUMBER OF RECPR. LAT. POINTS=',I3)
         CALL KPFIEN(II,NRECPT,NRPT,IPR)
         DO 21 K=1,3
            KB(K)=KK1M(K,I)*ICC2M(I)+KK2M(K,I)*ICC1M(I)
   21    CONTINUE
         IC=ICC1M(I)*ICC2M(I)*2
         CALL TSIREP(KB,IC,JD)
         CALL DGTRMD(JDO,NR,NH,NSTR,ND,MTR,JTR,IPA)
         IF(IPR.GE.2) WRITE(6,603) KB,IC
  603    FORMAT(' CMPTBL (',3I3,')/',I4)
         CALL CMPTBL(JD,II,KB,IC,NR,NH,MTR,IPR)   
         DO 22 IR=1,NR
            IF(JTR(IR).EQ.0.AND.IPA(IR).LT.IR) GO TO 22
            IIR=IR
            NAXEN=NAXEN+1
            IF(NAXEN.GT.MAXIRP) THEN
               WRITE(6,*) ' MAXIRP=',MAXIRP,' NAXEN=',NAXEN
               STOP
            END IF
            IF(NSPIN.NE.2) IUD=1
            IF(NSPIN.EQ.2) IUD=2 
            CALL INTPOR(JD,II,KB,IC,NH,IIR,IUD,IPR)
            JUD(NAXEN)=IUD
            IF(NSPIN.EQ.3) THEN
               NAXEN=NAXEN+1
               IF(NAXEN.GT.MAXIRP) THEN
                  WRITE(6,*) ' MAXIRP=',MAXIRP,' NAXEN=',NAXEN
                  STOP
               END IF
               IUD=2
               CALL INTPOR(JD,II,KB,IC,NH,IIR,IUD,IPR)
               JUD(NAXEN)=IUD
            END IF
   22    CONTINUE
   11 CONTINUE
      RETURN
      END
      SUBROUTINE RECFND(II,NRECPT,NRPT,IPR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)                
C
      COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
     &      ,ICC1M(NAXMAX),ICC2M(NAXMAX)
C
      DIMENSION KB(3),NRECPT(3,20),KG(3,10),KBB(3,10)
      JJ=0
      DO 10 J=1,4
         IF(J.EQ.1) THEN
            DO 11 K=1,3
               KB(K)=2*KK1M(K,II)*ICC2M(II)-KK2M(K,II)*ICC1M(II)
   11       CONTINUE
            ICC=ICC2M(II)*ICC1M(II)
         ELSE IF(J.EQ.2) THEN
            DO 12 K=1,3
               KB(K)=KK1M(K,II)
   12       CONTINUE
            ICC=ICC1M(II)
         ELSE IF(J.EQ.3) THEN
            DO 13 K=1,3
               KB(K)=KK2M(K,II)
   13       CONTINUE
            ICC=ICC2M(II)
         ELSE IF(J.EQ.4) THEN
            DO 14 K=1,3
               KB(K)=2*KK2M(K,II)*ICC1M(II)-KK1M(K,II)*ICC2M(II)
   14       CONTINUE
            ICC=ICC1M(II)*ICC2M(II)
         END IF
         CALL NEAREC(KB,ICC,KBB,KG,NNG)
         DO 16 N=1,NNG
            IF(JJ.EQ.0) GO TO 17
            DO 18 IJ=1,JJ
               DO 19 K=1,3
                  IF(KG(K,N).NE.NRECPT(K,IJ)) GO TO 18
   19          CONTINUE
               GO TO 16
   18       CONTINUE
   17       JJ=JJ+1
            DO 15 K=1,3
               NRECPT(K,JJ)=KG(K,N)
   15       CONTINUE
            IF(IPR.GE.4) WRITE(6,601) II,J,JJ,(NRECPT(K,JJ),K=1,3)
  601       FORMAT(6I5)
   16    CONTINUE
   10 CONTINUE
      NRPT=JJ
      RETURN
      END
      SUBROUTINE KPFIEN(II,NRECPT,NRPT,IPR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)                
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (MAXFKP=500,MAXIRP=33)
      CHARACTER*2 MRRM
C
      COMPLEX*16 CW
      COMMON/SPG2/IL,NG,IG(48),JV(2,3,48)
      COMMON/STK/KS(3,48),JS(48),NS,ICBB,CW(48,12)
C
      COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
     &      ,ICC1M(NAXMAX),ICC2M(NAXMAX)
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
      COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
     &   ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
C
      DIMENSION KB(3),KBC(3)
      DIMENSION NRECPT(3,20),JIND(MAXKPT),IMM(MAXKPT)
C
      CALL ADDINV(NTLATC)
      X1=DBLE(KK1M(1,II))/ICC1M(II)
      Y1=DBLE(KK1M(2,II))/ICC1M(II)
      Z1=DBLE(KK1M(3,II))/ICC1M(II)
      DX=DBLE(KK2M(1,II))/ICC2M(II)-X1
      DY=DBLE(KK2M(2,II))/ICC2M(II)-Y1
      DZ=DBLE(KK2M(3,II))/ICC2M(II)-Z1
C       write(6,*) x1,y1,z1
C       write(6,*) dx,dy,dz
      D6=1.0D-6
      NPOINT=0
      DO 10 I=1,NKPINT
         JIND(I)=0
   10 CONTINUE
      DO 1 I=1,NKPINT
         IF(JIND(I).NE.0) GO TO 1
         KB(1)=KXM(1,I)
         KB(2)=KXM(2,I)
         KB(3)=KXM(3,I)
         IC=ICM(I)
c          write(6,*) 'kb',kb,ic
         NJ=1
         IMM(1)=I
         IF(I.LE.NKPINT) THEN
            DO 4 JI=I+1,NKPINT
               IF(JIND(JI).NE.0) GO TO 4
               DO 5 K=1,3
                  IF(KB(K)*ICM(JI).NE.KXM(K,JI)*IC) GO TO 4
    5          CONTINUE
               NJ=NJ+1
               IMM(NJ)=JI
    4       CONTINUE
         END IF
C-----------------------------------------
C        STAR OF K IS OBTAINED
         CALL TSIREP(KB,IC,0)
         CALL ZZZY38
C------------------------------------------
         DO 2 J=1,NS
c          write(6,*) j,(ks(k,j),k=1,3),icbb
           DO 3 KGG=1,NRPT
              XK=DBLE(KS(1,J))/ICBB+DBLE(NRECPT(1,KGG))
              YK=DBLE(KS(2,J))/ICBB+DBLE(NRECPT(2,KGG))
              ZK=DBLE(KS(3,J))/ICBB+DBLE(NRECPT(3,KGG))
c           write(6,*) xk,yk,zk
              IF(DABS(DX).LE.D6.AND.DABS(XK-X1).GT.D6) GO TO 3 
              IF(DABS(DY).LE.D6.AND.DABS(YK-Y1).GT.D6) GO TO 3
              IF(DABS(DZ).LE.D6.AND.DABS(ZK-Z1).GT.D6) GO TO 3
              IF(DABS(DX).GT.D6) THEN
                 IF(DABS(DY).GT.D6) THEN
                    IF(DABS((XK-X1)/DX-(YK-Y1)/DY).GT.D6) GO TO 3
                 END IF
                 IF(DABS(DZ).GT.D6) THEN
                    IF(DABS((XK-X1)/DX-(ZK-Z1)/DZ).GT.D6) GO TO 3
                 END IF
              END IF
              IF(DABS(DY).GT.D6.AND.DABS(DZ).GT.D6) THEN
                 IF(DABS((YK-Y1)/DY-(ZK-Z1)/DZ).GT.D6) GO TO 3
              END IF 
              IF(IX(II).EQ.1.AND.
     &         ((XK-X1)/DX.LT.-1.0D0.OR.(XK-X1)/DX.GT.2.0D0)) GO TO 3
              IF(IX(II).EQ.2.AND.
     &         ((YK-Y1)/DY.LT.-1.0D0.OR.(YK-Y1)/DY.GT.2.0D0)) GO TO 3
              IF(IX(II).EQ.3.AND.
     &         ((ZK-Z1)/DZ.LT.-1.0D0.OR.(ZK-Z1)/DZ.GT.2.0D0)) GO TO 3
              KBC(1)=KS(1,J)+NRECPT(1,KGG)*ICBB
              KBC(2)=KS(2,J)+NRECPT(2,KGG)*ICBB
              KBC(3)=KS(3,J)+NRECPT(3,KGG)*ICBB
              IF(NPOINT.EQ.0) GO TO 20
              DO 21 N=1,NPOINT
                 DO 22 K=1,3
                    IF(KBC(K)*ICCB(N).NE.KKB(K,N)*ICBB) GO TO 21
   22            CONTINUE
                 IF(MRNM(IFKP(N)).EQ.MRNM(I)) GO TO 3
   21         CONTINUE
   20         DO 6 JI=1,NJ
                 JIND(IMM(JI))=I
                 NPOINT=NPOINT+1
                 IF(NPOINT.GT.MAXFKP) THEN
                    WRITE(6,*) ' MAXFKP=',MAXFKP,' NPOINT=',NPOINT
                    STOP
                 END IF
                 IFKP(NPOINT)=IMM(JI)
                 JSKP(NPOINT)=JS(J)
                 KKB(1,NPOINT)=KBC(1)
                 KKB(2,NPOINT)=KBC(2)
                 KKB(3,NPOINT)=KBC(3)
                 ICCB(NPOINT)=ICBB
                 IF(IPR.GE.4) THEN
                    N=NPOINT
                    WRITE(6,601) N,KBC,ICBB,IFKP(N)
     &                          ,KB,IC,MRRM(IFKP(N)),MRNM(IFKP(N))
     &                     ,j,kgg,(KS(k,J),K=1,3)
  601               FORMAT(I3,'(',3I4,')/',I4,I3,
     &                     '(',3I3,')/',I3,2X,A2,I3,' J='I3,' KGG='
     &                         ,I3,'KS=',3i4)
                 END IF
    6         CONTINUE
c             write(6,*) xk,yk,zk
c             write(6,*) npoint,kbc,icbb,i,j,kgg
    3      CONTINUE
    2   CONTINUE
    1 CONTINUE
      CALL REMINV
      RETURN
      END
      SUBROUTINE CMPTBL(JD,II,KBB,ICC,NRT,NH,MTR,IPR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)    
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (MAXFKP=500,MAXIRP=33)
      PARAMETER (N65=65)
      CHARACTER*2 MRRM
      COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
     &      ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
      COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
     &   ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
C
      DIMENSION KB(3),KBX(3),KBB(3)
      DIMENSION ITCR(12),ND(12),JTR(12),IPA(12)
      DIMENSION NDX(12),JTRX(12),IPAX(12)
      DIMENSION ICP(12,12),ICP0(12,12)
C
      JDOB=JD
      DO 11 I=1,NPOINT
         INDM(I)=0
   11 CONTINUE
      DO 1 I=1,NPOINT
         IF(INDM(I).NE.0) GO TO 1
         KB(1)=KKB(1,I)
         KB(2)=KKB(2,I)
         KB(3)=KKB(3,I)
         IC=ICCB(I)
         KBX(1)=KXM(1,IFKP(I))
         KBX(2)=KXM(2,IFKP(I))
         KBX(3)=KXM(3,IFKP(I))
         ICX=ICM(IFKP(I))
         CALL CORRES(KB,IC,KBX,ICX,JDOB,ITCR,NRR,IND)
         IF(IND.EQ.0) THEN
            write(6,*) ' kb=',kb,ic,' kbx=',kbx,icx
            WRITE(6,*) ' STOP AT CORRES IN CMPTBL'
            STOP
         END IF
         CALL TSIREP(KBX,ICX,JDOB)
         CALL DGTRMD(JDO,NRX,NHX,NSTRX,NDX,MTRX,JTRX,IPAX)
         CALL TSIREP(KB,IC,JDOB)
         CALL DGTRMD(JDO,NRR,NHH,NSTR,ND,MTRR,JTR,IPA)
         DO 4 IR1=1,NRX
            IF(JTR(IR1).EQ.0.AND.IR1.LT.IPAX(IR1).AND.
     &           ITCR(IR1).GT.IPA(ITCR(IR1))) THEN
                 IW1=ITCR(IR1)
                 ITCR(IR1)=ITCR(IPAX(IR1))
                 ITCR(IPAX(IR1))=IW1
            END IF
    4    CONTINUE
C          write(6,*) ' cor',kb,ic,kbx,icx,(itcr(L),L=1,nrr)
         IF(NH.NE.NHH.OR.MTR.NE.MTRR) THEN
            CALL COMPAT(KBB,ICC,NRTT,KB,IC,NRR,JDOB,ICP,INDC)
C            write(6,*) kbb,icc,kb,ic
C            do 1000 i1=1,nrt
C            write(6,*) ' icp',(icp(L,i1),L=1,nrr)
C 1000       continue
            IF(INDC.NE.0) THEN
               WRITE(6,*) 'STOP AT COMPAT IN COMTBL'
               STOP
            END IF
            CALL CMPTRV(KBB,ICC,KB,IC,JDOB,ICP,ICP0)
C             do 1001 i1=1,nrt
C             write(6,*) ' icp0',(icp0(L,i1),L=1,nrr)
C 1001        continue
            DO 2 IIR=1,NRT
               ICPTBL(IIR,I)=ICP0(ITCR(MRNM(IFKP(I))),IIR)
    2       CONTINUE
            INDM(I)=1
            IF(IPR.GE.4) WRITE(6,601) I,INDM(I),(ICPTBL(K,I),K=1,NRT)
  601       FORMAT(14I3)
            IF(I.EQ.NPOINT) GO TO 1
            DO 20 J=I+1,NPOINT
               DO 22 K=1,3
                  IF(KB(K)*ICCB(J).NE.KKB(K,J)*IC) GO TO 1
   22          CONTINUE
               DO 23 IIR=1,NRT
                  ICPTBL(IIR,J)=ICP0(ITCR(MRNM(IFKP(J))),IIR)
   23          CONTINUE
               INDM(J)=1
               IF(IPR.GE.4)
     &            WRITE(6,601) J,INDM(J),(ICPTBL(K,J),K=1,NRT)
   20       CONTINUE
         ELSE
            DO 3 IIR=1,NRT
               ICPTBL(IIR,I)=0
    3       CONTINUE
            ICPTBL(ITCR(MRNM(IFKP(I))),I)=1
            INDM(I)=-1
            IF(IPR.GE.4)
     &            WRITE(6,601) I,INDM(I),(ICPTBL(K,I),K=1,NRT)
            IF(I.EQ.NPOINT) GO TO 1
            DO 30 J=I+1,NPOINT
               DO 32 K=1,3
                  IF(KB(K)*ICCB(J).NE.KKB(K,J)*IC) GO TO 1
   32          CONTINUE
               DO 33 IIR=1,NRT
                  ICPTBL(IIR,J)=0
   33          CONTINUE
               ICPTBL(ITCR(MRNM(IFKP(J))),J)=1
               INDM(J)=-1
               IF(IPR.GE.4)
     &            WRITE(6,601) J,INDM(J),(ICPTBL(K,J),K=1,NRT)
   30       CONTINUE
         END IF
    1 CONTINUE
      RETURN
      END
      SUBROUTINE INTPOR(JD,II,KBB,ICC,NH,IIR,IUD,IPR)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NAXMAX=20)                
      PARAMETER (MAXKPT=1200,MAXEIG=120)
      PARAMETER (MAXFKP=500,MAXIRP=33)
      PARAMETER (N65=65)
      CHARACTER*2 MRRM,MRR
      COMMON/AXE/EAX(N65,MAXEIG,MAXIRP),NKPT(MAXIRP)
     &      ,NLIN(MAXIRP),JRR(MAXIRP),JUD(MAXIRP),NAXEN
      COMMON/AXP/IFKP(MAXFKP),JSKP(MAXFKP),KKB(3,MAXFKP)
     &   ,ICCB(MAXFKP),ICPTBL(12,MAXFKP),INDM(MAXFKP),NPOINT
      COMMON/ISP/IX(NAXMAX),RT(NAXMAX),PS(2,NAXMAX),ICON(NAXMAX),
     &          MK(3,NAXMAX),NAXM
      COMMON/KKEND/KK1M(3,NAXMAX),KK2M(3,NAXMAX)
     &      ,ICC1M(NAXMAX),ICC2M(NAXMAX)
      COMMON/ENR/KXM(3,MAXKPT),ICM(MAXKPT),IUDM(MAXKPT)
     &   ,MRRM(MAXKPT),MRNM(MAXKPT),MWEIM(MAXKPT),NSTM(MAXKPT)
     &   ,NEIGM(MAXKPT),EIGM(MAXEIG,MAXKPT),NKPINT
      COMMON/MRK/NSIT(MAXKPT),MSIT(3,MAXKPT),XIMS(3,MAXKPT)
C
      DIMENSION KBB(3)
      DIMENSION XX(N65),YY(N65),X(MAXFKP),YM(MAXEIG,MAXFKP)
      DIMENSION SM(MAXFKP),JMAX(MAXFKP)
      DIMENSION IXM(MAXFKP),Y(MAXFKP),XXI(MAXFKP)
C
      D6=1.0D-6
      IF(IPR.GE.2) THEN
        KX=KBB(1)
        KY=KBB(2)
        KZ=KBB(3)
        CALL KPNAME(KX,KY,KZ,ICC,MRR,KG)
        WRITE(6,600) MRR,IIR
  600   FORMAT(' INTERPORATION START FOR ',A2,I2)
      END IF
      JRR(NAXEN)=II
      JDOB=JD
      WIDE=DBLE(KK2M(IX(II),II))/ICC2M(II)
     &    -DBLE(KK1M(IX(II),II))/ICC1M(II)
      NXI=0
      DO 1 I=1,NPOINT
         IF(IUD.NE.IUDM(IFKP(I))) GO TO 1
         XI=(DBLE(KKB(IX(II),I))/ICCB(I)
     &      -DBLE(KK1M(IX(II),II))/ICC1M(II))/WIDE
C       write(6,*) i,xi,nxi
         IF(IIR.NE.1) GO TO 13
         IF(XI+D6.LT.0.0) GO TO 13
         IF(XI-D6.GT.1.0D0) GO TO 13
         IF(XI+D6.LT.1.0D0.OR.II.EQ.NAXM.OR.ICON(II).EQ.1) THEN
            NSIT(IFKP(I))=NSIT(IFKP(I))+1
            IF(NSIT(IFKP(I)).GT.3) THEN
               WRITE(6,*) ' NSIT MUST BE LESS THAN 3',NSIT(IFKP(I))
               STOP
            END IF
            MSIT(NSIT(IFKP(I)),IFKP(I))=II
            XIMS(NSIT(IFKP(I)),IFKP(I))=XI
         END IF
C         write(6,*) i,ifkp(i),nsit(ifkp(i)),ii,naxm,xi,'13up'
   13    CONTINUE
         IF(ICPTBL(IIR,I).EQ.0) GO TO 1
         IF(NXI.EQ.0) GO TO 11
         DO 15 NI=1,NXI
            IF(DABS(XI-X(NI)).LT.1.0D-6) GO TO 12
   15    CONTINUE
   11    NXI=NXI+1
         JMAX(NXI)=0
         X(NXI)=XI
         NI=NXI
   12    continue
C   12    IF(INDM(I).EQ.1) THEN
            DO 5 L=1,NEIGM(IFKP(I))
               DO 2 K=1,ICPTBL(IIR,I)
                  IF(JMAX(NI).EQ.0) GO TO 51
                  DO 52 JJ=1,JMAX(NI)
                     J1=JMAX(NI)-JJ+1
                     IF(YM(J1,NI).LT.EIGM(L,IFKP(I))) GO TO 53
                     YM(J1+1,NI)=YM(J1,NI)
   52             CONTINUE
   51             J1=0
   53             YM(J1+1,NI)=EIGM(L,IFKP(I))
                  JMAX(NI)=JMAX(NI)+1
    2          CONTINUE
    5       CONTINUE
C         ELSE
C            DO 54 L=1,NEIGM(IFKP(I))
C               YM(L,NI)=EIGM(L,IFKP(I))
C   54       CONTINUE
C            JMAX(NI)=NEIGM(IFKP(I))
C         END IF
    1 CONTINUE
      DO 6 I=1,NXI
         IXM(I)=I
    6 CONTINUE
      XXI(1)=X(1)
      DO 61 I=2,NXI
         DO 62 J=1,I-1
            JJ=I-J
            IF(XXI(JJ).LT.X(I)) GO TO 63
            XXI(JJ+1)=XXI(JJ)
            IXM(JJ+1)=IXM(JJ)
   62    CONTINUE
         JJ=0
   63    JJ=JJ+1
         XXI(JJ)=X(I)
         IXM(JJ)=I
   61 CONTINUE
C        write(6,*) ' xxi',(xxi(i),i=1,nxi)
      NLIN(NAXEN)=0
      NKPT(NAXEN)=N65
      IF(NXI.LT.3) GO TO 73
      JJMAX=0
      DO 21 I=1,NXI
         IF(JMAX(I).GT.JJMAX) JJMAX=JMAX(I)
   21 CONTINUE
C        write(6,*) 'jjmax=',jjmax
C        write(6,*) (jmax(ixm(i)),i=1,nxi)
      DO 70 J=1,JJMAX 
         NNN=0
         NCO=0
         DO 71 I=1,NXI
            IF(JMAX(IXM(I)).LT.J) GO TO 71
            NNN=NNN+1
            X(NNN)=XXI(I)
            IF(XXI(I)+D6.GE.0.0.AND.XXI(I)-D6.LT.1.0D0) NCO=NCO+1
            Y(NNN)=YM(J,IXM(I))
   71    CONTINUE
         IF(NCO.LT.3) GO TO 70
         IF(IPR.GE.4) THEN
            WRITE(6,*) ' BAND NO=',J
            WRITE(6,666) (X(I),I=1,NNN)
            WRITE(6,666) (Y(I),I=1,NNN)
  666       FORMAT(10F7.4)
         END IF
         C1=4.0*((Y(2)-Y(1))/(X(2)-X(1)))
         AMU1=2.0
         CN=4.0*((Y(NNN)-Y(NNN-1))/(X(NNN)-X(NNN-1)))
         ALMN=2.0
         NLIN(NAXEN)=J
         INIT=0
         ICO=0
         DO 100 I=1,N65
            DD=DBLE(I-1)/(N65-1)
            IF(DD.LT.X(1)) GO TO 100
            IF(DD.GT.X(NNN)) GO TO 101
            IF(INIT.EQ.0) INIT=I
            ICO=ICO+1
            XX(ICO)=DD
  100    CONTINUE
  101    CONTINUE
         CALL S3N(X,Y,SM,XX,YY,NNN,ICO,C1,CN,AMU1,ALMN)
         DO 72 I=1,N65
            EAX(I,J,NAXEN)=99.0
   72    CONTINUE
         DO 75 I=1,ICO
            EAX(INIT+I-1,J,NAXEN)=YY(I)
   75    CONTINUE
   70 CONTINUE
   73 RETURN
      END
      SUBROUTINE S3N(X,Y,SM,XX,YY,N,NN,C1,CN,AMU1,ALMN)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(99),Y(99),SM(99),XX(200),YY(200)
      DIMENSION H(99),ALM(99),AMU(99),C(99),P(99),Q(99),U(99)
      N1=N-1
      DO 110 I=2,N
      H(I)=X(I)-X(I-1)
  110 CONTINUE
      DO 120 I=2,N1
      ALM(I)=H(I+1)/(H(I)+H(I+1))
      AMU(I)=1.0-ALM(I)
  120 CONTINUE
      DO 130 I=2,N1
      C(I)=3.0*(ALM(I)*(Y(I)-Y(I-1))/H(I)+AMU(I)*(Y(I+1)-Y(I))/H(I+1))
  130 CONTINUE
      C(1)=C1
      C(N)=CN
      AMU(1)=AMU1
      ALM(N)=ALMN
      P(1)=2.0
      Q(1)=-AMU(1)/P(1)
      U(1)=C(1)/P(1)
      DO 140 K=2,N
      P(K)=ALM(K)*Q(K-1)+2.0
      Q(K)=-AMU(K)/P(K)
      U(K)=(C(K)-ALM(K)*U(K-1))/P(K)
  140 CONTINUE
      SM(N)=U(N)
      DO 150 K=1,N1
      K1=N1-K+1
      SM(K1)=Q(K1)*SM(K1+1)+U(K1)
  150 CONTINUE
      DO 160 I=1,NN
      XXI=XX(I)
      DO 170 K=2,N
      IF(XXI.GT.X(K)) GO TO 170
      J1=K
      GO TO 180
  170 CONTINUE
  180 J=J1-1
      SMJ=SM(J)
      SMJ1=SM(J1)
      YJ=Y(J)
      YJ1=Y(J1)
      HJ1=H(J1)
      XJ1=X(J1)-XXI
      XJ=XXI-X(J)
      HJ2=HJ1*HJ1
      HJ3=HJ2*HJ1
      YY(I)=SMJ*XJ1*XJ1*XJ/HJ2-SMJ1*XJ*XJ*XJ1/HJ2+YJ*XJ1*XJ1*(2.0*XJ+
     &      HJ1)/HJ3+YJ1*XJ*XJ*(2.0*XJ1+HJ1)/HJ3
  160 CONTINUE
      RETURN
      END
