  L=0  JL1=1  JL2=1    OVER= 1.000D+00
 
  L=1  JL1=1  JL2=1    OVER= 1.000D+00
  L=1  JL1=2  JL2=1    OVER=-3.590D-12
  L=1  JL1=2  JL2=2    OVER= 1.000D+00
  L=1  JL1=3  JL2=1    OVER=-3.733D-18
  L=1  JL1=3  JL2=2    OVER= 2.112D-17
  L=1  JL1=3  JL2=3    OVER= 1.000D+00
 
  L=2  JL1=1  JL2=1    OVER= 1.000D+00
  L=2  JL1=2  JL2=1    OVER= 6.597D-10
  L=2  JL1=2  JL2=2    OVER= 1.000D+00
  L=2  JL1=3  JL2=1    OVER= 2.073D-12
  L=2  JL1=3  JL2=2    OVER=-7.180D-12
  L=2  JL1=3  JL2=3    OVER= 1.000D+00
  L=2  JL1=4  JL2=1    OVER= 9.422D-18
  L=2  JL1=4  JL2=2    OVER= 3.052D-18
  L=2  JL1=4  JL2=3    OVER= 3.722D-18
  L=2  JL1=4  JL2=4    OVER= 1.000D+00
  L=2  JL1=5  JL2=1    OVER= 5.734D-18
  L=2  JL1=5  JL2=2    OVER= 1.267D-17
  L=2  JL1=5  JL2=3    OVER= 1.802D-17
  L=2  JL1=5  JL2=4    OVER=-3.590D-12
  L=2  JL1=5  JL2=5    OVER= 1.000D+00
 
  L=3  JL1=1  JL2=1    OVER= 1.000D+00
  L=3  JL1=2  JL2=1    OVER= 4.403D-18
  L=3  JL1=2  JL2=2    OVER= 1.000D+00
  L=3  JL1=3  JL2=1    OVER= 9.687D-18
  L=3  JL1=3  JL2=2    OVER= 6.282D-12
  L=3  JL1=3  JL2=3    OVER= 1.000D+00
  L=3  JL1=4  JL2=1    OVER=-8.310D-19
  L=3  JL1=4  JL2=2    OVER=-1.393D-10
  L=3  JL1=4  JL2=3    OVER=-4.866D-12
  L=3  JL1=4  JL2=4    OVER= 1.000D+00
  L=3  JL1=5  JL2=1    OVER= 2.234D-18
  L=3  JL1=5  JL2=2    OVER= 8.574D-12
  L=3  JL1=5  JL2=3    OVER=-8.268D-12
  L=3  JL1=5  JL2=4    OVER=-8.974D-13
  L=3  JL1=5  JL2=5    OVER= 1.000D+00
  L=3  JL1=6  JL2=1    OVER= 8.851D-10
  L=3  JL1=6  JL2=2    OVER=-1.494D-17
  L=3  JL1=6  JL2=3    OVER=-2.581D-18
  L=3  JL1=6  JL2=4    OVER= 1.187D-17
  L=3  JL1=6  JL2=5    OVER= 1.778D-17
  L=3  JL1=6  JL2=6    OVER= 1.000D+00
  L=3  JL1=7  JL2=1    OVER= 2.781D-12
  L=3  JL1=7  JL2=2    OVER= 3.244D-18
  L=3  JL1=7  JL2=3    OVER=-1.296D-18
  L=3  JL1=7  JL2=4    OVER=-4.245D-18
  L=3  JL1=7  JL2=5    OVER=-7.405D-18
  L=3  JL1=7  JL2=6    OVER=-7.180D-12
  L=3  JL1=7  JL2=7    OVER= 1.000D+00




	IMPLICIT REAL*8(A-H,O-Z), INTEGER*4(I-N)
      COMMON /CUBICY/Y0(0:3),NY0(0:3,7),TY0(0:3,7,7),NEXP(0:3,7,7,7)
      DATA PI /3.14159265D0/

      IW=7
      NT0=1000
	NF0=2000
      DO 10 L=0,3
      DO 11 JL1=1,2*L+1
      DO 12 JL2=1,JL1

      SUM1=0.D0
      DO 1 NT=1,NT0
      T=DFLOAT(NT)*PI/DFLOAT(NT0)

            SUM2=0.D0
            DO 2 NF=1,NF0
            F=DFLOAT(NF)*PI*2.D0/DFLOAT(NF0)

            X=DSIN(T)*DCOS(F)
	      Y=DSIN(T)*DSIN(F)
	      Z=DCOS(T)

               YL1=0.D0
               DO  K=1,NY0(L,JL1)
               YL1= YL1+Y0(L)*TY0(L,JL1,k)*
     |                    x**NEXP(L,JL1,k,1)*
     |                    y**NEXP(L,JL1,k,2)*
     |                    z**NEXP(L,JL1,k,3)
               ENDDO
               YL2=0.D0
               DO  K=1,NY0(L,JL2)
               YL2= YL2+Y0(L)*TY0(L,JL2,k)*
     |                    x**NEXP(L,JL2,k,1)*
     |                    y**NEXP(L,JL2,k,2)*
     |                    z**NEXP(L,JL2,k,3)
               ENDDO
    2        SUM2=SUM2+YL1*YL2*2.D0*PI/DFLOAT(NF0)
    1 SUM1=SUM1 + SUM2*DSIN(T)*PI/DFLOAT(NT0)
   
      WRITE(IW,600) L,JL1,JL2,SUM1
  600 FORMAT('  L=',I1,'  JL1=',I1,'  JL2=',I1,3X,' OVER=', 1PD10.3)
   12 CONTINUE

   11 CONTINUE
      WRITE(IW,*)
   10 CONTINUE
      STOP
	END

