      program clebsch
*---------------------------------------------------------------
*
*     Berechnet Clebsch-Gordan-Koeffizienten
*     
*---------------------------------------------------------------
      implicit none
      integer*4  N,DIM
      parameter  (N=2,DIM=1000)
      integer    sMal2MAT(N)
      integer*4  EVNumber,LastEVNumber,DimOfH
      integer*1  BasisElement(N,DIM)
      integer*1  TestElement(N),Element1(N),Element2(N)
      integer    Absenker,AbGes,MMaxMal2
      integer    AbsenkTabelle(100) ! Dim eigentlich Sum sMal2MAT(i)
      integer*4  AbGesListe(100)
      integer*4  k,L,m,i,j
      integer*4  Feldplatz
      integer*4  HDim
*     
      real*8     SSquareMatrix(Dim,Dim),SSquare,SumSSPlus1
      real*8     SSquareEigenWert(Dim),SSquareEigenvector(Dim,Dim)
      real*8     SGes,MGes
*
      real*8     wk(3*Dim+1)
      integer    Error
*
      real*8     Pi,C,HQuerC
      complex*16 ui
      character  Version*50
*
      common/const/Pi,C,HQuerC
*
      Pi     = 3.141592635d0
      HQuerC = 197.327053d0
      C      = 2.99792458d23
      ui     = DCMPLX(0.d0,1.d0)
      Version= 'clebsch-gordan.f 01 Juergen Schnack  25. 11. 2003'
*---------------------------------------------------------------
*
*     output von N und s
*
      write(6,'(a)') Version
      write(6,'(a,i4)') 'N = ',N
*---------------------------------------------------------------
*
*     Spinquantenzahlen
*
      sMal2MAT(1) = 2
      sMal2MAT(2) = 2
*     
      SumSSPlus1 = 0.d0
      do i=1,N
         SumSSPlus1 = SumSSPlus1 
     &        + 0.5d0*sMal2MAT(i)*(0.5d0*sMal2MAT(i)+1.d0)
      enddo
*
      DimOfH = 1
      MMaxMal2   = 0
      write(6,*) 'Spinquantenzahlen 2*s(i)'
      do i=1,N
         write(6,*) i,sMal2MAT(i)
         DimOfH = DimOfH * (sMal2MAT(i)+1)
         MMaxMal2   = MMaxMal2 + sMal2MAT(i)
      enddo
      write(6,*) ' '
      write(6,*) 'Dimension des Hilbertraumes = ',DimOfH
      write(6,*) 'max. M bzw. S               = ',0.5d0*MMaxMal2
*---------------------------------------------------------------
*
*
      write(6,'(a)') '***********************************************'
      write(6,'(a,a)') '   s_1-m_1  s_2-m_2  SMax-S  MMax-M  CGC'
*---------------------------------------------------------------
*
*     äußerste Schleife läuft über Absenkungen
*
*
*     Magnonenvakuum, d.h. Absenkung = 0 zuerst analytisch
*
      EVNumber = 1
      LastEVNumber = EVNumber

*
      write(6,'(4i8,f15.8)')
     &               0   ! s_1-m_1
     &              ,0   ! s_2-m_2
     &              ,0   ! SMax-S
     &              ,0   ! MMax-M
     &              ,1.d0  ! Clebsch-Gordon-Koefizient
*
      do AbGes=1,MMaxMal2
         MGes = 0.5d0*MMaxMal2 - AbGes
*
         do j=1,Abges
            AbsenkTabelle(j) = 1 ! alle Absenker auf Platz 1
         enddo
         do j=Abges+1,MMaxMal2
            AbsenkTabelle(j) = 0
         enddo
         FeldPlatz = 0
 1000    continue
         do j=1,N
            TestElement(j) = 0
         enddo
         do j=1,Abges
            TestElement(AbsenkTabelle(j)) =
     &           TestElement(AbsenkTabelle(j)) + 1
         enddo
         do j=1,N
            if (TestElement(j).gt.sMal2MAT(j)) goto 1100 ! kein gültiger Zustand
         enddo
         if(FeldPlatz.ge.DIM) goto 9999
         FeldPlatz = FeldPlatz + 1
         do j=1,N
            BasisElement(j,FeldPlatz) = TestElement(j)
         enddo
 1100    continue
*
         Absenker = 1
         AbsenkTabelle(Absenker) = AbsenkTabelle(Absenker) + 1
         if(AbsenkTabelle(Absenker).le.N) goto 1000  ! neuer Absenker o.k.
 1200    continue
         Absenker = Absenker + 1
         if(Absenker.gt.Abges) goto 1300
         AbsenkTabelle(Absenker) = AbsenkTabelle(Absenker) + 1
         if(AbsenkTabelle(Absenker).le.N) then
            do j=1,Absenker-1
               AbsenkTabelle(j) = AbsenkTabelle(Absenker)
            enddo
            goto 1000
         else
            goto 1200  
         endif
 1300    continue ! keine weiteren Absenkungen
