subroutine dipd(dipc,r1,r2,xcos,nu)
 !transform generalised coordinates to those for particular c
 !system. this version transforms to ab2 bondlength-bondangle c
 !coordinates. allowance must be made for the numbering of the atoms c
 !additionally, the zbisc option is included. c
      implicit double precision (a-h,o-y), logical (z)
      !
      double precision,intent(out) :: dipc
      double precision,intent(in)  :: r1,r2,xcos
      integer,intent(in)           :: nu
      double precision             :: dipx,dipy
      !
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      !
      !  (r = r . s = r'. t = theta) c
      data x1/1.0d0/,x0/0.0d0/,tiny/9.0d-15/,x2/2.0d0/,pi/3.14159265358979323846264d0/
      !
      if (g1 .eq. x0) then
         !  bondlength bondangle coordinates: atom 1 = atom 2
         q1 = r1
         q2 = r2
         theta = acos(xcos)
         cost = xcos
         !
      else if (g2 .eq. x0) then
         !  scattering coordinates: atom 2 = atom 3
         xx = r1 * g1
         yy = r1 * (x1 - g1)
         alpha= acos(xcos)
         if (r2 .eq. x0 .or. xcos .ge. (x1 - tiny)) then
            q1 = abs(xx - r2)
            q2 = (yy + r2)
            cost = -x1
         else if (xcos .le. (tiny - x1)) then
            q1 = (xx + r2)
            q2 = abs(yy + r2)
            cost = x1
         else
            q1 = sqrt(xx*xx + r2*r2 - x2*xx*r2*xcos)
            q2 = sqrt(yy*yy + r2*r2 + x2*yy*r2*xcos)
            cost = (q1**2 + q2**2 - r1**2) / (x2 * q1 * q2)
         endif
         xsin= sqrt(1.0d0 - xcos*xcos)
         theta = acos(cost)
         beta= asin(xsin*yy/q2)
      else
         !  general coordinates (including radau): atom 1 = atom 2
         f1= x1/g1
         f2= x1/g2
         f12= x1 - f1*f2
         p1= r1*(x1-f1)/(g2*f12)
         p2= r2*(x1-f2)/(g1*f12)
         s1= r1-p1
         s2= r2-p2
         q1= sqrt(p1*p1 + s2*s2 + x2*p1*s2*xcos)/(x1-g1)
         q2= sqrt(p2*p2 + s1*s1 + x2*p2*s1*xcos)/(x1-g2)
         q3= sqrt(p1*p1 + p2*p2 - x2*p1*p2*xcos)
         cost = (q1*q1 + q2*q2 - q3*q3)/(x2*q1*q2)
         theta = acos(cost)
      endif
      !
      !call dms_h2s(dipx,dipy,q1,q2,theta)
      !
      call dms_h2s_qz(dipx,dipy,q1,q2,theta)
      !
      !   bondlength-bondangle co-ordinates c
      !
      if (g1.eq.x0) then
         gamma= theta/x2
         ycos= cos(gamma)
         ysin= sin(gamma)
         if (zembed) then
            if (nu.eq.0) then
               dipc= +dipy*ycos - dipx*ysin
            else
               dipc= +dipx*ycos + dipy*ysin
            endif
         else
            if (nu.eq.0) then
               dipc= +dipy*ycos + dipx*ysin
            else
               dipc= -dipx*ycos + dipy*ysin
            endif
         endif
         !
         !  scattering co-ordinates c
      else if (g2.eq.x0) then
         !
         gamma= beta - theta/x2
         if (zembed) then
            ycos= cos(gamma)
            ysin= sin(gamma)
            if (nu.eq.0) then
               dipc= -dipx*ysin - dipy*ycos
            else
               dipc= +dipx*ycos - dipy*ysin
            endif
         else
            delta= alpha - gamma
            ycos= cos(delta)
            ysin= sin(delta)
            if (nu.eq.0) then
               dipc= +dipx*ysin + dipy*ycos
            else
               dipc= +dipx*ycos - dipy*ysin
            endif
         endif
         !
         !  all other co-ordinates c
      else
         !
         if (zbisc) then
            h1 = g1*q1
            cosa = (r2*r2 + q2*q2 - h1*h1)/(x2*r2*q2)
            alpha = (acos(xcos)-theta)/x2 - acos(cosa)
            ycos= - cos(alpha)
            ysin= + sin(alpha)
            if (nu.eq.1) then
               dipc= +dipx*ysin - dipy*ycos
            else
               dipc= -dipx*ycos - dipy*ysin
            endif
         elseif (zembed) then
            h1 = g1*q1
            cosa = (r2*r2 + h1*h1 - q2*q2)/(x2*r2*h1)
            alpha= acos(cosa)
            ycos= - cos(alpha + theta/x2)
            ysin= + sin(alpha + theta/x2)
            if (nu.eq.0) then
               dipc= -dipx*ysin + dipy*ycos
            else
               dipc= +dipx*ycos + dipy*ysin
            endif
         else
            h2 = g2*q2
            cosa = (r1*r1 + h2*h2 - q1*q1)/(x2*r1*h2)
            alpha = acos(cosa)
            ycos= - cos(alpha + theta/x2)
            ysin= + sin(alpha + theta/x2)
            if (nu.eq.0) then
               dipc= +dipx*ysin + dipy*ycos
            else
               dipc= -dipx*ycos + dipy*ysin
            endif
         endif
      endif
      return
      !
  end subroutine dipd


!     NOTES:
!     -----
!     CCSD[T] DMS of H2S copmuted by Cours, Rosmus, and Tyuterev, JCP 117, 5192 (2002), 
!     referred to as set II in this paper.
!     S is at the origin, the H2 atom lies in the positive P-Q quadrant.  The cartesian coordinate
!     system is such that the Q-axis bisects the included bond angle for
!     any geometry.  This subroutine is intended to be called by the
!     DIPS routine in the triatom program suite of Tennyson et al. (1993).
!
!
!                                ^  P
!                                l        H2
!                                l
!                                l
!                - - - - - - - - S - - - - - - - - >  Q
!                                l
!                                l
!                                l         H1

  subroutine dms_h2s(dipp,dipq,q1,q2,alpha)
   !
   double precision,intent(in)  :: q1,q2,alpha
   double precision,intent(out) :: dipp,dipq
   !
   double precision :: r1,r2,req,alphaeq
   double precision :: f(3),y(3)

   double precision,parameter :: pi=3.14159265358979323846264d0,bohr = 0.529177249d0,autode = 2.5417662D0
   !
   integer  :: iq(3,37),ip(3,26),i
   double precision :: p(26),q(37)
   
   ip = transpose(reshape((/ &
            1,1,1,1,0,0,0,0,2,2,2,0,0,0,3,3,0,0,2,2,1,1,4,3,0,1,&
            0,0,0,0,1,1,1,1,0,0,0,2,2,2,0,0,3,3,1,1,2,2,0,1,4,3,&
            0,1,2,3,0,1,2,3,0,1,2,0,1,2,0,1,0,1,0,1,0,1,0,0,0,0/),(/26,3/)))
   
   iq = transpose(reshape((/ &
             0,0,0,0,0,0,0,1,1,1,1,2,2,2,1,1,1,3,3,2,2,4,3,2,0,0,0,0,0,0,0,0,0,1,1,0,1,&
             0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,0,1,2,1,1,1,1,2,2,2,3,3,2,2,4,3,&
             0,1,2,3,4,5,6,0,1,2,3,0,1,2,0,1,2,0,1,0,1,0,0,0,0,1,2,3,0,1,2,0,1,0,1,0,0/),(/37,3/)))
    !
    p = (/-0.0011949d0,0.7494089d0,-0.2577043d0,0.1158388d0,0.0011949d0,-0.7494089d0,0.2577043d0,-0.1158388d0,0.1744552d0,-0.3865805d0,0.3602985d0,-0.1744552d0,0.3865805d0,-0.3602985d0,0.0871393d0,-0.5725051d0,-0.0871393d0,0.5725051d0,0.1968289d0,0.3250446d0,-0.1968289d0,-0.3250446d0,0.0418678d0,0.0803807d0,-0.0418678d0,-0.0803807d0/)
    q = (/0.9616806d0,-0.1584264d0,0.6139368d0,0.1716701d0,0.7225179d0,0.2645379d0,0.0000000d0,-0.0173490d0,0.4847617d0,-0.1531882d0,-0.0107459d0,-0.1192017d0,0.0109174d0,-1.3908281d0,-0.2454039d0,0.7981090d0,1.0311684d0,-0.1130682d0,-0.8015431d0,0.2283408d0,-3.0729349d0,0.1638799d0,-0.0684265d0,-0.7991572d0,-0.0173490d0,0.4847617d0,-0.1531882d0,-0.0107459d0,-0.1192017d0,0.0109174d0,-1.3908281d0,-0.1130682d0,-0.8015431d0,0.2283408d0,-3.0729349d0,0.1638799d0,-0.0684265d0/)
    !
    r1 = q1*bohr
    r2 = q2*bohr
    !
    req = 1.336429d0
    alphaeq = 92.196125D0*pi/180.0d0
    !
    y(1) = r1 - req
    y(2) = r2 - req
    y(3) = cos(alpha) - cos(alphaeq)
    !
    dipp = 0 
    !
    do i=1,size(p)
       !
       dipp = dipp + y(1)**ip(1,i)*y(2)**ip(2,i)*y(3)**ip(3,i)*p(i)
       !
       !if (ip(1,i)/=ip(2,i)) then
       !  !
       !  dipp = dipp - y(1)**ip(2,i)*y(2)**ip(1,i)*y(3)**ip(3,i)*p(i)
       !  !
       !endif
       !
    enddo
    !
    dipq = 0 
    !
    do i=1,size(q)
       !
       dipq = dipq + y(1)**iq(1,i)*y(2)**iq(2,i)*y(3)**iq(3,i)*q(i)
       !
       !if (iq(1,i)/=iq(2,i)) then
       !  !
       !  dipq = dipq + y(1)**iq(2,i)*y(2)**iq(1,i)*y(3)**iq(3,i)*q(i)
       !  !
       !endif
       !
    enddo
    !
    dipq = dipq * sin(alpha)
    !
    dipp = -dipp/autode   
    dipq =  dipq/autode
    !
  end subroutine dms_h2s
  !
!     NOTES:
!     -----
!     CCSD(T)/aug-cc-pV(Q+d)Z DMS of H2S copmuted by Yurchenko, 2011, 
!     S is at the origin, the H2 atom lies in the positive P-Q quadrant.  The cartesian coordinate
!     system is such that the Q-axis bisects the included bond angle for
!     any geometry.  This subroutine is intended to be called by the
!     DIPS routine in the triatom program suite of Tennyson et al. (1993).
  !
  subroutine dms_h2s_qz(dipp,dipq,q1,q2,alpha)
   !
   double precision,intent(in)  :: q1,q2,alpha
   double precision,intent(out) :: dipp,dipq
   !
   double precision             :: r1,r2
   double precision,parameter   :: bohr = 0.529177249d0,autode = 2.5417662D0
   !
   r1 = q1*bohr
   r2 = q2*bohr
   !
   call dipol_xy2_p_qz(r1,r2,alpha,dipp)
   call dipol_xy2_q_qz(r1,r2,alpha,dipq)
   !
   dipp =-dipp/autode  
   dipq = dipq/autode
   !
  end subroutine dms_h2s_qz
  !
  subroutine dipol_xy2_p_qz(r1,r2,alpha,f)
      !
      double precision,intent(in):: r1,r2,alpha
      integer,parameter          :: N = 72
      double precision           :: p(3:N),f
      double precision           :: y1,y2,y3
      double precision           :: ae,re,v0,v1,v2,v3,v4,v5,v6,v7,v8
      double precision,parameter :: pi=3.14159265358979323846264d0

!CCSD(T)/aug-cc-pV(6+d)Z,after adding the corrections,dump=1 (the complete surface up to 10000cm-1) 
p=(/0.00478832298768,-0.76979371155700,-0.235102597053,0.221487070349,0.392103566418,-0.235313238562,-0.201920584750,-0.04406670432450,-0.08155801815960,-0.22732589567600,-0.42901497958100,-0.08342893145700,0.32742586865100,0.35835788269200,0.01101719774270,0.86537852828000,-0.74141092287700,-0.31440493326500,0.00000000000000,0.84582980447000,-0.34837144518800,0.47161811756700,-0.46659525419400,-0.45517814081000,0.85332644680100,0.34423282322100,1.45102983805000,1.07338903107000,-1.75612011876000,-3.34559352047000,0.00000000000000,2.02849627363000,1.19866670605000,-9.26719338566000,0.00000000000000,0.00000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000,0.0000000000000/)
p=-p
b0     = 1.0d0

      re     = 1.336d0
      ae     = 92.2d0/180.0d0*pi
      !
      y1 = (r1 - re)*exp(-b0*(r1-re)**2)
      y2 = (r2 - re)*exp(-b0*(r2-re)**2)
      y3 = cos(alpha) - cos(ae)
      !
 v1 = p( 3)*y1**1*y2**0*y3**0& 
    - p( 3)*y1**0*y2**1*y3**0
 v2 = p( 4)*y1**1*y2**0*y3**1& 
    - p( 4)*y1**0*y2**1*y3**1&  
    + p( 5)*y1**2*y2**0*y3**0& 
    - p( 5)*y1**0*y2**2*y3**0

 v3 = p( 6)*y1**1*y2**0*y3**2& 
    - p( 6)*y1**0*y2**1*y3**2& 
    + p( 7)*y1**2*y2**0*y3**1& 
    - p( 7)*y1**0*y2**2*y3**1& 
    + p( 8)*y1**2*y2**1*y3**0& 
    - p( 8)*y1**1*y2**2*y3**0& 
    + p( 9)*y1**3*y2**0*y3**0& 
    - p( 9)*y1**0*y2**3*y3**0

  v4 =p(10)*y1**1*y2**0*y3**3& 
    - p(10)*y1**0*y2**1*y3**3& 
    + p(11)*y1**2*y2**0*y3**2& 
    - p(11)*y1**0*y2**2*y3**2& 
    + p(12)*y1**2*y2**1*y3**1& 
    - p(12)*y1**1*y2**2*y3**1& 
    + p(13)*y1**3*y2**0*y3**1& 
    - p(13)*y1**0*y2**3*y3**1& 
    + p(14)*y1**3*y2**1*y3**0& 
    - p(14)*y1**1*y2**3*y3**0& 
    + p(15)*y1**4*y2**0*y3**0& 
    - p(15)*y1**0*y2**4*y3**0

  v5 =p(16)*y1**1*y2**0*y3**4& 
    - p(16)*y1**0*y2**1*y3**4& 
    + p(17)*y1**2*y2**0*y3**3& 
    - p(17)*y1**0*y2**2*y3**3& 
    + p(18)*y1**2*y2**1*y3**2& 
    - p(18)*y1**1*y2**2*y3**2& 
    + p(19)*y1**3*y2**0*y3**2& 
    - p(19)*y1**0*y2**3*y3**2& 
    + p(20)*y1**3*y2**1*y3**1& 
    - p(20)*y1**1*y2**3*y3**1& 
    + p(21)*y1**3*y2**2*y3**0& 
    - p(21)*y1**2*y2**3*y3**0& 
    + p(22)*y1**4*y2**0*y3**1& 
    - p(22)*y1**0*y2**4*y3**1& 
    + p(23)*y1**4*y2**1*y3**0& 
    - p(23)*y1**1*y2**4*y3**0& 
    + p(24)*y1**5*y2**0*y3**0& 
    - p(24)*y1**0*y2**5*y3**0

  v6 =p(25)*y1**1*y2**0*y3**5& 
    - p(25)*y1**0*y2**1*y3**5& 
    + p(26)*y1**2*y2**0*y3**4& 
    - p(26)*y1**0*y2**2*y3**4& 
    + p(27)*y1**2*y2**1*y3**3& 
    - p(27)*y1**1*y2**2*y3**3& 
    + p(28)*y1**3*y2**0*y3**3& 
    - p(28)*y1**0*y2**3*y3**3& 
    + p(29)*y1**3*y2**1*y3**2& 
    - p(29)*y1**1*y2**3*y3**2& 
    + p(30)*y1**3*y2**2*y3**1& 
    - p(30)*y1**2*y2**3*y3**1& 
    + p(31)*y1**4*y2**0*y3**2& 
    - p(31)*y1**0*y2**4*y3**2& 
    + p(32)*y1**4*y2**1*y3**1& 
    - p(32)*y1**1*y2**4*y3**1& 
    + p(33)*y1**4*y2**2*y3**0& 
    - p(33)*y1**2*y2**4*y3**0& 
    + p(34)*y1**5*y2**0*y3**1& 
    - p(34)*y1**0*y2**5*y3**1& 
    + p(35)*y1**5*y2**1*y3**0& 
    - p(35)*y1**1*y2**5*y3**0& 
    + p(36)*y1**6*y2**0*y3**0& 
    - p(36)*y1**0*y2**6*y3**0

 v7 = p(37)*y1**1*y2**0*y3**6& 
    - p(37)*y1**0*y2**1*y3**6& 
    + p(38)*y1**2*y2**0*y3**5& 
    - p(38)*y1**0*y2**2*y3**5& 
    + p(39)*y1**2*y2**1*y3**4& 
    - p(39)*y1**1*y2**2*y3**4& 
    + p(40)*y1**3*y2**0*y3**4& 
    - p(40)*y1**0*y2**3*y3**4& 
    + p(41)*y1**3*y2**1*y3**3& 
    - p(41)*y1**1*y2**3*y3**3& 
    + p(42)*y1**3*y2**2*y3**2& 
    - p(42)*y1**2*y2**3*y3**2& 
    + p(43)*y1**4*y2**0*y3**3& 
    - p(43)*y1**0*y2**4*y3**3& 
    + p(44)*y1**4*y2**1*y3**2& 
    - p(44)*y1**1*y2**4*y3**2& 
    + p(45)*y1**4*y2**2*y3**1& 
    - p(45)*y1**2*y2**4*y3**1& 
    + p(46)*y1**4*y2**3*y3**0& 
    - p(46)*y1**3*y2**4*y3**0& 
    + p(47)*y1**5*y2**0*y3**2& 
    - p(47)*y1**0*y2**5*y3**2& 
    + p(48)*y1**5*y2**1*y3**1& 
    - p(48)*y1**1*y2**5*y3**1& 
    + p(49)*y1**5*y2**2*y3**0& 
    - p(49)*y1**2*y2**5*y3**0& 
    + p(50)*y1**6*y2**0*y3**1& 
    - p(50)*y1**0*y2**6*y3**1& 
    + p(51)*y1**6*y2**1*y3**0& 
    - p(51)*y1**1*y2**6*y3**0& 
    + p(52)*y1**7*y2**0*y3**0& 
    - p(52)*y1**0*y2**7*y3**0

 v8 = p(53)*y1**1*y2**0*y3**7& 
    - p(53)*y1**0*y2**1*y3**7& 
    + p(54)*y1**2*y2**0*y3**6& 
    - p(54)*y1**0*y2**2*y3**6& 
    + p(55)*y1**2*y2**1*y3**5& 
    - p(55)*y1**1*y2**2*y3**5& 
    + p(56)*y1**3*y2**0*y3**5& 
    - p(56)*y1**0*y2**3*y3**5& 
    + p(57)*y1**3*y2**1*y3**4& 
    - p(57)*y1**1*y2**3*y3**4& 
    + p(58)*y1**3*y2**2*y3**3& 
    - p(58)*y1**2*y2**3*y3**3& 
    + p(59)*y1**4*y2**0*y3**4& 
    - p(59)*y1**0*y2**4*y3**4& 
    + p(60)*y1**4*y2**1*y3**3& 
    - p(60)*y1**1*y2**4*y3**3& 
    + p(61)*y1**4*y2**2*y3**2& 
    - p(61)*y1**2*y2**4*y3**2& 
    + p(62)*y1**4*y2**3*y3**1& 
    - p(62)*y1**3*y2**4*y3**1& 
    + p(63)*y1**5*y2**0*y3**3& 
    - p(63)*y1**0*y2**5*y3**3& 
    + p(64)*y1**5*y2**1*y3**2& 
    - p(64)*y1**1*y2**5*y3**2& 
    + p(65)*y1**5*y2**2*y3**1& 
    - p(65)*y1**2*y2**5*y3**1& 
    + p(66)*y1**5*y2**3*y3**0& 
    - p(66)*y1**3*y2**5*y3**0& 
    + p(67)*y1**6*y2**0*y3**2& 
    - p(67)*y1**0*y2**6*y3**2& 
    + p(68)*y1**6*y2**1*y3**1& 
    - p(68)*y1**1*y2**6*y3**1& 
    + p(69)*y1**6*y2**2*y3**0& 
    - p(69)*y1**2*y2**6*y3**0& 
    + p(70)*y1**7*y2**0*y3**1& 
    - p(70)*y1**0*y2**7*y3**1& 
    + p(71)*y1**7*y2**1*y3**0& 
    - p(71)*y1**1*y2**7*y3**0& 
    + p(72)*y1**8*y2**0*y3**0& 
    - p(72)*y1**0*y2**8*y3**0

    f=v1+v2+v3+v4+v5+v6+v7+v8
    !
  end subroutine dipol_xy2_p_qz

  !
  subroutine dipol_xy2_q_qz(r1,r2,alpha,f)
      !
      double precision,intent(in) :: r1,r2,alpha
      integer,parameter :: N = 99
      double precision :: q(5:N),f
      double precision            :: y1,y2,y3
      double precision            :: ae,re,v0,v1,v2,v3,v4,v5,v6,v7,v8
      double precision,parameter :: pi=3.14159265358979323846264d0
      !

!CCSD(T)/aug-cc-pV(6+d)Z,after adding the corrections,dump=1, (the complete surface up to 10000cm-1)
q=(/0.971586586978,-0.154669269244,-0.0103547035951,0.631678225719,0.468306603182,-0.199744184765,-0.160945068349,0.14924506451,-0.193136744056,0.253560036277,0.0749427103131,0.097845862367,-0.247322531599,0.497899435495,0.139727431856,0.000000000000,-0.296478903793,0.0375506517156,0.11894339122,0.150249679713,-0.380932782156,0.398238925603,0.419436974602,-0.079723123993,0.844421821396,-0.181669708916,0.395662875017,0.539606808952,-0.713817995961,-0.773713256945,-0.118918435926,-0.119030652505,0.522999467924,0.144946313652,0.835619288737,-0.110287239252,-0.993936795627,-0.813442843351,-0.267447777307,1.39085188424,0.194713595783,0.367885601675,0.75201022025,1.78047067626,-2.35559011703,-2.47076517208,0.000000000000,4.59039539209,0.517852842654,-7.83392690144,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000,0.000000000000/)
b0 =1.0d0
      !
      re     = 1.336d0
      ae     = 92.2d0/180.0d0*pi
      !
      y1 = (r1 - re)*exp(-b0*(r1-re)**2)
      y2 = (r2 - re)*exp(-b0*(r2-re)**2)
      y3 = cos(alpha) - cos(ae)
      !
 v0 = q(5)*y1**0*y2**0*y3**0
 v1 = q(6)*y1**0*y2**0*y3**1& 
    + q(7)*y1**1*y2**0*y3**0& 
    + q(7)*y1**0*y2**1*y3**0
 v2 = q(8)*y1**0*y2**0*y3**2& 
    + q(9)*y1**1*y2**0*y3**1& 
    + q(9)*y1**0*y2**1*y3**1& 
    + q(10)*y1**1*y2**1*y3**0& 
    + q(11)*y1**2*y2**0*y3**0& 
    + q(11)*y1**0*y2**2*y3**0
 v3 = q(12)*y1**0*y2**0*y3**3& 
    + q(13)*y1**1*y2**0*y3**2& 
    + q(13)*y1**0*y2**1*y3**2& 
    + q(14)*y1**1*y2**1*y3**1& 
    + q(15)*y1**2*y2**0*y3**1& 
    + q(15)*y1**0*y2**2*y3**1& 
    + q(16)*y1**2*y2**1*y3**0& 
    + q(16)*y1**1*y2**2*y3**0& 
    + q(17)*y1**3*y2**0*y3**0& 
    + q(17)*y1**0*y2**3*y3**0

  v4 = q(18)*y1**0*y2**0*y3**4& 
    + q(19)*y1**1*y2**0*y3**3& 
    + q(19)*y1**0*y2**1*y3**3& 
    + q(20)*y1**1*y2**1*y3**2& 
    + q(21)*y1**2*y2**0*y3**2& 
    + q(21)*y1**0*y2**2*y3**2& 
    + q(22)*y1**2*y2**1*y3**1& 
    + q(22)*y1**1*y2**2*y3**1& 
    + q(23)*y1**2*y2**2*y3**0& 
    + q(24)*y1**3*y2**0*y3**1& 
    + q(24)*y1**0*y2**3*y3**1& 
    + q(25)*y1**3*y2**1*y3**0& 
    + q(25)*y1**1*y2**3*y3**0& 
    + q(26)*y1**4*y2**0*y3**0& 
    + q(26)*y1**0*y2**4*y3**0

  v5 = q(27)*y1**0*y2**0*y3**5& 
    + q(28)*y1**1*y2**0*y3**4& 
    + q(28)*y1**0*y2**1*y3**4& 
    + q(29)*y1**1*y2**1*y3**3& 
    + q(30)*y1**2*y2**0*y3**3& 
    + q(30)*y1**0*y2**2*y3**3& 
    + q(31)*y1**2*y2**1*y3**2& 
    + q(31)*y1**1*y2**2*y3**2& 
    + q(32)*y1**2*y2**2*y3**1& 
    + q(33)*y1**3*y2**0*y3**2& 
    + q(33)*y1**0*y2**3*y3**2& 
    + q(34)*y1**3*y2**1*y3**1& 
    + q(34)*y1**1*y2**3*y3**1& 
    + q(35)*y1**3*y2**2*y3**0& 
    + q(35)*y1**2*y2**3*y3**0& 
    + q(36)*y1**4*y2**0*y3**1& 
    + q(36)*y1**0*y2**4*y3**1& 
    + q(37)*y1**4*y2**1*y3**0& 
    + q(37)*y1**1*y2**4*y3**0& 
    + q(38)*y1**5*y2**0*y3**0& 
    + q(38)*y1**0*y2**5*y3**0

 v6 = q(39)*y1**0*y2**0*y3**6& 
    + q(40)*y1**1*y2**0*y3**5& 
    + q(40)*y1**0*y2**1*y3**5& 
    + q(41)*y1**1*y2**1*y3**4& 
    + q(42)*y1**2*y2**0*y3**4& 
    + q(42)*y1**0*y2**2*y3**4& 
    + q(43)*y1**2*y2**1*y3**3& 
    + q(43)*y1**1*y2**2*y3**3& 
    + q(44)*y1**2*y2**2*y3**2& 
    + q(45)*y1**3*y2**0*y3**3& 
    + q(45)*y1**0*y2**3*y3**3& 
    + q(46)*y1**3*y2**1*y3**2& 
    + q(46)*y1**1*y2**3*y3**2& 
    + q(47)*y1**3*y2**2*y3**1& 
    + q(47)*y1**2*y2**3*y3**1& 
    + q(48)*y1**3*y2**3*y3**0& 
    + q(49)*y1**4*y2**0*y3**2& 
    + q(49)*y1**0*y2**4*y3**2& 
    + q(50)*y1**4*y2**1*y3**1& 
    + q(50)*y1**1*y2**4*y3**1& 
    + q(51)*y1**4*y2**2*y3**0& 
    + q(51)*y1**2*y2**4*y3**0& 
    + q(52)*y1**5*y2**0*y3**1& 
    + q(52)*y1**0*y2**5*y3**1& 
    + q(53)*y1**5*y2**1*y3**0& 
    + q(53)*y1**1*y2**5*y3**0& 
    + q(54)*y1**6*y2**0*y3**0& 
    + q(54)*y1**0*y2**6*y3**0


 v7 = q(55)*y1**0*y2**0*y3**7& 
    + q(56)*y1**1*y2**0*y3**6& 
    + q(56)*y1**0*y2**1*y3**6& 
    + q(57)*y1**1*y2**1*y3**5& 
    + q(58)*y1**2*y2**0*y3**5& 
    + q(58)*y1**0*y2**2*y3**5& 
    + q(59)*y1**2*y2**1*y3**4& 
    + q(59)*y1**1*y2**2*y3**4& 
    + q(60)*y1**2*y2**2*y3**3& 
    + q(61)*y1**3*y2**0*y3**4& 
    + q(61)*y1**0*y2**3*y3**4& 
    + q(62)*y1**3*y2**1*y3**3& 
    + q(62)*y1**1*y2**3*y3**3& 
    + q(63)*y1**3*y2**2*y3**2& 
    + q(63)*y1**2*y2**3*y3**2& 
    + q(64)*y1**3*y2**3*y3**1& 
    + q(65)*y1**4*y2**0*y3**3& 
    + q(65)*y1**0*y2**4*y3**3& 
    + q(66)*y1**4*y2**1*y3**2& 
    + q(66)*y1**1*y2**4*y3**2& 
    + q(67)*y1**4*y2**2*y3**1& 
    + q(67)*y1**2*y2**4*y3**1& 
    + q(68)*y1**4*y2**3*y3**0& 
    + q(68)*y1**3*y2**4*y3**0& 
    + q(69)*y1**5*y2**0*y3**2& 
    + q(69)*y1**0*y2**5*y3**2& 
    + q(70)*y1**5*y2**1*y3**1& 
    + q(70)*y1**1*y2**5*y3**1& 
    + q(71)*y1**5*y2**2*y3**0& 
    + q(71)*y1**2*y2**5*y3**0& 
    + q(72)*y1**6*y2**0*y3**1& 
    + q(72)*y1**0*y2**6*y3**1& 
    + q(73)*y1**6*y2**1*y3**0& 
    + q(73)*y1**1*y2**6*y3**0& 
    + q(74)*y1**7*y2**0*y3**0& 
    + q(74)*y1**0*y2**7*y3**0

 v8 = q(75)*y1**0*y2**0*y3**8& 
    + q(76)*y1**1*y2**0*y3**7& 
    + q(76)*y1**0*y2**1*y3**7& 
    + q(77)*y1**1*y2**1*y3**6& 
    + q(78)*y1**2*y2**0*y3**6& 
    + q(78)*y1**0*y2**2*y3**6& 
    + q(79)*y1**2*y2**1*y3**5& 
    + q(79)*y1**1*y2**2*y3**5& 
    + q(80)*y1**2*y2**2*y3**4& 
    + q(81)*y1**3*y2**0*y3**5& 
    + q(81)*y1**0*y2**3*y3**5& 
    + q(82)*y1**3*y2**1*y3**4& 
    + q(82)*y1**1*y2**3*y3**4& 
    + q(83)*y1**3*y2**2*y3**3& 
    + q(83)*y1**2*y2**3*y3**3& 
    + q(84)*y1**3*y2**3*y3**2& 
    + q(85)*y1**4*y2**0*y3**4& 
    + q(85)*y1**0*y2**4*y3**4& 
    + q(86)*y1**4*y2**1*y3**3& 
    + q(86)*y1**1*y2**4*y3**3& 
    + q(87)*y1**4*y2**2*y3**2& 
    + q(87)*y1**2*y2**4*y3**2& 
    + q(88)*y1**4*y2**3*y3**1& 
    + q(88)*y1**3*y2**4*y3**1& 
    + q(89)*y1**4*y2**4*y3**0& 
    + q(90)*y1**5*y2**0*y3**3& 
    + q(90)*y1**0*y2**5*y3**3& 
    + q(91)*y1**5*y2**1*y3**2& 
    + q(91)*y1**1*y2**5*y3**2& 
    + q(92)*y1**5*y2**2*y3**1& 
    + q(92)*y1**2*y2**5*y3**1& 
    + q(93)*y1**5*y2**3*y3**0& 
    + q(93)*y1**3*y2**5*y3**0& 
    + q(94)*y1**6*y2**0*y3**2& 
    + q(94)*y1**0*y2**6*y3**2& 
    + q(95)*y1**6*y2**1*y3**1& 
    + q(95)*y1**1*y2**6*y3**1& 
    + q(96)*y1**6*y2**2*y3**0& 
    + q(96)*y1**2*y2**6*y3**0& 
    + q(97)*y1**7*y2**0*y3**1& 
    + q(97)*y1**0*y2**7*y3**1& 
    + q(98)*y1**7*y2**1*y3**0& 
    + q(98)*y1**1*y2**7*y3**0& 
    + q(99)*y1**8*y2**0*y3**0& 
    + q(99)*y1**0*y2**8*y3**0


    f=(v0+v1+v2+v3+v4+v5+v6+v7+v8)*sin(alpha)
 
    !
  end subroutine dipol_xy2_q_qz