C =======================================================================
C =======================================================================
C =======================================================================
C =======================================================================
c structure of the cubic harmonics (=Real Spherical harmonics in 
c                                         carstessian coordinates)

c  CUBIC HARMONICS_[L,JL](x,y,z)= Y0(L)*Sum[  TY0(L,JL,k)*
c          x^NEXP(L,JL,k,1)*
C          y^NEXP(L,JL,k,2)*
C          z^NEXP(L,JL,k,3),  {k, 1, NY0[L,JL]}]

C  JL=1, 2,.....2*L+1
C   L=0, 1, 2, 3 (S, P, D & F)
C miraglia 15/03/2011
C-----------------------------------------------------------------------
c
      BLOCKDATA CUBICYLM
	IMPLICIT REAL*8(A-H,O-Z), INTEGER*4(I-N)
      COMMON /CUBICY/Y0(0:3),NY0(0:3,7),TY0(0:3,7,7),NEXP(0:3,7,7,7)
c---------------------------------------------------------
c-----s state---------------------------------------------
      DATA Y0(0) /0.2820947917738781D0/ !1/Sqrt[4*Pi]
c     total of 2*0+1=1 states
c.........................................................
	DATA  NY0(0,1) /1/    ! 1 summands

      DATA  TY0(0,1,1)/ 1.D0/ !1st term
      DATA NEXP(0,1,1,1)/0/ !1st term power of x
      DATA NEXP(0,1,1,2)/0/ !                  y
      DATA NEXP(0,1,1,3)/0/ !                  z

c---------------------------------------------------------
c-----P states -------------------------------------------
      DATA Y0(1) /0.4886025119029199D0/ !Sqrt[3/(4*Pi)]
c     total of 2*1+1=3 states
c.........................................................
 
	DATA  NY0(1,1) /1/    ! 1 summands

      DATA  TY0(1,1,1)/ 1.D0/ !1st term
      DATA NEXP(1,1,1,1)/1/ !1st term power of x
      DATA NEXP(1,1,1,2)/0/ !                  y
      DATA NEXP(1,1,1,3)/0/ !                  z

c.........................................................
 
	DATA  NY0(1,2) /1/    ! 1 summands

      DATA  TY0(1,2,1)/ 1.D0/ !1st term
      DATA NEXP(1,2,1,1)/0/ !1st term power of x
      DATA NEXP(1,2,1,2)/1/ !                  y
      DATA NEXP(1,2,1,3)/0/ !                  z

c.........................................................
 
	DATA  NY0(1,3) /1/    ! 1 summands

      DATA  TY0(1,3,1)/ 1.D0/ !1st term
      DATA NEXP(1,3,1,1)/0/ !1st term power of x
      DATA NEXP(1,3,1,2)/0/ !                  y
      DATA NEXP(1,3,1,3)/1/ !                  z

