o1 = a**2
        o2 = -2.00000000000000d0*mass*rBL
        o3 = rBL**2
        o4 = o1 + o2 + o3
        o5 = sqrt(o4)
        o6 = cos(qgrd(j))
        o7 = o6**2
        o8 = o1*o7
        o9 = o3 + o8
        o10 = o9**2
        o11 = 1/o10
        o12 = 1.0!psi2dv(eta,q)
        o13 = o12**2
        o14 = 1/o13
        o15 = sin(qgrd(j))
        o16 = o15**2
        o17 = o1 + o3
        o18 = 2.00000000000000d0*o17*o3
        o19 = -o1
        o20 = o19 + o3
        o21 = o20*o9
        o22 = o18 + o21
        o23 = mass**2
        o24 = o10**2
        o25 = 1/o24
        o26 = o22**2
        o27 = o13**2
        o28 = o27**2
        o29 = o27*o28
        o30 = 1/o29
        o31 = o1**2
        o32 = o16**2
c        exc33(i,j) = 0
c        exc32(i,j) = -2.00000000000000d0*a*mass*o1*o11*o14*o15*o16*o5*o6
c     &  *rBL
c        exc31(i,j) = a*mass*o11*o14*o16*o22
c        exc22(i,j) = 0
c        exc21(i,j) = 0
c        exc11(i,j) = 0
        ksq(i,j) = 2.00000000000000d0*o1*o16*o23*o25*o26*o30 + 8.0000000
     &  000000d0*o1*o23*o25*o3*o30*o31*o32*o4*o7