*
         AbGesListe(AbGes) = FeldPlatz
*
*     SSquarematrix für Unterraum konstanter Absenkung berechnen
*
         do L=1,AbGesListe(AbGes)
            do j=1,N
               Element2(j) = BasisElement(j,L)
            enddo
            do m=1,AbGesListe(AbGes)
               do j=1,N
                  Element1(j) = BasisElement(j,m)
               enddo
               SSquareMatrix(m,L) = 
     &              2.d0*SSquare(N,sMal2MAT,Element1,Element2)
            enddo
         enddo
         do L=1,AbGesListe(AbGes) ! Diagonale dazuaddieren
            SSquareMatrix(L,L) = SSquareMatrix(L,L)
     &           + SumSSPlus1
         enddo
*     
*     Diagonalisieren
*
*     DSYEV
*
*
         call DSYEV('V','U',AbGesListe(AbGes),SSquareMatrix,Dim
     &        ,SSquareEigenwert,wk,3*Dim+1,Error)
*
*
*     Eigenvektor zum Eigenwert i: H(j,i) sind die Komponenten des EV i
*
         write(6,'(a)')
     &        '-----------------------------------------------'
         do i=1,AbGesListe(AbGes)
            do j=1,AbGesListe(AbGes)
               write(6,'(4i8,f15.8)') (BasisElement(k,j),k=1,2)
     &              ,(MMaxMal2
     &               -2*NINT(DSQRT(SSquareEigenwert(i)+0.25d0)-0.5d0))/2
     &              ,AbGes,SSquareMatrix(j,i)
c$$$               write(6,'(2i8,f7.2,i8,f15.8)') (BasisElement(k,j),k=1,2)
c$$$     &              ,DSQRT(SSquareEigenwert(i)+0.25d0)-0.5d0
c$$$     &              ,AbGes,SSquareMatrix(j,i)
            enddo
         enddo
*
         LastEVNumber = EVNumber
*
      enddo                     ! hier geht's zur nächsten Absenkung
*---------------------------------------------------------------
*
*     Schluss
*
 9999 continue
*---------------------------------------------------------------
*
      stop
      end
*===============================================================
*
*     Unterprogramme
*
*===============================================================
      real*8 function Overlap(N,Element1,Element2)
      implicit none
      integer*4  N,j
      integer*1    Element1(N),Element2(N)
*
      real*8     Test
*
      Test=0.0d0
      do j=1,N
         Test=Test+DABS(1.d0*(Element1(j)-Element2(j)))
      enddo
      if(Test.lt.1.d-6) then 
         Overlap=1.d0
      else 
         Overlap=0.d0
      endif
*
      return
      end
*===============================================================
      real*8 function SSquare(N,sMal2MAT,Element1,Element2)
*---------------------------------------------------------------
*
*---------------------------------------------------------------
      implicit none
      integer*4    N,i,j,k
      integer*1    TestElement(N),Element1(N),Element2(N)
      integer      sMal2MAT(N)
*
      real*8     s(N)
      real*8     SSq,Overlap
*
      SSq = 0.d0
*
      do i=1,N
         s(i) = 0.5d0 * sMal2MAT(i)
      enddo
*
*     Schleife über Spinplätze
*
      i=1
      j=2
*---------------------------------------------------------------
*     z - Komponente
      if (Overlap(N,Element1,Element2).gt.0.999d0) then
         SSq = SSq + (s(i)-Element1(i))*(s(j)-Element2(j))
      endif
*     x-, y-Komponenten
      do k=1,N
         TestElement(k)=Element2(k)
      enddo
      TestElement(j)=TestElement(j)+1
      TestElement(i)=TestElement(i)-1
      SSq = SSq + 0.5d0*Overlap(N,Element1,TestElement)
     &     * DSQRT((TestElement(j))*(2.d0*s(j)-TestElement(j)+1))
     &     * DSQRT((2.d0*s(i)-TestElement(i))*(TestElement(i)+1))
*
      do k=1,N
         TestElement(k)=Element2(k)
      enddo
      TestElement(j)=TestElement(j)-1
      TestElement(i)=TestElement(i)+1
      SSq = SSq + 0.5d0*Overlap(N,Element1,TestElement)
     &     * DSQRT((2.d0*s(j)-TestElement(j))*(TestElement(j)+1))
     &     * DSQRT((TestElement(i))*(2.d0*s(i)-TestElement(i)+1))
*---------------------------------------------------------------
*
      SSquare = SSq
*
      return
      end