c---------------------------------------------------------
c----D states---------------------------------------------
      DATA   Y0(2) /0.3153915652525200D0/ !Sqrt[5/(16*Pi)]
c     total of 2*2+1=5 states
c.........................................................
 
	DATA  NY0(2,1) /3/    ! 3 summands

      DATA  TY0(2,1,1)/-1.D0/ !1st term
      DATA  TY0(2,1,2)/-1.D0/ !2nd term
      DATA  TY0(2,1,3)/ 2.D0/ !3rd term

      DATA NEXP(2,1,1,1)/2/ !1st term power of x
      DATA NEXP(2,1,1,2)/0/ !                  y
      DATA NEXP(2,1,1,3)/0/ !                  z

      DATA NEXP(2,1,2,1)/0/ !2nd               x
      DATA NEXP(2,1,2,2)/2/ !                  y
      DATA NEXP(2,1,2,3)/0/ !                  z

      DATA NEXP(2,1,3,1)/0/ !3rd               x
      DATA NEXP(2,1,3,2)/0/ !                  y
      DATA NEXP(2,1,3,3)/2/ !                  z
c......................................................... 

	DATA  NY0(2,2) /2/    ! 2 summands
      DATA  TY0(2,2,1)/ 1.732050807568877D0/ ! Sqrt[3] !1st term
      DATA  TY0(2,2,2)/-1.732050807568877D0/ !-Sqrt[3] !2nd term
      DATA NEXP(2,2,1,1)/2/ !1st term power of x
      DATA NEXP(2,2,1,2)/0/ !                  y
      DATA NEXP(2,2,1,3)/0/ !                  z

      DATA NEXP(2,2,2,1)/0/ !2nd               x
      DATA NEXP(2,2,2,2)/2/ !                  y
      DATA NEXP(2,2,2,3)/0/ !                  z
c......................................................... 

	DATA  NY0(2,3) /1/    ! 1 summand
      DATA  TY0(2,3,1)/ 3.464101615137755D0/ ! 2*Sqrt[3] !1st term
      DATA NEXP(2,3,1,1)/1/ !1st term power of x
      DATA NEXP(2,3,1,2)/1/ !                  y
      DATA NEXP(2,3,1,3)/0/ !                  z
c......................................................... 

	DATA  NY0(2,4) /1/  ! 1 summand
      DATA  TY0(2,4,1)/ 3.464101615137755D0/ ! 2*Sqrt[3]
      DATA NEXP(2,4,1,1)/0/ !1st term power of x
      DATA NEXP(2,4,1,2)/1/ !                  y
      DATA NEXP(2,4,1,3)/1/ !                  z
c......................................................... 

	DATA  NY0(2,5) /1/    ! 1 summand
      DATA  TY0(2,5,1)/ 3.464101615137755D0/ ! 2*Sqrt[3]
      DATA NEXP(2,5,1,1)/1/ !1st term power of x
      DATA NEXP(2,5,1,2)/0/ !                  y
      DATA NEXP(2,5,1,3)/1/ !                  z

c---------------------------------------------------------
c----F states---------------------------------------------
      DATA   Y0(3) /0.3731763325901154D0/ !Sqrt[7/(16*Pi)]
c     total of 2*3+1=7 states
c......................................................... 

	DATA  NY0(3,1) /3/    ! 3 summands

      DATA  TY0(3,1,1)/-3.D0/ !1st term
      DATA  TY0(3,1,2)/-3.D0/ !2nd term
      DATA  TY0(3,1,3)/ 2.D0/ !3rd term

      DATA NEXP(3,1,1,1)/2/ !1st term power of x
      DATA NEXP(3,1,1,2)/0/ !                  y
      DATA NEXP(3,1,1,3)/1/ !                  z

      DATA NEXP(3,1,2,1)/0/ !2nd               x
      DATA NEXP(3,1,2,2)/2/ !                  y
      DATA NEXP(3,1,2,3)/1/ !                  z

      DATA NEXP(3,1,3,1)/0/ !3rd               x
      DATA NEXP(3,1,3,2)/0/ !                  y
      DATA NEXP(3,1,3,3)/3/ !                  z
c......................................................... 

	DATA  NY0(3,2) /3/    ! 3 summands

      DATA  TY0(3,2,1)/ 2.D0/ !1st term
      DATA  TY0(3,2,2)/-3.D0/ !2nd term
      DATA  TY0(3,2,3)/-3.D0/ !3rd term

      DATA NEXP(3,2,1,1)/3/ !1st term power of x
      DATA NEXP(3,2,1,2)/0/ !                  y
      DATA NEXP(3,2,1,3)/0/ !                  z

      DATA NEXP(3,2,2,1)/1/ !2nd               x
      DATA NEXP(3,2,2,2)/2/ !                  y
      DATA NEXP(3,2,2,3)/0/ !                  z

      DATA NEXP(3,2,3,1)/1/ !3rd               x
      DATA NEXP(3,2,3,2)/0/ !                  y
      DATA NEXP(3,2,3,3)/2/ !                  z
 
c......................................................... 

	DATA  NY0(3,3) /3/    ! 3 summands

      DATA  TY0(3,3,1)/-3.D0/ !1st term
      DATA  TY0(3,3,2)/ 2.D0/ !2nd term
      DATA  TY0(3,3,3)/-3.D0/ !3rd term

      DATA NEXP(3,3,1,1)/2/ !1st term power of x
      DATA NEXP(3,3,1,2)/1/ !                  y
      DATA NEXP(3,3,1,3)/0/ !                  z

      DATA NEXP(3,3,2,1)/0/ !2nd               x
      DATA NEXP(3,3,2,2)/3/ !                  y
      DATA NEXP(3,3,2,3)/0/ !                  z

      DATA NEXP(3,3,3,1)/0/ !3rd               x
      DATA NEXP(3,3,3,2)/1/ !                  y
      DATA NEXP(3,3,3,3)/2/ !                  z

c......................................................... 

	DATA  NY0(3,4) /2/    ! 2 summands

      DATA  TY0(3,4,1)/ 3.872983346207417D0/ !1st term !Sqrt[15]
      DATA  TY0(3,4,2)/-3.872983346207417D0/ !2nd term !-Sqrt[15]

      DATA NEXP(3,4,1,1)/1/ !1st term power of x
      DATA NEXP(3,4,1,2)/2/ !                  y
      DATA NEXP(3,4,1,3)/0/ !                  z

      DATA NEXP(3,4,2,1)/1/ !2nd               x
      DATA NEXP(3,4,2,2)/0/ !                  y
      DATA NEXP(3,4,2,3)/2/ !                  z

c......................................................... 

	DATA  NY0(3,5) /2/    ! 2 summands

      DATA  TY0(3,5,1)/ 3.872983346207417D0/ !1st term !Sqrt[15]
      DATA  TY0(3,5,2)/-3.872983346207417D0/ !2nd term !-Sqrt[15]

      DATA NEXP(3,5,1,1)/0/ !1st term power of x
      DATA NEXP(3,5,1,2)/1/ !                  y
      DATA NEXP(3,5,1,3)/2/ !                  z

      DATA NEXP(3,5,2,1)/2/ !2nd               x
      DATA NEXP(3,5,2,2)/1/ !                  y
      DATA NEXP(3,5,2,3)/0/ !                  z

c......................................................... 

	DATA  NY0(3,6) /2/    ! 2 summands

      DATA  TY0(3,6,1)/ 3.872983346207417D0/ !1st term !Sqrt[15]
      DATA  TY0(3,6,2)/-3.872983346207417D0/ !2nd term !-Sqrt[15]

      DATA NEXP(3,6,1,1)/2/ !1st term power of x
      DATA NEXP(3,6,1,2)/0/ !                  y
      DATA NEXP(3,6,1,3)/1/ !                  z

      DATA NEXP(3,6,2,1)/0/ !2nd               x
      DATA NEXP(3,6,2,2)/2/ !                  y
      DATA NEXP(3,6,2,3)/1/ !                  z

c......................................................... 

	DATA  NY0(3,7) /1/    ! 1 summand1

      DATA  TY0(3,7,1)/ 7.745966692414834D0/ !1st term !2*Sqrt[15]

      DATA NEXP(3,7,1,1)/1/ !1st term power of x
      DATA NEXP(3,7,1,2)/1/ !                  y
      DATA NEXP(3,7,1,3)/1/ !                  z

c---------------------------------------------------------
c---------------------------------------------------------
      END