module xy3_alpha
use accuracy
implicit none

integer(ik), parameter :: maxnelem_sym=15, maxnterms=1000, maxnelem=3**4

type ext_func
  integer(ik) :: rank = 1
  integer(ik) :: nelem = 3
  integer(ik) :: nelem_sym = 3
  integer(ik) :: nterms(maxnelem_sym) = 0
  real(ark) :: coef(maxnterms,maxnelem_sym) = 0
end type ext_func


contains


subroutine xy3_polariz_symadap(nelem, extf, natoms, xyz, f)

  integer(ik),intent(in) ::  nelem, natoms
  real(ark),intent(in)   ::  xyz(natoms,3)
  type(ext_func), intent(in) :: extf
  real(ark),intent(out)  ::  f(nelem)

  integer(ik) :: iatom, iterm, imu, nsym, ind(100), i
  real(ark) :: xyz0(3), xyz_(natoms,3), r1, r2, r3, alpha1, alpha2, alpha3, v12(3), v23(3), v31(3), e1(3), e2(3), e3(3), n3(3), rho, tmat(4,3), &!
               coords(7), mu_sym(6), smat(10,6), smat_inv(6,6), smat_10(10,6)

  if (nelem/=6) then
    write(out, '(/a,1x,i3,1x,a)') 'xy3_polariz_symadap error: number of elements for the polarizability tensor =', nelem, ', expected 6'
    stop
  endif

  xyz0 = xyz(1,:)
  do iatom=1, natoms
    xyz_(iatom,:) = xyz(iatom,:) - xyz0(:)
  enddo

  r1 = sqrt(sum(xyz_(2,:)**2))
  r2 = sqrt(sum(xyz_(3,:)**2))
  r3 = sqrt(sum(xyz_(4,:)**2))

  e1 = xyz_(2,:)/r1
  e2 = xyz_(3,:)/r2
  e3 = xyz_(4,:)/r3

  alpha1 = acos(sum(e2*e3))
  alpha2 = acos(sum(e1*e3))
  alpha3 = acos(sum(e1*e2))

  v12 = vector_product(e1,e2)
  v23 = vector_product(e2,e3)
  v31 = vector_product(e3,e1)

  n3 = v12 + v23 + v31

  n3 = n3/sqrt(sum(n3(:)**2))

  rho = acos(sum(n3*e1))

  tmat(1,:) = e1
  tmat(2,:) = e2
  tmat(3,:) = e3
  tmat(4,:) = n3

  coords = (/r1,r2,r3,alpha1,alpha2,alpha3,rho/)

  ! compute six (A, E, E, A) symmetry-adapted components from the input coefficients

  !  1   sqrt(3)*(a11 + a22 + a33)/Rational(3.0)                 sym= A1'
  !  2   sqrt(6)*(Rational(2.0)*a12 - a13 - a23)/Rational(6.0)   sym= Ea'
  !  3   sqrt(2)*(a13 - a23)/Rational(2.0)                       sym= Eb'
  !  4   sqrt(6)*(a14 + a24 - Rational(2.0)*a34)/Rational(6.0)   sym= Ea"
  !  5   sqrt(2)*(a14 - a24)/Rational(2.0)                       sym= Eb"
  !  6   a44                                                     sym= A1'

  mu_sym(1) = fit_xy3_dipole_A_rhosym(extF%nterms(1), extF%coef(1:extF%nterms(1),1), coords)
  mu_sym(2) = fit_xy3_dipole_E_rhosym(1, extF%nterms(2), extF%coef(1:extF%nterms(2),2), coords)
  mu_sym(3) = fit_xy3_dipole_E_rhosym(2, extF%nterms(2), extF%coef(1:extF%nterms(2),2), coords)
  mu_sym(4) = fit_xy3_dipole_E_rhoAsym(1, extF%nterms(3), extF%coef(1:extF%nterms(3),3), coords)
  mu_sym(5) = fit_xy3_dipole_E_rhoAsym(2, extF%nterms(3), extF%coef(1:extF%nterms(3),3), coords)
  mu_sym(6) = fit_xy3_dipole_A_rhosym(extF%nterms(4), extF%coef(1:extF%nterms(4),4), coords)

  nsym = 10
  call symadap_d3h_polariz_10( tmat(1:4,1:3), smat_10(1:nsym,1:6) )

  ! select six required symmetry-adapted combinations

  nsym = 6
  ind(1:nsym) = (/1,5,6,8,9,10/)
  do i=1, nsym
    smat(i,1:6) = smat_10(ind(i),1:6)
  enddo

  ! invert "smat" to obtain symmetry-adapted-to-Cartesian transformation matrix

  call matrix_pseudoinverse_ark( 6, 6, smat(1:6,1:6), smat_inv(1:6,1:6) )

  ! compute six Cartesian components

  f(1:6) = matmul(smat_inv(1:6,1:6), mu_sym(1:6))

  contains

  function vector_product(v1,v2) result (v)
    real(ark), intent(in) :: v1(3), v2(3)
    real(ark) :: v(3)
    v(1) = v1(2)*v2(3)-v1(3)*v2(2)
    v(2) = v1(3)*v2(1)-v1(1)*v2(3)
    v(3) = v1(1)*v2(2)-v1(2)*v2(1)
  end function vector_product

end subroutine xy3_polariz_symadap



function fit_xy3_dipole_A_rhosym(nparams, params, coords) result(f)

  integer(ik), intent(in) :: nparams
  real(ark), intent(in) :: params(nparams), coords(7)
  real(ark) :: f

  real(ark) :: r1, r2, r3, alpha1, alpha2, alpha3, rad, rho, s4, s5

  r1     = coords(1)
  r2     = coords(2)
  r3     = coords(3)
  alpha1 = coords(4)
  alpha2 = coords(5)
  alpha3 = coords(6)
  rho    = coords(7)

  s4 = (2.0_ark*alpha1-alpha2-alpha3)/sqrt(6.0_ark)
  s5 = (alpha2-alpha3)/sqrt(2.0_ark)

  f = xy3_dms_A(nparams, params, (/r1,r2,r3,s4,s5,rho/))

end function fit_xy3_dipole_A_rhosym



function fit_xy3_dipole_E_rhosym(icomp, nparams, params, coords) result(f)

  integer(ik), intent(in) :: icomp, nparams
  real(ark), intent(in) :: params(nparams), coords(7)
  real(ark) :: f

  real(ark) :: r1, r2, r3, alpha1, alpha2, alpha3, rad, rho, s4, s5

  r2     = coords(1)
  r3     = coords(2)
  r1     = coords(3)
  alpha2 = coords(4)
  alpha3 = coords(5)
  alpha1 = coords(6)
  rho    = coords(7)

  s4 = (2.0_ark*alpha1-alpha2-alpha3)/sqrt(6.0_ark)
  s5 = (alpha2-alpha3)/sqrt(2.0_ark)

  f = xy3_dms_E(icomp, nparams, params, (/r1,r2,r3,s4,s5,rho/))

  if (icomp==1) f = -f

end function fit_xy3_dipole_E_rhosym



function fit_xy3_dipole_E_rhoasym(icomp, nparams, params, coords) result(f)

  integer(ik), intent(in) :: icomp, nparams
  real(ark), intent(in) :: params(nparams), coords(7)
  real(ark) :: f

  real(ark) :: r1, r2, r3, alpha1, alpha2, alpha3, rad, rho, s4, s5

  r2     = coords(1)
  r3     = coords(2)
  r1     = coords(3)
  alpha2 = coords(4)
  alpha3 = coords(5)
  alpha1 = coords(6)
  rho    = coords(7)

  s4 = (2.0_ark*alpha1-alpha2-alpha3)/sqrt(6.0_ark)
  s5 = (alpha2-alpha3)/sqrt(2.0_ark)

  f = xy3_dms_E(icomp, nparams, params, (/r1,r2,r3,s4,s5,rho/))

  if (icomp==1) f = -f

  f = f * cos(rho)

end function fit_xy3_dipole_E_rhoasym



subroutine symadap_d3h_polariz_10(tmat, smat)

  real(ark), intent(in) :: tmat(:,:)
  real(ark), intent(out) :: smat(:,:)

  if (any(shape(tmat)/=(/4,3/))) then
    write(out, '(/a,1x,2(1x,i4),1x,a)') 'symadap_d3h_polariz_10 error: dimensions of t-matrix =', shape(tmat), '/= (/4,3/)'
    stop 'STOP'
  endif

  if (any(shape(smat)/=(/10,6/))) then
    write(out, '(/a,1x,2(1x,i4),1x,a)') 'symadap_d3h_polariz_10 error: dimensions of s-matrix =', shape(smat), '/= (/10,6/)'
    stop 'STOP'
  endif

smat( 1 , 1 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)**2 + tmat(2,1)**2 + tmat(3,1)**2)
smat( 1 , 2 )= (2.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(1,2) + tmat(2,1)*tmat(2,2) + &!
      tmat(3,1)*tmat(3,2))
smat( 1 , 3 )= (2.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(1,3) + tmat(2,1)*tmat(2,3) + &!
      tmat(3,1)*tmat(3,3))
smat( 1 , 4 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,2)**2 + tmat(2,2)**2 + tmat(3,2)**2)
smat( 1 , 5 )= (2.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,2)*tmat(1,3) + tmat(2,2)*tmat(2,3) + &!
      tmat(3,2)*tmat(3,3))
smat( 1 , 6 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,3)**2 + tmat(2,3)**2 + tmat(3,3)**2)
smat( 2 , 1 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,1)**2 + tmat(2,1)**2 - 2*tmat(3,1)**2)
smat( 2 , 2 )= (1.0_ark/3.0_ark)*sqrt(6.0_ark)*(tmat(1,1)*tmat(1,2) + tmat(2,1)*tmat(2,2) - 2 &!
      *tmat(3,1)*tmat(3,2))
smat( 2 , 3 )= (1.0_ark/3.0_ark)*sqrt(6.0_ark)*(tmat(1,1)*tmat(1,3) + tmat(2,1)*tmat(2,3) - 2 &!
      *tmat(3,1)*tmat(3,3))
smat( 2 , 4 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,2)**2 + tmat(2,2)**2 - 2*tmat(3,2)**2)
smat( 2 , 5 )= (1.0_ark/3.0_ark)*sqrt(6.0_ark)*(tmat(1,2)*tmat(1,3) + tmat(2,2)*tmat(2,3) - 2 &!
      *tmat(3,2)*tmat(3,3))
smat( 2 , 6 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,3)**2 + tmat(2,3)**2 - 2*tmat(3,3)**2)
smat( 3 , 1 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,1)**2 - tmat(2,1)**2)
smat( 3 , 2 )= sqrt(2.0_ark)*(tmat(1,1)*tmat(1,2) - tmat(2,1)*tmat(2,2))
smat( 3 , 3 )= sqrt(2.0_ark)*(tmat(1,1)*tmat(1,3) - tmat(2,1)*tmat(2,3))
smat( 3 , 4 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,2)**2 - tmat(2,2)**2)
smat( 3 , 5 )= sqrt(2.0_ark)*(tmat(1,2)*tmat(1,3) - tmat(2,2)*tmat(2,3))
smat( 3 , 6 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,3)**2 - tmat(2,3)**2)
smat( 4 , 1 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(2,1) + tmat(1,1)*tmat(3,1) + &!
      tmat(2,1)*tmat(3,1))
smat( 4 , 2 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(2,2) + tmat(1,1)*tmat(3,2) + &!
      tmat(1,2)*tmat(2,1) + tmat(1,2)*tmat(3,1) + tmat(2,1)*tmat(3,2) + &!
      tmat(2,2)*tmat(3,1))
smat( 4 , 3 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(2,3) + tmat(1,1)*tmat(3,3) + &!
      tmat(1,3)*tmat(2,1) + tmat(1,3)*tmat(3,1) + tmat(2,1)*tmat(3,3) + &!
      tmat(2,3)*tmat(3,1))
smat( 4 , 4 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,2)*tmat(2,2) + tmat(1,2)*tmat(3,2) + &!
      tmat(2,2)*tmat(3,2))
smat( 4 , 5 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,2)*tmat(2,3) + tmat(1,2)*tmat(3,3) + &!
      tmat(1,3)*tmat(2,2) + tmat(1,3)*tmat(3,2) + tmat(2,2)*tmat(3,3) + &!
      tmat(2,3)*tmat(3,2))
smat( 4 , 6 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,3)*tmat(2,3) + tmat(1,3)*tmat(3,3) + &!
      tmat(2,3)*tmat(3,3))
smat( 5 , 1 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,1)*tmat(2,1) - tmat(1,1)*tmat(3,1) - &!
      tmat(2,1)*tmat(3,1))
smat( 5 , 2 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,1)*tmat(2,2) - tmat(1,1)*tmat(3,2) + &!
      2*tmat(1,2)*tmat(2,1) - tmat(1,2)*tmat(3,1) - tmat(2,1)*tmat(3,2 &!
      ) - tmat(2,2)*tmat(3,1))
smat( 5 , 3 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,1)*tmat(2,3) - tmat(1,1)*tmat(3,3) + &!
      2*tmat(1,3)*tmat(2,1) - tmat(1,3)*tmat(3,1) - tmat(2,1)*tmat(3,3 &!
      ) - tmat(2,3)*tmat(3,1))
smat( 5 , 4 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,2)*tmat(2,2) - tmat(1,2)*tmat(3,2) - &!
      tmat(2,2)*tmat(3,2))
smat( 5 , 5 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,2)*tmat(2,3) - tmat(1,2)*tmat(3,3) + &!
      2*tmat(1,3)*tmat(2,2) - tmat(1,3)*tmat(3,2) - tmat(2,2)*tmat(3,3 &!
      ) - tmat(2,3)*tmat(3,2))
smat( 5 , 6 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(2*tmat(1,3)*tmat(2,3) - tmat(1,3)*tmat(3,3) - &!
      tmat(2,3)*tmat(3,3))
smat( 6 , 1 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(3,1)*(tmat(1,1) - tmat(2,1))
smat( 6 , 2 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,1)*tmat(3,2) + tmat(1,2)*tmat(3,1) - &!
      tmat(2,1)*tmat(3,2) - tmat(2,2)*tmat(3,1))
smat( 6 , 3 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,1)*tmat(3,3) + tmat(1,3)*tmat(3,1) - &!
      tmat(2,1)*tmat(3,3) - tmat(2,3)*tmat(3,1))
smat( 6 , 4 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(3,2)*(tmat(1,2) - tmat(2,2))
smat( 6 , 5 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,2)*tmat(3,3) + tmat(1,3)*tmat(3,2) - &!
      tmat(2,2)*tmat(3,3) - tmat(2,3)*tmat(3,2))
smat( 6 , 6 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(3,3)*(tmat(1,3) - tmat(2,3))
smat( 7 , 1 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*tmat(4,1)*(tmat(1,1) + tmat(2,1) + tmat(3,1))
smat( 7 , 2 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(4,2) + tmat(1,2)*tmat(4,1) + &!
      tmat(2,1)*tmat(4,2) + tmat(2,2)*tmat(4,1) + tmat(3,1)*tmat(4,2) + &!
      tmat(3,2)*tmat(4,1))
smat( 7 , 3 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,1)*tmat(4,3) + tmat(1,3)*tmat(4,1) + &!
      tmat(2,1)*tmat(4,3) + tmat(2,3)*tmat(4,1) + tmat(3,1)*tmat(4,3) + &!
      tmat(3,3)*tmat(4,1))
smat( 7 , 4 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*tmat(4,2)*(tmat(1,2) + tmat(2,2) + tmat(3,2))
smat( 7 , 5 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*(tmat(1,2)*tmat(4,3) + tmat(1,3)*tmat(4,2) + &!
      tmat(2,2)*tmat(4,3) + tmat(2,3)*tmat(4,2) + tmat(3,2)*tmat(4,3) + &!
      tmat(3,3)*tmat(4,2))
smat( 7 , 6 )= (1.0_ark/3.0_ark)*sqrt(3.0_ark)*tmat(4,3)*(tmat(1,3) + tmat(2,3) + tmat(3,3))
smat( 8 , 1 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*tmat(4,1)*(tmat(1,1) + tmat(2,1) - 2*tmat(3,1 &!
      ))
smat( 8 , 2 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,1)*tmat(4,2) + tmat(1,2)*tmat(4,1) + &!
      tmat(2,1)*tmat(4,2) + tmat(2,2)*tmat(4,1) - 2*tmat(3,1)*tmat(4,2 &!
      ) - 2*tmat(3,2)*tmat(4,1))
smat( 8 , 3 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,1)*tmat(4,3) + tmat(1,3)*tmat(4,1) + &!
      tmat(2,1)*tmat(4,3) + tmat(2,3)*tmat(4,1) - 2*tmat(3,1)*tmat(4,3 &!
      ) - 2*tmat(3,3)*tmat(4,1))
smat( 8 , 4 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*tmat(4,2)*(tmat(1,2) + tmat(2,2) - 2*tmat(3,2 &!
      ))
smat( 8 , 5 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*(tmat(1,2)*tmat(4,3) + tmat(1,3)*tmat(4,2) + &!
      tmat(2,2)*tmat(4,3) + tmat(2,3)*tmat(4,2) - 2*tmat(3,2)*tmat(4,3 &!
      ) - 2*tmat(3,3)*tmat(4,2))
smat( 8 , 6 )= (1.0_ark/6.0_ark)*sqrt(6.0_ark)*tmat(4,3)*(tmat(1,3) + tmat(2,3) - 2*tmat(3,3 &!
      ))
smat( 9 , 1 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(4,1)*(tmat(1,1) - tmat(2,1))
smat( 9 , 2 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,1)*tmat(4,2) + tmat(1,2)*tmat(4,1) - &!
      tmat(2,1)*tmat(4,2) - tmat(2,2)*tmat(4,1))
smat( 9 , 3 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,1)*tmat(4,3) + tmat(1,3)*tmat(4,1) - &!
      tmat(2,1)*tmat(4,3) - tmat(2,3)*tmat(4,1))
smat( 9 , 4 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(4,2)*(tmat(1,2) - tmat(2,2))
smat( 9 , 5 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*(tmat(1,2)*tmat(4,3) + tmat(1,3)*tmat(4,2) - &!
      tmat(2,2)*tmat(4,3) - tmat(2,3)*tmat(4,2))
smat( 9 , 6 )= (1.0_ark/2.0_ark)*sqrt(2.0_ark)*tmat(4,3)*(tmat(1,3) - tmat(2,3))
smat( 10 , 1 )= tmat(4,1)**2
smat( 10 , 2 )= 2*tmat(4,1)*tmat(4,2)
smat( 10 , 3 )= 2*tmat(4,1)*tmat(4,3)
smat( 10 , 4 )= tmat(4,2)**2
smat( 10 , 5 )= 2*tmat(4,2)*tmat(4,3)
smat( 10 , 6 )= tmat(4,3)**2

end subroutine symadap_d3h_polariz_10



subroutine matrix_pseudoinverse_ark(m, n, mat, invmat, info)

  integer(ik), intent(in) :: m, n
  real(ark), intent(in) :: mat(m,n)
  real(ark), intent(out) :: invmat(n,m)
  integer(ik), intent(out), optional :: info

  integer(ik) lwork, info_, i, j
  double precision work1(1), matd(m,n), matu(m,m), matvt(n,n), invmatd(n,m), mat_d(n,m), sv(n), tmat(n,m), tol
  double precision, allocatable :: work(:)

  tol = 1.0d-08 !epsilon(1.0d0)

  matd = dble(mat)

  lwork = -1
  call dgesvd('A', 'A', m, n, matd, m, sv, matu, m, matvt, n, work1, lwork, info_)
  lwork = int(work1(1), kind=ik)

  allocate(work(lwork), stat=info_)
  if (info_/=0) then
    write(out, '(/a,1x,i6)') 'matrix_pseudoinverse_ark error: failed to allocate workspace array required for SVD, size =', lwork
    stop
  endif

  call dgesvd('A', 'A', m, n, matd, m, sv, matu, m, matvt, n, work, lwork, info_)

  if (info_/=0) then
    write(out, '(/a,1x,i6)') 'matrix_pseudoinverse_ark error: SVD failed, info =', info_
    stop
  endif

  if (present(info)) info = 0

  mat_d = 0.0
  do i=1, n
    if (sv(i)>=tol) then
      mat_d(i,i) = 1.0d0/sv(i)
    else
      if (present(info)) then
        info = i
        mat_d(i,i) = 0.0
      else
        write(out, '(/a)') 'matrix_pseudoinverse_ark error: matrix is singular'
        write(out, '(a)') 'matrix:'
        do j=1, m
          write(out, '(100(1x,f10.6))') mat(j,1:n)
        enddo
        write(out, '(a)') 'singular elements:'
        do j=1, n
          write(out, '(1x,f12.8)') sv(j)
        enddo
        stop 'STOP'
      endif
    endif
  enddo

  call dgemm('N', 'T', n, m, m, 1.0d0, mat_d(1:n,1:m), n, matu(1:m,1:m), m, 0.0d0, tmat(1:n,1:m), n)
  call dgemm('T', 'N', n, m, n, 1.0d0, matvt(1:n,1:n), n, tmat(1:n,1:m), n, 0.0d0, invmatd(1:n,1:m), n)

  invmat = real(invmatd,kind=ark)

  deallocate(work)

end subroutine matrix_pseudoinverse_ark



function xy3_dms_A(parmax, param, local) result(f)

    integer, intent(in)   ::  parmax
    real(ark), intent(in) ::  param(parmax)
    real(ark), intent(in) ::  local(6)
    real(ark)             ::  f

    real(ark) r14, r24, r34, alpha1, alpha2, alpha3, rhobar
    real(ark) y1,y2,y3,y4,y5,alpha,sinrho,cosrho,drho
    real(ark) v0,v1,v2,v3,v4,v5,v6
    real(ark) rhoe, re14, beta, de

    real(ark)                                        &!
      fea    ,fea1  ,                                    &!
      fea11  ,fea12  ,fea14  ,fea44  ,                   &!
      fea111 ,fea112 ,fea114 ,fea123 ,                   &!
      fea124 ,fea144 ,fea155 ,fea455 ,                   &!
      fea1111,fea1112,fea1114,fea1122,                   &!
      fea1123,fea1124,fea1125,fea1144,                   &!
      fea1155,fea1244,fea1255,fea1444,                   &!
      fea1455,fea4444

    real(ark) ::   &!
      fea44444 ,fea33455 ,fea33445 ,fea33345 ,fea33344 ,&!
      fea33334 ,fea33333 ,fea25555 ,fea24455 ,fea24445 ,fea23333 ,&!
      fea13455 ,fea13445 ,fea13345 ,fea12355 ,fea11334 ,fea11333 ,&!
      fea11255 ,fea11245 ,fea11234 ,fea11233 ,fea11135 ,fea11134 ,&!
      fea11123 ,fea555555,fea444444,fea335555,fea334455,fea334445,&!
      fea333555,fea333333,fea244555,fea244455,fea233445,fea233444,&!
      fea233345,fea233344,fea233335,fea223355,fea222335,fea222334,&!
      fea222333,fea222255,fea222245,fea222233,fea222224,fea145555,&!
      fea134444,fea133444,fea133345,fea133334,fea133333,fea124555,&!
      fea124455,fea123455,fea123345,fea113555,fea113345,fea112355,&!
      fea112335,fea112233,fea111444,fea111234,fea111233,fea111123

    real(ark) :: s1,s2,s3,s4,s5,s6

    real(ark) :: &!
      f1a,f2a,f3a,f4a,f5a,f6a,f7a,f0a1, &!
      f1a1,f2a1,f3a1,f4a1,f5a1,f6a1,  &!
      f0a11,f1a11,f2a11,f3a11,f4a11, &!
      f0a12,f1a12,f2a12,f3a12,f4a12, &!
      f0a14,f1a14,f2a14,f3a14,f4a14, &!
      f0a44,f1a44,f2a44,f3a44,f4a44, &!
      f0a111,f1a111,f2a111,f3a111  , &!
      f0a112,f1a112,f2a112,f3a112  , &!
      f0a114,f1a114,f2a114,f3a114  , &!
      f0a123,f1a123,f2a123,f3a123  , &!
      f0a124,f1a124,f2a124,f3a124  , &!
      f0a144,f1a144,f2a144,f3a144  , &!
      f0a155,f1a155,f2a155,f3a155  , &!
      f0a455,f1a455,f2a455,f3a455  , &!
      f0a1111,f1a1111,f2a1111      , &!
      f0a1112,f1a1112,f2a1112      , &!
      f0a1114,f1a1114,f2a1114      , &!
      f0a1122,f1a1122,f2a1122      , &!
      f0a1123,f1a1123,f2a1123      , &!
      f0a1124,f1a1124,f2a1124      , &!
      f0a1125,f1a1125,f2a1125      , &!
      f0a1144,f1a1144,f2a1144      , &!
      f0a1155,f1a1155,f2a1155      , &!
      f0a1244,f1a1244,f2a1244      , &!
      f0a1255,f1a1255,f2a1255      , &!
      f0a1444,f1a1444,f2a1444      , &!
      f0a1455,f1a1455,f2a1455      , &!
      f0a4444,f1a4444,f2a4444      , &!
      f0a44444 ,f1a44444 ,           &!
      f2a44444 ,f0a33455 ,f1a33455 ,f2a33455 ,f0a33445 ,f1a33445 ,&!
      f2a33445 ,f0a33345 ,f1a33345 ,f2a33345 ,f0a33344 ,f1a33344 ,&!
      f2a33344 ,f0a33334 ,f1a33334 ,f2a33334 ,f0a33333 ,f1a33333 ,&!
      f2a33333 ,f0a25555 ,f1a25555 ,f2a25555 ,f0a24455 ,f1a24455 ,&!
      f2a24455 ,f0a24445 ,f1a24445 ,f2a24445 ,f0a23333 ,f1a23333 ,&!
      f2a23333 ,f0a13455 ,f1a13455 ,f2a13455 ,f0a13445 ,f1a13445 ,&!
      f2a13445 ,f0a13345 ,f1a13345 ,f2a13345 ,f0a12355 ,f1a12355 ,&!
      f2a12355 ,f0a11334 ,f1a11334 ,f2a11334 ,f0a11333 ,f1a11333 ,&!
      f2a11333 ,f0a11255 ,f1a11255 ,f2a11255 ,f0a11245 ,f1a11245 ,&!
      f2a11245 ,f0a11234 ,f1a11234 ,f2a11234 ,f0a11233 ,f1a11233 ,&!
      f2a11233 ,f0a11135 ,f1a11135 ,f2a11135 ,f0a11134 ,f1a11134 ,&!
      f2a11134 ,f0a11123 ,f1a11123 ,f2a11123 ,f0a555555,f1a555555,&!
      f2a555555,f0a444444,f1a444444,f2a444444,f0a335555,f1a335555,&!
      f2a335555,f0a334455,f1a334455,f2a334455,f0a334445,f1a334445,&!
      f2a334445,f0a333555,f1a333555,f2a333555,f0a333333,f1a333333,&!
      f2a333333,f0a244555,f1a244555,f2a244555,f0a244455,f1a244455,&!
      f2a244455,f0a233445,f1a233445,f2a233445,f0a233444,f1a233444,&!
      f2a233444,f0a233345,f1a233345,f2a233345,f0a233344,f1a233344,&!
      f2a233344,f0a233335,f1a233335,f2a233335,f0a223355,f1a223355,&!
      f2a223355,f0a222335,f1a222335,f2a222335,f0a222334,f1a222334,&!
      f2a222334,f0a222333,f1a222333,f2a222333,f0a222255,f1a222255,&!
      f2a222255,f0a222245,f1a222245,f2a222245,f0a222233,f1a222233,&!
      f2a222233,f0a222224,f1a222224,f2a222224,f0a145555,f1a145555,&!
      f2a145555,f0a134444,f1a134444,f2a134444,f0a133444,f1a133444,&!
      f2a133444,f0a133345,f1a133345,f2a133345,f0a133334,f1a133334,&!
      f2a133334,f0a133333,f1a133333,f2a133333,f0a124555,f1a124555,&!
      f2a124555,f0a124455,f1a124455,f2a124455,f0a123455,f1a123455,&!
      f2a123455,f0a123345,f1a123345,f2a123345,f0a113555,f1a113555,&!
      f2a113555,f0a113345,f1a113345,f2a113345,f0a112355,f1a112355,&!
      f2a112355,f0a112335,f1a112335,f2a112335,f0a112233,f1a112233,&!
      f2a112233,f0a111444,f1a111444,f2a111444,f0a111234,f1a111234,&!
      f2a111234,f0a111233,f1a111233,f2a111233,f0a111123,f1a111123,&!
      f2a111123



   rhoe       = param(  1)* pi / 180.0_ark
   re14       = param( 2)
   beta       = param( 3)**2
   de         = param( 4)
   f1a        = param(  5)
   f2a        = param(  6)
   f3a        = param(  7)
   f4a        = param(  8)
   f5a        = param(  9)
   f6a        = param( 10)
   f7a        = param( 11)
   f0a1       = param( 12)
   f1a1       = param( 13)
   f2a1       = param( 14)
   f3a1       = param( 15)
   f4a1       = param( 16)
   f5a1       = param( 17)
   f6a1       = param( 18)
   f0a11      = param( 19)
   f1a11      = param( 20)
   f2a11      = param( 21)
   f3a11      = param( 22)
   f4a11      = param( 23)
   f0a12      = param( 24)
   f1a12      = param( 25)
   f2a12      = param( 26)
   f3a12      = param( 27)
   f4a12      = param( 28)
   f0a14      = param( 29)
   f1a14      = param( 30)
   f2a14      = param( 31)
   f3a14      = param( 32)
   f4a14      = param( 33)
   f0a44      = param( 34)
   f1a44      = param( 35)
   f2a44      = param( 36)
   f3a44      = param( 37)
   f4a44      = param( 38)
   f0a111     = param( 39)
   f1a111     = param( 40)
   f2a111     = param( 41)
   f3a111     = param( 42)
   f0a112     = param( 43)
   f1a112     = param( 44)
   f2a112     = param( 45)
   f3a112     = param( 46)
   f0a114     = param( 47)
   f1a114     = param( 48)
   f2a114     = param( 49)
   f3a114     = param( 50)
   f0a123     = param( 51)
   f1a123     = param( 52)
   f2a123     = param( 53)
   f3a123     = param( 54)
   f0a124     = param( 55)
   f1a124     = param( 56)
   f2a124     = param( 57)
   f3a124     = param( 58)
   f0a144     = param( 59)
   f1a144     = param( 60)
   f2a144     = param( 61)
   f3a144     = param( 62)
   f0a155     = param( 63)
   f1a155     = param( 64)
   f2a155     = param( 65)
   f3a155     = param( 66)
   f0a455     = param( 67)
   f1a455     = param( 68)
   f2a455     = param( 69)
   f3a455     = param( 70)
   f0a1111    = param( 71)
   f1a1111    = param( 72)
   f2a1111    = param( 73)
   f0a1112    = param( 74)
   f1a1112    = param( 75)
   f2a1112    = param( 76)
   f0a1114    = param( 77)
   f1a1114    = param( 78)
   f2a1114    = param( 79)
   f0a1122    = param( 80)
   f1a1122    = param( 81)
   f2a1122    = param( 82)
   f0a1123    = param( 83)
   f1a1123    = param( 84)
   f2a1123    = param( 85)
   f0a1124    = param( 86)
   f1a1124    = param( 87)
   f2a1124    = param( 88)
   f0a1125    = param( 89)
   f1a1125    = param( 90)
   f2a1125    = param( 91)
   f0a1144    = param( 92)
   f1a1144    = param( 93)
   f2a1144    = param( 94)
   f0a1155    = param( 95)
   f1a1155    = param( 96)
   f2a1155    = param( 97)
   f0a1244    = param( 98)
   f1a1244    = param( 99)
   f2a1244    = param(100)
   f0a1255    = param(101)
   f1a1255    = param(102)
   f2a1255    = param(103)
   f0a1444    = param(104)
   f1a1444    = param(105)
   f2a1444    = param(106)
   f0a1455    = param(107)
   f1a1455    = param(108)
   f2a1455    = param(109)
   f0a4444    = param(110)
   f1a4444    = param(111)
   f2a4444    = param(112)
   !
   !
   if (parmax>112) then
     !
     f0a44444   = param(113)
     f1a44444   = param(114)
     f2a44444   = param(115)
     f0a33455   = param(116)
     f1a33455   = param(117)
     f2a33455   = param(118)
     f0a33445   = param(119)
     f1a33445   = param(120)
     f2a33445   = param(121)
     f0a33345   = param(122)
     f1a33345   = param(123)
     f2a33345   = param(124)
     f0a33344   = param(125)
     f1a33344   = param(126)
     f2a33344   = param(127)
     f0a33334   = param(128)
     f1a33334   = param(129)
     f2a33334   = param(130)
     f0a33333   = param(131)
     f1a33333   = param(132)
     f2a33333   = param(133)
     f0a25555   = param(134)
     f1a25555   = param(135)
     f2a25555   = param(136)
     f0a24455   = param(137)
     f1a24455   = param(138)
     f2a24455   = param(139)
     f0a24445   = param(140)
     f1a24445   = param(141)
     f2a24445   = param(142)
     f0a23333   = param(143)
     f1a23333   = param(144)
     f2a23333   = param(145)
     f0a13455   = param(146)
     f1a13455   = param(147)
     f2a13455   = param(148)
     f0a13445   = param(149)
     f1a13445   = param(150)
     f2a13445   = param(151)
     f0a13345   = param(152)
     f1a13345   = param(153)
     f2a13345   = param(154)
     f0a12355   = param(155)
     f1a12355   = param(156)
     f2a12355   = param(157)
     f0a11334   = param(158)
     f1a11334   = param(159)
     f2a11334   = param(160)
     f0a11333   = param(161)
     f1a11333   = param(162)
     f2a11333   = param(163)
     f0a11255   = param(164)
     f1a11255   = param(165)
     f2a11255   = param(166)
     f0a11245   = param(167)
     f1a11245   = param(168)
     f2a11245   = param(169)
     f0a11234   = param(170)
     f1a11234   = param(171)
     f2a11234   = param(172)
     f0a11233   = param(173)
     f1a11233   = param(174)
     f2a11233   = param(175)
     f0a11135   = param(176)
     f1a11135   = param(177)
     f2a11135   = param(178)
     f0a11134   = param(179)
     f1a11134   = param(180)
     f2a11134   = param(181)
     f0a11123   = param(182)
     f1a11123   = param(183)
     f2a11123   = param(184)
     f0a555555  = param(185)
     f1a555555  = param(186)
     f2a555555  = param(187)
     f0a444444  = param(188)
     f1a444444  = param(189)
     f2a444444  = param(190)
     f0a335555  = param(191)
     f1a335555  = param(192)
     f2a335555  = param(193)
     f0a334455  = param(194)
     f1a334455  = param(195)
     f2a334455  = param(196)
     f0a334445  = param(197)
     f1a334445  = param(198)
     f2a334445  = param(199)
     f0a333555  = param(200)
     f1a333555  = param(201)
     f2a333555  = param(202)
     f0a333333  = param(203)
     f1a333333  = param(204)
     f2a333333  = param(205)
     f0a244555  = param(206)
     f1a244555  = param(207)
     f2a244555  = param(208)
     f0a244455  = param(209)
     f1a244455  = param(210)
     f2a244455  = param(211)
     f0a233445  = param(212)
     f1a233445  = param(213)
     f2a233445  = param(214)
     f0a233444  = param(215)
     f1a233444  = param(216)
     f2a233444  = param(217)
     f0a233345  = param(218)
     f1a233345  = param(219)
     f2a233345  = param(220)
     f0a233344  = param(221)
     f1a233344  = param(222)
     f2a233344  = param(223)
     f0a233335  = param(224)
     f1a233335  = param(225)
     f2a233335  = param(226)
     f0a223355  = param(227)
     f1a223355  = param(228)
     f2a223355  = param(229)
     f0a222335  = param(230)
     f1a222335  = param(231)
     f2a222335  = param(232)
     f0a222334  = param(233)
     f1a222334  = param(234)
     f2a222334  = param(235)
     f0a222333  = param(236)
     f1a222333  = param(237)
     f2a222333  = param(238)
     f0a222255  = param(239)
     f1a222255  = param(240)
     f2a222255  = param(241)
     f0a222245  = param(242)
     f1a222245  = param(243)
     f2a222245  = param(244)
     f0a222233  = param(245)
     f1a222233  = param(246)
     f2a222233  = param(247)
     f0a222224  = param(248)
     f1a222224  = param(249)
     f2a222224  = param(250)
     f0a145555  = param(251)
     f1a145555  = param(252)
     f2a145555  = param(253)
     f0a134444  = param(254)
     f1a134444  = param(255)
     f2a134444  = param(256)
     f0a133444  = param(257)
     f1a133444  = param(258)
     f2a133444  = param(259)
     f0a133345  = param(260)
     f1a133345  = param(261)
     f2a133345  = param(262)
     f0a133334  = param(263)
     f1a133334  = param(264)
     f2a133334  = param(265)
     f0a133333  = param(266)
     f1a133333  = param(267)
     f2a133333  = param(268)
     f0a124555  = param(269)
     f1a124555  = param(270)
     f2a124555  = param(271)
     f0a124455  = param(272)
     f1a124455  = param(273)
     f2a124455  = param(274)
     f0a123455  = param(275)
     f1a123455  = param(276)
     f2a123455  = param(277)
     f0a123345  = param(278)
     f1a123345  = param(279)
     f2a123345  = param(280)
     f0a113555  = param(281)
     f1a113555  = param(282)
     f2a113555  = param(283)
     f0a113345  = param(284)
     f1a113345  = param(285)
     f2a113345  = param(286)
     f0a112355  = param(287)
     f1a112355  = param(288)
     f2a112355  = param(289)
     f0a112335  = param(290)
     f1a112335  = param(291)
     f2a112335  = param(292)
     f0a112233  = param(293)
     f1a112233  = param(294)
     f2a112233  = param(295)
     f0a111444  = param(296)
     f1a111444  = param(297)
     f2a111444  = param(298)
     f0a111234  = param(299)
     f1a111234  = param(300)
     f2a111234  = param(301)
     f0a111233  = param(302)
     f1a111233  = param(303)
     f2a111233  = param(304)
     f0a111123  = param(305)
     f1a111123  = param(306)
     f2a111123  = param(307)
     !
   endif

  r14    = local(1) ;  r24    = local(2) ;  r34    = local(3)
  y4 = local(4) ;  y5 = local(5) ;  rhobar = local(6)

  sinrho = sin(rhobar)
  cosrho = cos(rhobar)

  drho=(sin(rhoe)-sinrho)

  y1=1.00_ark*(r14-re14) *exp(-beta*(r14-re14)**2)
  y2=1.00_ark*(r24-re14) *exp(-beta*(r24-re14)**2)
  y3=1.00_ark*(r34-re14) *exp(-beta*(r34-re14)**2)

  v0=de+f1a*drho+f2a*drho**2+f3a*drho**3+f4a*drho**4+f5a*drho**5 &!
        +f6a*drho**6+f7a*drho**7

  fea1= f0a1+f1a1*drho+f2a1*drho**2+f3a1*drho**3+f4a1*drho**4+f5a1*drho**5+f6a1*drho**6

  fea11=   f0a11+f1a11*drho+f2a11*drho**2+f3a11*drho**3+f4a11*drho**4
  fea12=   f0a12+f1a12*drho+f2a12*drho**2+f3a12*drho**3+f4a12*drho**4
  fea14=   f0a14+f1a14*drho+f2a14*drho**2+f3a14*drho**3+f4a14*drho**4
  fea44=   f0a44+f1a44*drho+f2a44*drho**2+f3a44*drho**3+f4a44*drho**4

  fea111= f0a111+f1a111*drho+f2a111*drho**2+f3a111*drho**3
  fea112= f0a112+f1a112*drho+f2a112*drho**2+f3a112*drho**3
  fea114= f0a114+f1a114*drho+f2a114*drho**2+f3a114*drho**3
  fea123= f0a123+f1a123*drho+f2a123*drho**2+f3a123*drho**3
  fea124= f0a124+f1a124*drho+f2a124*drho**2+f3a124*drho**3
  fea144= f0a144+f1a144*drho+f2a144*drho**2+f3a144*drho**3
  fea155= f0a155+f1a155*drho+f2a155*drho**2+f3a155*drho**3
  fea455= f0a455+f1a455*drho+f2a455*drho**2+f3a455*drho**3

  fea1111= f0a1111+f1a1111*drho+f2a1111*drho**2
  fea1112= f0a1112+f1a1112*drho+f2a1112*drho**2
  fea1114= f0a1114+f1a1114*drho+f2a1114*drho**2
  fea1122= f0a1122+f1a1122*drho+f2a1122*drho**2
  fea1123= f0a1123+f1a1123*drho+f2a1123*drho**2
  fea1124= f0a1124+f1a1124*drho+f2a1124*drho**2
  fea1125= f0a1125+f1a1125*drho+f2a1125*drho**2
  fea1144= f0a1144+f1a1144*drho+f2a1144*drho**2
  fea1155= f0a1155+f1a1155*drho+f2a1155*drho**2
  fea1244= f0a1244+f1a1244*drho+f2a1244*drho**2
  fea1255= f0a1255+f1a1255*drho+f2a1255*drho**2
  fea1444= f0a1444+f1a1444*drho+f2a1444*drho**2
  fea1455= f0a1455+f1a1455*drho+f2a1455*drho**2
  fea4444= f0a4444+f1a4444*drho+f2a4444*drho**2


  if (parmax>112) then
     !
     fea44444 = f0a44444  + f1a44444 *drho+ f2a44444 *drho**2
     fea33455 = f0a33455  + f1a33455 *drho+ f2a33455 *drho**2
     fea33445 = f0a33445  + f1a33445 *drho+ f2a33445 *drho**2
     fea33345 = f0a33345  + f1a33345 *drho+ f2a33345 *drho**2
     fea33344 = f0a33344  + f1a33344 *drho+ f2a33344 *drho**2
     fea33334 = f0a33334  + f1a33334 *drho+ f2a33334 *drho**2
     fea33333 = f0a33333  + f1a33333 *drho+ f2a33333 *drho**2
     fea25555 = f0a25555  + f1a25555 *drho+ f2a25555 *drho**2
     fea24455 = f0a24455  + f1a24455 *drho+ f2a24455 *drho**2
     fea24445 = f0a24445  + f1a24445 *drho+ f2a24445 *drho**2
     fea23333 = f0a23333  + f1a23333 *drho+ f2a23333 *drho**2
     fea13455 = f0a13455  + f1a13455 *drho+ f2a13455 *drho**2
     fea13445 = f0a13445  + f1a13445 *drho+ f2a13445 *drho**2
     fea13345 = f0a13345  + f1a13345 *drho+ f2a13345 *drho**2
     fea12355 = f0a12355  + f1a12355 *drho+ f2a12355 *drho**2
     fea11334 = f0a11334  + f1a11334 *drho+ f2a11334 *drho**2
     fea11333 = f0a11333  + f1a11333 *drho+ f2a11333 *drho**2
     fea11255 = f0a11255  + f1a11255 *drho+ f2a11255 *drho**2
     fea11245 = f0a11245  + f1a11245 *drho+ f2a11245 *drho**2
     fea11234 = f0a11234  + f1a11234 *drho+ f2a11234 *drho**2
     fea11233 = f0a11233  + f1a11233 *drho+ f2a11233 *drho**2
     fea11135 = f0a11135  + f1a11135 *drho+ f2a11135 *drho**2
     fea11134 = f0a11134  + f1a11134 *drho+ f2a11134 *drho**2
     fea11123 = f0a11123  + f1a11123 *drho+ f2a11123 *drho**2
     fea555555= f0a555555 + f1a555555*drho+ f2a555555*drho**2
     fea444444= f0a444444 + f1a444444*drho+ f2a444444*drho**2
     fea335555= f0a335555 + f1a335555*drho+ f2a335555*drho**2
     fea334455= f0a334455 + f1a334455*drho+ f2a334455*drho**2
     fea334445= f0a334445 + f1a334445*drho+ f2a334445*drho**2
     fea333555= f0a333555 + f1a333555*drho+ f2a333555*drho**2
     fea333333= f0a333333 + f1a333333*drho+ f2a333333*drho**2
     fea244555= f0a244555 + f1a244555*drho+ f2a244555*drho**2
     fea244455= f0a244455 + f1a244455*drho+ f2a244455*drho**2
     fea233445= f0a233445 + f1a233445*drho+ f2a233445*drho**2
     fea233444= f0a233444 + f1a233444*drho+ f2a233444*drho**2
     fea233345= f0a233345 + f1a233345*drho+ f2a233345*drho**2
     fea233344= f0a233344 + f1a233344*drho+ f2a233344*drho**2
     fea233335= f0a233335 + f1a233335*drho+ f2a233335*drho**2
     fea223355= f0a223355 + f1a223355*drho+ f2a223355*drho**2
     fea222335= f0a222335 + f1a222335*drho+ f2a222335*drho**2
     fea222334= f0a222334 + f1a222334*drho+ f2a222334*drho**2
     fea222333= f0a222333 + f1a222333*drho+ f2a222333*drho**2
     fea222255= f0a222255 + f1a222255*drho+ f2a222255*drho**2
     fea222245= f0a222245 + f1a222245*drho+ f2a222245*drho**2
     fea222233= f0a222233 + f1a222233*drho+ f2a222233*drho**2
     fea222224= f0a222224 + f1a222224*drho+ f2a222224*drho**2
     fea145555= f0a145555 + f1a145555*drho+ f2a145555*drho**2
     fea134444= f0a134444 + f1a134444*drho+ f2a134444*drho**2
     fea133444= f0a133444 + f1a133444*drho+ f2a133444*drho**2
     fea133345= f0a133345 + f1a133345*drho+ f2a133345*drho**2
     fea133334= f0a133334 + f1a133334*drho+ f2a133334*drho**2
     fea133333= f0a133333 + f1a133333*drho+ f2a133333*drho**2
     fea124555= f0a124555 + f1a124555*drho+ f2a124555*drho**2
     fea124455= f0a124455 + f1a124455*drho+ f2a124455*drho**2
     fea123455= f0a123455 + f1a123455*drho+ f2a123455*drho**2
     fea123345= f0a123345 + f1a123345*drho+ f2a123345*drho**2
     fea113555= f0a113555 + f1a113555*drho+ f2a113555*drho**2
     fea113345= f0a113345 + f1a113345*drho+ f2a113345*drho**2
     fea112355= f0a112355 + f1a112355*drho+ f2a112355*drho**2
     fea112335= f0a112335 + f1a112335*drho+ f2a112335*drho**2
     fea112233= f0a112233 + f1a112233*drho+ f2a112233*drho**2
     fea111444= f0a111444 + f1a111444*drho+ f2a111444*drho**2
     fea111234= f0a111234 + f1a111234*drho+ f2a111234*drho**2
     fea111233= f0a111233 + f1a111233*drho+ f2a111233*drho**2
     fea111123= f0a111123 + f1a111123*drho+ f2a111123*drho**2
    !
  endif



      v1 = (y3+y2+y1)*fea1

      v2 = (y2*y3+y1*y3+y1*y2)*fea12                                                                 &!
       +(y2**2+y3**2+y1**2)*fea11                                                                    &!
       +(-sqrt(3.0_ark)*y3*y5/2.0_ark-y3*y4/2.0_ark+y1*y4+sqrt(3.0_ark)*y2*y5/2.0_ark-y2*y4/2.0_ark)*fea14 &!
       +(y5**2+y4**2)*fea44

      v3 = (y1*y3*y4+y1*y2*y4-2.0_ark*y2*y3*y4+sqrt(3.0_ark)*y1*y2*y5-sqrt(3.0_ark)*y1*y3*y5)*fea124   &!
       +(3.0_ark/4.0_ark*y3*y4**2-sqrt(3.0_ark)*y3*y4*y5/2.0_ark+y1*y5**2+y2*y5**2/4.0_ark               &!
       +3.0_ark/4.0_ark*y2*y4**2+sqrt(3.0_ark)*y2*y4*y5/2.0_ark+y3*y5**2/4.0_ark)*fea155                 &!
       +(y2*y3**2+y1*y3**2+y1**2*y3+y1*y2**2+y2**2*y3+y1**2*y2)*fea112+                             &!
       (-y4**3/3.0_ark+y4*y5**2)*fea455+fea123*y1*y2*y3                                              &!
       +(y1*y4**2+3.0_ark/4.0_ark*y3*y5**2+3.0_ark/4.0_ark*y2*y5**2+y2*y4**2/4.0_ark                     &!
       -sqrt(3.0_ark)*y2*y4*y5/2.0_ark+sqrt(3.0_ark)*y3*y4*y5/2.0_ark+y3*y4**2/4.0_ark)*fea144           &!
       +(y3**3+y2**3+y1**3)*fea111+(-y2**2*y4/2.0_ark-y3**2*y4/2.0_ark+sqrt(3.0_ark)*y2**2*y5/2.0_ark   &!
       +y1**2*y4-sqrt(3.0_ark)*y3**2*y5/2.0_ark)*fea114
       !
      s2 = (y4**4+y5**4+2.0_ark*y4**2*y5**2)*fea4444+(3.0_ark/8.0_ark*sqrt(3.0_ark)*&!
       y2*y5**3-3.0_ark/8.0_ark*sqrt(3.0_ark)*y3*y4**2*y5-3.0_ark/8.0_ark*sqrt(3.0_ark)*y3*&!
       y5**3-9.0_ark/8.0_ark*y2*y4*y5**2-y3*y4**3/8.0_ark-y2*y4**3/8.0_ark-9.0_ark/8.0_ark*&!
       y3*y4*y5**2+y1*y4**3+3.0_ark/8.0_ark*sqrt(3.0_ark)*y2*y4**2*y5)*fea1444 &!
       +(3.0_ark/4.0_ark*y2**2*y4**2+3.0_ark/4.0_ark*y3**2*y4**2+y1**2*y5**2+y3**2*y5**2/4.0_ark &!
       -sqrt(3.0_ark)*y3**2*y4*y5/2.0_ark+sqrt(3.0_ark)*y2**2*y4*y5/2.0_ark+y2**2&!
       *y5**2/4.0_ark)*fea1155
       s1 = s2+(y3**2*y4**2/4.0_ark+3.0_ark/4.0_ark*y3**2*y5**2+y1**2*y4**2+y2**2*&!
       y4**2/4.0_ark+sqrt(3.0_ark)*y3**2*y4*y5/2.0_ark-sqrt(3.0_ark)*y2**2*y4*y5/2.0_ark&!
       +3.0_ark/4.0_ark*y2**2*y5**2)*fea1144+(y1**3*y4+sqrt(3.0_ark)*y2**3*y5/2.0_ark&!
       -sqrt(3.0_ark)*y3**3*y5/2.0_ark-y2**3*y4/2.0_ark-y3**3*y4/2.0_ark)*fea1114&!
       +(y2**4+y1**4+y3**4)*fea1111+(sqrt(3.0_ark)*y1*y3*y4*y5+3.0_ark/2.0_ark*y2*y3*y5**2&!
       -y2*y3*y4**2/2.0_ark+y1*y2*y4**2-sqrt(3.0_ark)*y1*y2*y4*y5+y1*y3*y4**2)*fea1244
       !
      s2 = s1+(y1*y3*y5**2+y1*y2*y5**2-sqrt(3.0_ark)*y1*y3*y4*y5-y2*y3*y5**&!
       2/2.0_ark+3.0_ark/2.0_ark*y2*y3*y4**2+sqrt(3.0_ark)*y1*y2*y4*y5)*fea1255+&!
       (-y1*y3**2*y4/2.0_ark+y1**2*y3*y4-sqrt(3.0_ark)*y1*y3**2*y5/2.0_ark-sqrt(3.0_ark)*y2&!
       *y3**2*y5/2.0_ark+y1**2*y2*y4+sqrt(3.0_ark)*y2**2*y3*y5/2.0_ark-y2**2*y3*y4&!
       /2.0_ark+sqrt(3.0_ark)*y1*y2**2*y5/2.0_ark-y2*y3**2*y4/2.0_ark-y1*y2**2*y4/2.0_ark&!
       )*fea1124+(y1**2*y2*y5+sqrt(3.0_ark)*y1*y3**2*y4/2.0_ark+sqrt(3.0_ark)*y1*&!
       y2**2*y4/2.0_ark-sqrt(3.0_ark)*y2*y3**2*y4/2.0_ark-sqrt(3.0_ark)*y2**2*y3*y4/2.0_ark&!
       -y2**2*y3*y5/2.0_ark+y2*y3**2*y5/2.0_ark-y1*y3**2*y5/2.0_ark+y1*y2**2*y5&!
       /2.0_ark-y1**2*y3*y5)*fea1125
       !
      v4 = s2+(y2*y3**3+y1**3*y3+y1**3*y2+y1*y2**3+y1*y3**3+y2**3*y3)*fea1112+&!
       (y2**2*y3**2+y1**2*y3**2+y1**2*y2**2)*fea1122+(y1*y2**2*y3&!
       +y1**2*y2*y3+y1*y2*y3**2)*fea1123+(5.0_ark/8.0_ark*y2*y4*y5**2+sqrt(3.0_ark)*&!
       y2*y5**3/8.0_ark-sqrt(3.0_ark)*y3*y4**2*y5/8.0_ark+sqrt(3.0_ark)*y2*y4**2*y5/8.0_ark&!
       -3.0_ark/8.0_ark*y2*y4**3+y1*y4*y5**2-sqrt(3.0_ark)*y3*y5**3/8.0_ark&!
       +5.0_ark/8.0_ark*y3*y4*y5**2-3.0_ark/8.0_ark*y3*y4**3)*fea1455

      v5 = 0.0_ark ;  v6 = 0.0_ark
      !
      if (parmax>112) then
        !
        s3 = (y4**5-2.0_ark*y4**3*y5**2-3.0_ark*y4*y5**4)*fea44444+(-4.0_ark*y3*y4*&!
        y5**3*sqrt(3.0_ark)+9.0_ark*y1*y4**2*y5**2-3.0_ark/2.0_ark*y1*y4**4+4.0_ark*y2*y4&!
        *y5**3*sqrt(3.0_ark)+3.0_ark*y2*y4**4+5.0_ark/2.0_ark*y1*y5**4+3.0_ark*y3*y4**4+&!
        y2*y5**4+y3*y5**4)*fea25555+(-y2*y4**4+y3*y4**2*y5**2-2.0_ark*y2*y4*y5&!
        **3*sqrt(3.0_ark)-y3*y4**4-7.0_ark/2.0_ark*y1*y4**2*y5**2-3.0_ark/4.0_ark*y1*y5**4&!
        +2.0_ark*y3*y4*y5**3*sqrt(3.0_ark)+y2*y4**2*y5**2+5.0_ark/4.0_ark*y1*y4**4)*fea24455
        !
        s2 = s3+(y2*y4**3*y5-3.0_ark*y3*y4*y5**3+2.0_ark/3.0_ark*y3*y4**4*sqrt(3.0_ark&!
        )+3.0_ark/4.0_ark*y1*y5**4*sqrt(3.0_ark)+3.0_ark*y2*y4*y5**3-&!
        7.0_ark/12.0_ark*y1*y4**4*sqrt(3.0_ark)+3.0_ark/2.0_ark*y1*y4**2*y5**2*sqrt(3.0_ark)-y3*y4**3*y5&!
        +2.0_ark/3.0_ark*y2*y4**4*sqrt(3.0_ark))*fea24445+(-y2**2*y5**3+y3**2*y4**2*y5+ &!
        y3**2*y5**3+4.0_ark/9.0_ark*y2**2*y4**3*sqrt(3.0_ark)-5.0_ark/9.0_ark*y1**2*y4**3*&!
        sqrt(3.0_ark)+4.0_ark/9.0_ark*y3**2*y4**3*sqrt(3.0_ark)-y2**2*y4**2*y5&!
        -y1**2*y4*y5**2*sqrt(3.0_ark))*fea33445+(y3**2*y4*y5**2-y1**2*y4**3/3.0_ark&!
        -y3**2*y4**3/3.0_ark+y1**2*y4*y5**2+y2**2*y4*y5**2-y2**2*y4**3/3.0_ark)*fea33455
        !
        s1 = s2+(-y2**3*y4*y5+y3**3*y4*y5+y2**3*y5**2*sqrt(3.0_ark)/3.0_ark+y1**&!
        3*y4**2*sqrt(3.0_ark)/2.0_ark+y3**3*y5**2*sqrt(3.0_ark)/3.0_ark- &!
        y1**3*y5**2*sqrt(3.0_ark)/6.0_ark)*fea33345+(y3**3*y4**2+y3**3*y5**2+y2**3*y4**2+y2**3&!
        *y5**2+y1**3*y5**2+y1**3*y4**2)*fea33344+(y3**4*y4+sqrt(3.0_ark)*y3**&!
        4*y5+y2**4*y4-2.0_ark*y1**4*y4-sqrt(3.0_ark)*y2**4*y5)*fea33334+(y2**5+ &!
        y3**5+y1**5)*fea33333+(-4.0_ark/9.0_ark*y1*y2*y4**3*sqrt(3.0_ark)-y1*y2*y5**3+ &!
        y1*y3*y4**2*y5+y2*y3*y4*y5**2*sqrt(3.0_ark)-y1*y2*y4**2*y5+5.0_ark/9.0_ark&!
        *y2*y3*y4**3*sqrt(3.0_ark)-4.0_ark/9.0_ark*y1*y3*y4**3*sqrt(3.0_ark)+y1*y3*y5&!
        **3)*fea13445+(y2*y3*y4*y5**2+y1*y2*y4*y5**2-y2*y3*y4**3/3.0_ark- &!
        y1*y2*y4**3/3.0_ark-y1*y3*y4**3/3.0_ark+y1*y3*y4*y5**2)*fea13455

        s3 = s1+(y1**2*y3*y5**2+y2**2*y3*y4**2+y2**2*y3*y5**2+y1*y2**2*y5**2+&!
        y1**2*y2*y5**2+y1*y2**2*y4**2+y2*y3**2*y4**2+y1*y3**2*y4**2+&!
        y1**2*y3*y4**2+y1**2*y2*y4**2+y1*y3**2*y5**2+y2*y3**2*y5**2)*fea11255&!
        +(2.0_ark/3.0_ark*y1**2*y3*y4**2*sqrt(3.0_ark)+y1*y3**2*y5**2*sqrt(3.0_ark)/2.0_ark+&!
        y1*y2**2*y5**2*sqrt(3.0_ark)/2.0_ark+y2**2*y3*y5**2*sqrt(3.0_ark)/2.0_ark-&!
        y1*y2**2*y4*y5+y2*y3**2*y4*y5+y1*y3**2*y4*y5-y2**2*y3*y4*y5+y2*y3**&!
        2*y4**2*sqrt(3.0_ark)/6.0_ark+y1*y3**2*y4**2*sqrt(3.0_ark)/6.0_ark+y1*y2**2*y4&!
        **2*sqrt(3.0_ark)/6.0_ark+2.0_ark/3.0_ark*y1**2*y2*y4**2*sqrt(3.0_ark)+&!
        y2*y3**2*y5**2*sqrt(3.0_ark)/2.0_ark+y2**2*y3*y4**2*sqrt(3.0_ark)/6.0_ark)*fea13345
        s4 = s3+(y1**2*y2*y4*y5+y1**2*y3*y4**2*sqrt(3.0_ark)/3.0_ark+y1**2*y2*y4&!
        **2*sqrt(3.0_ark)/3.0_ark-y1*y2**2*y4**2*sqrt(3.0_ark)/6.0_ark+y2*y3**2*y4*y5-&!
        y2**2*y3*y4*y5-y1**2*y3*y4*y5+y2*y3**2*y4**2*sqrt(3.0_ark)/3.0_ark+y1*y2&!
        **2*y5**2*sqrt(3.0_ark)/2.0_ark-y1*y3**2*y4**2*sqrt(3.0_ark)/6.0_ark+y2**2*y3*&!
        y4**2*sqrt(3.0_ark)/3.0_ark+y1*y3**2*y5**2*sqrt(3.0_ark)/2.0_ark)*fea11245
        s2 = s4+(-y1**3*y2*y5+y1**3*y3*y5+y2**3*y3*y5/2.0_ark-y1*y2**3*y4*sqrt(3.0_ark)/2.0_ark-&!
        y1*y2**3*y5/2.0_ark-y2*y3**3*y5/2.0_ark+y1*y3**3*y5/2.0_ark+y2&!
        **3*y3*y4*sqrt(3.0_ark)/2.0_ark+y2*y3**3*y4*sqrt(3.0_ark)/2.0_ark-y1*y3**3*y4*&!
        sqrt(3.0_ark)/2.0_ark)*fea11135+(y1**3*y3*y4-y2**3*y3*y4/2.0_ark+y1**3*y2*y4-&!
        y2*y3**3*y4/2.0_ark-y1*y3**3*y4/2.0_ark+y1*y2**3*y5*sqrt(3.0_ark)/2.0_ark+y2&!
        **3*y3*y5*sqrt(3.0_ark)/2.0_ark-y2*y3**3*y5*sqrt(3.0_ark)/2.0_ark-y1*y2**3*y4/&!
        2.0_ark-y1*y3**3*y5*sqrt(3.0_ark)/2.0_ark)*fea11134

        v5 = s2+(y1*y2**4+y1**4*y3+y1**4*y2+y2**4*y3+y2*y3**4+y1*y3**4)*fea23333+&!
        (-2.0_ark*y2**2*y3**2*y4+y1**2*y2**2*y4-sqrt(3.0_ark)*y1**2*y3**2&!
        *y5+sqrt(3.0_ark)*y1**2*y2**2*y5+y1**2*y3**2*y4)*fea11334+(y1**2*y3**&!
        3+y1**3*y3**2+y2**2*y3**3+y1**2*y2**3+y1**3*y2**2+y2**3*y3**2)*fea11333+&!
        (y1*y2*y3*y4**2+y1*y2*y3*y5**2)*fea12355+(-y1*y2*y3**2*y4/2.0_ark-&!
        y1*y2**2*y3*y4/2.0_ark-sqrt(3.0_ark)*y1*y2*y3**2*y5/2.0_ark+y1**2*y2*y3*y4+&!
        sqrt(3.0_ark)*y1*y2**2*y3*y5/2.0_ark)*fea11234+(y1*y2**3*y3+y1*y2*y3**3+&!
        y1**3*y2*y3)*fea11123+(y1**2*y2**2*y3+y1*y2**2*y3**2+y1**2*y2*y3**2)*fea11233

        s3 = (y2**3*y4**3*sqrt(3.0_ark)-y2**3*y4**2*y5+y3**3*y4**2*y5-&!
        5.0_ark/3.0_ark*y2**3*y4*y5**2*sqrt(3.0_ark)+y3**3*y4**3*sqrt(3.0_ark)-5.0_ark/3.0_ark*y3**&!
        3*y4*y5**2*sqrt(3.0_ark)-y2**3*y5**3+y3**3*y5**3-8.0_ark/3.0_ark*y1**3*y4*y5**2*sqrt(3.0_ark))*fea333555+&!
        (y1**4*y5**2*sqrt(3.0_ark)/2.0_ark+y2**4*y4*y5+y2**4*y4**2*sqrt(3.0_ark)/3.0_ark+&!
        y3**4*y4**2*sqrt(3.0_ark)/3.0_ark-y3**4*y4&!
        *y5-y1**4*y4**2*sqrt(3.0_ark)/6.0_ark)*fea222245+(y1*y3**5+y1*y2**5+y2**&!
        5*y3+y1**5*y3+y1**5*y2+y2*y3**5)*fea133333+(y1**4*y3*y4-2.0_ark*y2**4&!
        *y3*y4+y1**4*y2*y4+y1*y2**4*y5*sqrt(3.0_ark)+y1*y3**4*y4-2.0_ark*y2*y3**&!
        4*y4+y1**4*y2*y5*sqrt(3.0_ark)-y1*y3**4*y5*sqrt(3.0_ark)-y1**4*y3*y5*sqrt(3.0_ark)+&!
        y1*y2**4*y4)*fea133334+(-y1*y2*y3*y4**3/3.0_ark+y1*y2*y3*y4*y5**2)*fea123455

        s4 = s3+(2.0_ark/3.0_ark*sqrt(3.0_ark)*y1*y2**2*y3**2*y4-y1**2*y2**2*y3*y5-&!
        sqrt(3.0_ark)*y1**2*y2**2*y3*y4/3.0_ark+y1**2*y2*y3**2*y5-&!
        sqrt(3.0_ark)*y1**2*y2*y3**2*y4/3.0_ark)*fea112335+(y1*y2**2*y3*y5**2+y1*y2*y3**2*y5**2+&!
        y1*y2*y3**2*y4**2+y1*y2**2*y3*y4**2+y1**2*y2*y3*y4**2+y1**2*y2*y3*y5**2)*fea112355

        s2 = s4+(y2**3*y3**2*y5-y1**3*y2**2*y5/2.0_ark-y1**2*y3**3*y5/2.0_ark-y2&!
        **2*y3**3*y5+y1**3*y2**2*y4*sqrt(3.0_ark)/2.0_ark-y1**2*y2**3*y4*sqrt(3.0_ark)/2.0_ark+&!
        y1**3*y3**2*y5/2.0_ark+y1**2*y2**3*y5/2.0_ark+y1**3*y3**2*y4*sqrt(3.0_ark)/2.0_ark-&!
        y1**2*y3**3*y4*sqrt(3.0_ark)/2.0_ark)*fea222335+(-y1**2*y2&!
        **2*y5**2*sqrt(3.0_ark)/2.0_ark-y1**2*y3**2*y5**2*sqrt(3.0_ark)/2.0_ark-y1**2*&!
        y2**2*y4**2*sqrt(3.0_ark)/6.0_ark-y1**2*y2**2*y4*y5-2.0_ark/3.0_ark*y2**2*y3**&!
        2*y4**2*sqrt(3.0_ark)+y1**2*y3**2*y4*y5-y1**2*y3**2*y4**2*sqrt(3.0_ark)/&!
        6.0_ark)*fea113345+(y2**2*y3**2*y5**2+y2**2*y3**2*y4**2+y1**2*y2**2*y5**2+&!
        y1**2*y3**2*y4**2+y1**2*y3**2*y5**2+y1**2*y2**2*y4**2)*fea223355

        s3 = s2+(y1*y2*y3**2*y4**2*sqrt(3.0_ark)/6.0_ark+y1*y2*y3**2*y4*y5+y1*y2&!
        *y3**2*y5**2*sqrt(3.0_ark)/2.0_ark+2.0_ark/3.0_ark*y1**2*y2*y3*y4**2*sqrt(3.0_ark&!
        )-y1*y2**2*y3*y4*y5+y1*y2**2*y3*y4**2*sqrt(3.0_ark)/6.0_ark+y1*y2**2*y3*&!
        y5**2*sqrt(3.0_ark)/2.0_ark)*fea123345+(-y1**3*y2**2*y5*sqrt(3.0_ark)/2.0_ark-&!
        y1**3*y2**2*y4/2.0_ark-y1**3*y3**2*y4/2.0_ark-y1**2*y2**3*y4/2.0_ark+y1**3*&!
        y3**2*y5*sqrt(3.0_ark)/2.0_ark-y1**2*y3**3*y4/2.0_ark+y2**3*y3**2*y4-y1**2*&!
        y2**3*y5*sqrt(3.0_ark)/2.0_ark+y2**2*y3**3*y4+y1**2*y3**3*y5*sqrt(3.0_ark)/&!
        2.0_ark)*fea222334+(3.0_ark*y3**2*y4**4+5.0_ark/2.0_ark*y1**2*y5**4+y2**2*y5**&!
        4+3.0_ark*y2**2*y4**4-4.0_ark*y3**2*y4*y5**3*sqrt(3.0_ark)+y3**2*y5**4+9.0_ark&!
        *y1**2*y4**2*y5**2-3.0_ark/2.0_ark*y1**2*y4**4+4.0_ark*y2**2*y4*y5**3*sqrt(&!
        3.0_ark))*fea335555+(y1**3*y2**3+y1**3*y3**3+y2**3*y3**3)*fea222333

        s4 = s3+(y3*y4**5/5.0_ark-y2*y4**4*y5*sqrt(3.0_ark)/2.0_ark-2.0_ark/5.0_ark*y1*y4&!
        **5-2.0_ark*y1*y4**3*y5**2-3.0_ark/10.0_ark*y2*y5**5*sqrt(3.0_ark)+y3*y4**3*y5&!
        **2+y3*y4**4*y5*sqrt(3.0_ark)/2.0_ark+y2*y4**3*y5**2+3.0_ark/10.0_ark*y3*y5**5&!
        *sqrt(3.0_ark)+y2*y4**5/5.0_ark)*fea244455+(y2**5*y4-2.0_ark*y1**5*y4-sqrt(&!
        3.0_ark)*y2**5*y5+y3**5*y4+sqrt(3.0_ark)*y3**5*y5)*fea222224

        s5 = s4+(-y3*y5**5*sqrt(3.0_ark)/5.0_ark+y2*y5**5*sqrt(3.0_ark)/5.0_ark+y1*y4*&!
        y5**4-7.0_ark/15.0_ark*y2*y4**5+y2*y4**4*y5*sqrt(3.0_ark)/3.0_ark-y3*y4**4*y5*&!
        sqrt(3.0_ark)/3.0_ark+y3*y4*y5**4+y2*y4*y5**4+2.0_ark*y1*y4**3*y5**2-7.0_ark/15.0_ark*y3*y4**5-&!
        y1*y4**5/15.0_ark)*fea145555

        s1 = s5+(-sqrt(3.0_ark)*y1*y2*y3**3*y5/2.0_ark+y1**3*y2*y3*y4+sqrt(3.0_ark)&!
        *y1*y2**3*y3*y5/2.0_ark-y1*y2**3*y3*y4/2.0_ark-y1*y2*y3**3*y4/2.0_ark)*fea111234+&!
        (y3*y4**4*y5/3.0_ark+y3*y4**5*sqrt(3.0_ark)/18.0_ark-y2*y4**4*y5/3.0_ark&!
        -y2*y4*y5**4*sqrt(3.0_ark)/2.0_ark-y3*y4**2*y5**3+2.0_ark/9.0_ark*y1*y4**5*sqrt(3.0_ark)+&!
        y2*y4**5*sqrt(3.0_ark)/18.0_ark+y2*y4**2*y5**3-2.0_ark/3.0_ark*y1*y4**&!
        3*y5**2*sqrt(3.0_ark)-y3*y4*y5**4*sqrt(3.0_ark)/2.0_ark)*fea244555+(y1*y2*y4**2*y5**2-&!
        3.0_ark/4.0_ark*y2*y3*y4**4-y1*y2*y5**4-y1*y3*y5**4+5.0_ark/4.0_ark*y2*y3*y5**4+&!
        y1*y3*y4**2*y5**2-7.0_ark/2.0_ark*y2*y3*y4**2*y5**2-2.0_ark*y1&!
        *y2*y4**3*y5*sqrt(3.0_ark)+2.0_ark*y1*y3*y4**3*y5*sqrt(3.0_ark))*fea124455

        s3 = s1+(y2**6+y1**6+y3**6)*fea333333+(y1*y2**4*y3+y1**4*y2*y3+y1*&!
        y2*y3**4)*fea111123+fea112233*y1**2*y2**2*y3**2+(y1**4*y4**2+y2**4&!
        *y4**2+y2**4*y5**2+y3**4*y4**2+y1**4*y5**2+y3**4*y5**2)*fea222255
        s4 = s3+(3.0_ark*y1*y3*y5**4+y1*y3*y4**4+9.0_ark*y2*y3*y4**2*y5**2-3.0_ark/&!
        2.0_ark*y2*y3*y5**4-4.0_ark*y1*y3*y4**3*y5*sqrt(3.0_ark)+y1*y2*y4**4+&!
        4.0_ark*y1*y2*y4**3*y5*sqrt(3.0_ark)+3.0_ark*y1*y2*y5**4+5.0_ark/2.0_ark*y2*y3*y4**4)*fea134444+&!
        (-y1*y3**2*y5**3*sqrt(3.0_ark)/3.0_ark-7.0_ark/3.0_ark*y1**2*y3*y4*y5&!
        **2+5.0_ark/3.0_ark*y1*y2**2*y4**2*y5*sqrt(3.0_ark)-13.0_ark/3.0_ark*y2**2*y3*y4*&!
        y5**2-4.0_ark/3.0_ark*y2*y3**2*y5**3*sqrt(3.0_ark)-7.0_ark/3.0_ark*y1**2*y2*y4*y5&!
        **2-16.0_ark/3.0_ark*y1*y3**2*y4*y5**2+4.0_ark/3.0_ark*y1**2*y3*y4**2*y5*sqrt(&!
        3.0_ark)+4.0_ark/3.0_ark*y2**2*y3*y5**3*sqrt(3.0_ark)+3.0_ark*y1**2*y2*y4**3+&!
        y2*y3**2*y4**3+y1*y2**2*y5**3*sqrt(3.0_ark)/3.0_ark+y2**2*y3*y4**3-13.0_ark/3.0_ark&!
        *y2*y3**2*y4*y5**2-5.0_ark/3.0_ark*y1*y3**2*y4**2*y5*sqrt(3.0_ark)-&!
        4.0_ark/3.0_ark*y1**2*y2*y4**2*y5*sqrt(3.0_ark)+3.0_ark*y1**2*y3*y4**3-16.0_ark/3.0_ark*y1*&!
        y2**2*y4*y5**2)*fea233444

        s5 = s4+(2.0_ark*y1*y3**2*y5**3+4.0_ark*y2*y3**2*y5**3+4.0_ark*y2**2*y3*y4*&!
        y5**2*sqrt(3.0_ark)-2.0_ark*y1*y2**2*y5**3+y1**2*y3*y4*y5**2*sqrt(3.0_ark)+&!
        6.0_ark*y1*y3**2*y4**2*y5-6.0_ark*y1*y2**2*y4**2*y5-3.0_ark*y1**2*y3*y4**2*&!
        y5+y1**2*y2*y4*y5**2*sqrt(3.0_ark)+4.0_ark*y1*y3**2*y4*y5**2*sqrt(3.0_ark)-&!
        3.0_ark*y1**2*y2*y4**3*sqrt(3.0_ark)-4.0_ark*y2**2*y3*y5**3+3.0_ark*y1**2*y2*y4**2*y5-&!
        y1**2*y2*y5**3+y1**2*y3*y5**3-3.0_ark*y1**2*y3*y4**3*sqrt(3.0_ark&!
        )+4.0_ark*y2*y3**2*y4*y5**2*sqrt(3.0_ark)+4.0_ark*y1*y2**2*y4*y5**2*sqrt(3.0_ark))*fea113555

        s2 = s5+(-2.0_ark/3.0_ark*y3**2*y4**4*sqrt(3.0_ark)-3.0_ark/2.0_ark*y1**2*y4**2*y5**2*sqrt(3.0_ark)-&!
        3.0_ark/4.0_ark*y1**2*y5**4*sqrt(3.0_ark)-y2**2*y4**3*y5+&!
        7.0_ark/12.0_ark*y1**2*y4**4*sqrt(3.0_ark)+y3**2*y4**3*y5+3.0_ark*y3**2*y4*y5**3&!
        -2.0_ark/3.0_ark*y2**2*y4**4*sqrt(3.0_ark)-3.0_ark*y2**2*y4*y5**3)*fea334445+(&!
        -3.0_ark*y1*y3*y4**3*y5+2.0_ark/3.0_ark*y1*y2*y5**4*sqrt(3.0_ark)-y1*y3*y4*y5**3+&!
        2.0_ark/3.0_ark*y1*y3*y5**4*sqrt(3.0_ark)+3.0_ark*y1*y2*y4**3*y5-7.0_ark/12.0_ark&!
        *y2*y3*y5**4*sqrt(3.0_ark)+3.0_ark/2.0_ark*y2*y3*y4**2*y5**2*sqrt(3.0_ark)+y1*&!
        y2*y4*y5**3+3.0_ark/4.0_ark*y2*y3*y4**4*sqrt(3.0_ark))*fea124555+(2.0_ark*y3**&!
        2*y4*y5**3*sqrt(3.0_ark)-7.0_ark/2.0_ark*y1**2*y4**2*y5**2+y2**2*y4**2*y5**&!
        2-y2**2*y4**4-y3**2*y4**4-2.0_ark*y2**2*y4*y5**3*sqrt(3.0_ark)-3.0_ark/4.0_ark&!
        *y1**2*y5**4+5.0_ark/4.0_ark*y1**2*y4**4+y3**2*y4**2*y5**2)*fea334455
        s3 = s2+(-6.0_ark*y4**2*y5**4+9.0_ark*y4**4*y5**2+y5**6)*fea555555+(y2*y3**3*y4**2+&!
        y2*y3**3*y5**2+y1*y3**3*y4**2+y1*y2**3*y4**2+y1**3*y2*y4**2+&!
        y1*y2**3*y5**2+y1**3*y3*y5**2+y1**3*y3*y4**2+y1**3*y2*y5**2+y2**3*y3*y4**2+&!
        y1*y3**3*y5**2+y2**3*y3*y5**2)*fea233344+(y1*y2**3*y5**2*sqrt(3.0_ark)/6.0_ark-&!
        y2**3*y3*y5**2*sqrt(3.0_ark)/3.0_ark-y2*y3**3*y5**2&!
        *sqrt(3.0_ark)/3.0_ark+y1**3*y2*y4*y5-y1**3*y2*y5**2*sqrt(3.0_ark)/3.0_ark-&!
        y1**3*y3*y4*y5-y1**3*y3*y5**2*sqrt(3.0_ark)/3.0_ark-y1*y3**3*y4**2*sqrt(3.0_ark&!
        )/2.0_ark+y1*y3**3*y5**2*sqrt(3.0_ark)/6.0_ark-y2**3*y3*y4*y5+y2*y3**3*y4*&!
        y5-y1*y2**3*y4**2*sqrt(3.0_ark)/2.0_ark)*fea233345+(-3.0_ark*y2**3*y4*y5**2&!
        +y3**3*y4**3-3.0_ark*y3**3*y4*y5**2-3.0_ark*y1**3*y4*y5**2+y2**3*y4**3+&!
        y1**3*y4**3)*fea111444+(y1*y2**3*y3**2+y1**3*y2**2*y3+y1**2*y2**3*y3+&!
        y1*y2**2*y3**3+y1**2*y2*y3**3+y1**3*y2*y3**2)*fea111233

        s4 = s3+(9.0_ark*y4**2*y5**4-6.0_ark*y4**4*y5**2+y4**6)*fea444444+(-5.0_ark&!
        /3.0_ark*y1*y2**2*y4**2*y5*sqrt(3.0_ark)+y1*y2**2*y4**3-4.0_ark/3.0_ark*y1**2*&!
        y3*y4**2*y5*sqrt(3.0_ark)-2.0_ark*y1**2*y2*y4**3-y1*y2**2*y5**3*sqrt(3.0_ark&!
        )/3.0_ark+4.0_ark/3.0_ark*y2**2*y3*y4*y5**2-4.0_ark/3.0_ark*y2**2*y3*y5**3*sqrt(&!
        3.0_ark)-2.0_ark*y1**2*y3*y4**3+7.0_ark/3.0_ark*y1*y2**2*y4*y5**2-2.0_ark/3.0_ark*y1&!
        **2*y3*y4*y5**2+y1*y3**2*y4**3+4.0_ark/3.0_ark*y2*y3**2*y5**3*sqrt(3.0_ark)&!
        +y1*y3**2*y5**3*sqrt(3.0_ark)/3.0_ark+4.0_ark/3.0_ark*y1**2*y2*y4**2*y5*sqrt(3.0_ark)&!
        +4.0_ark/3.0_ark*y2*y3**2*y4*y5**2+5.0_ark/3.0_ark*y1*y3**2*y4**2*y5*sqrt(&!
        3.0_ark)-2.0_ark/3.0_ark*y1**2*y2*y4*y5**2+7.0_ark/3.0_ark*y1*y3**2*y4*y5**2)*fea133444

        s5 = s4+(-y1**3*y2*y4*y5+2.0_ark/3.0_ark*y2**3*y3*y5**2*sqrt(3.0_ark)+y1*y3&!
        **3*y4**2*sqrt(3.0_ark)/2.0_ark+y1**3*y3*y4**2*sqrt(3.0_ark)/2.0_ark+y1**3*y3*&!
        y5**2*sqrt(3.0_ark)/6.0_ark+y1**3*y2*y5**2*sqrt(3.0_ark)/6.0_ark+y1**3*y3*y4*y5+&!
        y1*y2**3*y5**2*sqrt(3.0_ark)/6.0_ark+y1**3*y2*y4**2*sqrt(3.0_ark)/2.0_ark+&!
        2.0_ark/3.0_ark*y2*y3**3*y5**2*sqrt(3.0_ark)-y1*y2**3*y4*y5+y1*y2**3*y4**2*sqrt(3.0_ark)/2.0_ark+&!
        y1*y3**3*y5**2*sqrt(3.0_ark)/6.0_ark+y1*y3**3*y4*y5)*fea133345

        v6 = s5+(-y2**2*y3*y4**2*y5+y1**2*y3*y4*y5**2*sqrt(3.0_ark)/3.0_ark+y2*y3**2*y4**2*y5+&!
        y2*y3**2*y5**3-y1*y2**2*y5**3+4.0_ark/3.0_ark*y2**2*y3*y4*&!
        y5**2*sqrt(3.0_ark)+4.0_ark/3.0_ark*y2*y3**2*y4*y5**2*sqrt(3.0_ark)-y1*y2**2*y4**2*y5+&!
        4.0_ark/3.0_ark*y1*y3**2*y4*y5**2*sqrt(3.0_ark)-y2**2*y3*y5**3+y1*y3**2*y5**3+&!
        y1**2*y2*y4*y5**2*sqrt(3.0_ark)/3.0_ark-y1**2*y2*y4**3*sqrt(3.0_ark)&!
        +y1*y3**2*y4**2*y5-y1**2*y3*y4**3*sqrt(3.0_ark)+4.0_ark/3.0_ark*y1*y2**&!
        2*y4*y5**2*sqrt(3.0_ark))*fea233445+(y2*y3**4*y4*sqrt(3.0_ark)-y1**4*y2*&!
        y5+y2**4*y3*y4*sqrt(3.0_ark)-y1**4*y3*y4*sqrt(3.0_ark)+y2*y3**4*y5-2.0_ark*&!
        y1*y2**4*y5+2.0_ark*y1*y3**4*y5-y1**4*y2*y4*sqrt(3.0_ark)+y1**4*y3*y5-y2&!
        **4*y3*y5)*fea233335+(y2**2*y3**4+y1**4*y3**2+y1**2*y2**4+y2**4*y3&!
        **2+y1**2*y3**4+y1**4*y2**2)*fea222233
        !
       endif

       f = ( v0+(v1+v2+v3+v4+v5+v6) )

end function xy3_dms_A


function xy3_dms_E(ix, parmax, param, local) result(f)

  integer, intent(in)   ::  ix, parmax
  real(ark), intent(in) ::  param(parmax)
  real(ark), intent(in) ::  local(6)
  real(ark)             ::  f

  real(ark) r14, r24, r34, alpha1, alpha2, alpha3, rhobar
  real(ark) y1,y2,y3,y4,y5,alpha,sinrho,cosrho,coro
  real(ark) v0,v1,v2,v3,v4,v5,v6
  real(ark) d56,s1,s2,s3,s4,s5,s6,t4x,t4y,t56x,t56y,s4a,s4b,drho
  real(ark) rhoe, re14, beta

  real(ark)          &!
  fea1, fea4, fea11, fea12, fea14, fea24, fea44, fea111, &!
  fea112, fea122, fea114, fea124, fea224, fea144, fea244, fea444, &!
  fea135, fea155, fea1111, fea1112, fea1122, fea1222, fea1123, fea1114, &!
  fea1124, fea1224, fea2224, fea1234, fea1144, fea1244, fea1444, fea2444, &!
  fea4444, fea1125, fea1225, fea1245, fea2445, fea1155, fea1255, fea3355, &!
  fea1455, fea4455, fea11111, fea11112, fea11122, fea11123, fea11223, fea11114, &!
  fea11124, fea11224, fea11234, fea22344, fea11135, fea22235, fea11245, fea12245, &!
  fea22245, fea11155, fea11255, fea12255, fea22255, fea12355, fea11455, fea12455, &!
  fea22455, fea23455, fea24455, fea12555, fea22555, fea24555, fea15555, fea45555, &!
  fea111111, fea111112, fea111122, fea122333, fea123333, fea122334, fea123334, fea223334, &!
  fea133334, fea233334, fea333334, fea111344, fea112344, fea223344, fea133344, fea233344, &!
  fea333344, fea122444, fea222444, fea113444, fea123444, fea124444, fea444444, fea233335, &!
  fea233345, fea113445, fea133445, fea233445, fea235555, fea335555, fea145555, hea11112, &!
  hea11122, hea11244, hea13444, hea33444, hea12335, hea22335, hea13335, hea33335, &!
  hea14445, hea34555, hea55555, hea111333, hea122333, hea223333, hea233333, hea111334, &!
  hea113334, hea223334, hea233334, hea113344, hea123344, hea133344, hea233344, hea133444, &!
  hea233444, hea333444, hea122335, hea123335, hea333335, hea111145, hea123345, hea223345, &!
  hea234445, hea334445, hea144445, hea344445, hea444445, hea333355, hea333455, hea334455, &!
  hea344455, hea222555, hea223555, hea114555, hea234555, hea334555, hea135555, hea345555, hea355555

    real(ark) :: &!
    f0a1     , f1a1     , f2a1     , f3a1  ,f4a1,  &!
    f0a4     , f1a4     , f2a4     , f3a4  ,f4a4,  &!
    f0a11    , f1a11    , f2a11    , f3a11 ,       &!
    f0a12    , f1a12    , f2a12    , f3a12 ,       &!
    f0a14    , f1a14    , f2a14    , f3a14 ,       &!
    f0a24    , f1a24    , f2a24    , f3a24 ,       &!
    f0a44    , f1a44    , f2a44    , f3a44 ,       &!
    f0a111   , f1a111   , f2a111   , f3a111,       &!
    f0a112   , f1a112   , f2a112   , f3a112,       &!
    f0a122   , f1a122   , f2a122   , f3a122,       &!
    f0a114   , f1a114   , f2a114   , f3a114,       &!
    f0a124   , f1a124   , f2a124   , f3a124,       &!
    f0a224   , f1a224   , f2a224   , f3a224,       &!
    f0a144   , f1a144   , f2a144   , f3a144,       &!
    f0a244   , f1a244   , f2a244   , f3a244,       &!
    f0a444   , f1a444   , f2a444   , f3a444,       &!
    f0a135   , f1a135   , f2a135   , f3a135,       &!
    f0a155   , f1a155   , f2a155   , f3a155,       &!
    f0a1111  , f1a1111  , f2a1111  ,               &!
    f0a1112  , f1a1112  , f2a1112  ,               &!
    f0a1122  , f1a1122  , f2a1122  ,               &!
    f0a1222  , f1a1222  , f2a1222  ,               &!
    f0a1123  , f1a1123  , f2a1123  ,               &!
    f0a1114  , f1a1114  , f2a1114  ,               &!
    f0a1124  , f1a1124  , f2a1124  ,               &!
    f0a1224  , f1a1224  , f2a1224  ,               &!
    f0a2224  , f1a2224  , f2a2224  ,               &!
    f0a1234  , f1a1234  , f2a1234  ,               &!
    f0a1144  , f1a1144  , f2a1144  ,               &!
    f0a1244  , f1a1244  , f2a1244  ,               &!
    f0a1444  , f1a1444  , f2a1444  ,               &!
    f0a2444  , f1a2444  , f2a2444  ,               &!
    f0a4444  , f1a4444  , f2a4444  ,               &!
    f0a1125  , f1a1125  , f2a1125  ,               &!
    f0a1225  , f1a1225  , f2a1225  ,               &!
    f0a1245  , f1a1245  , f2a1245  ,               &!
    f0a2445  , f1a2445  , f2a2445  ,               &!
    f0a1155  , f1a1155  , f2a1155  ,               &!
    f0a1255  , f1a1255  , f2a1255  ,               &!
    f0a3355  , f1a3355  , f2a3355  ,               &!
    f0a1455  , f1a1455  , f2a1455  ,               &!
    f0a4455  , f1a4455  , f2a4455  ,               &!
    f0a11111 , f1a11111 ,                          &!
    f0a11112 , f1a11112 ,                          &!
    f0a11122 , f1a11122 ,                          &!
    f0a11123 , f1a11123 ,                          &!
    f0a11223 , f1a11223 ,                          &!
    f0a11114 , f1a11114 ,                          &!
    f0a11124 , f1a11124 ,                          &!
    f0a11224 , f1a11224 ,                          &!
    f0a11234 , f1a11234 ,                          &!
    f0a22344 , f1a22344 ,                          &!
    f0a11135 , f1a11135 ,                          &!
    f0a22235 , f1a22235 ,                          &!
    f0a11245 , f1a11245 ,                          &!
    f0a12245 , f1a12245 ,                          &!
    f0a22245 , f1a22245 ,                          &!
    f0a11155 , f1a11155 ,                          &!
    f0a11255 , f1a11255 ,                          &!
    f0a12255 , f1a12255 ,                          &!
    f0a22255 , f1a22255 ,                          &!
    f0a12355 , f1a12355 ,                          &!
    f0a11455 , f1a11455 ,                          &!
    f0a12455 , f1a12455 ,                          &!
    f0a22455 , f1a22455 ,                          &!
    f0a23455 , f1a23455 ,                          &!
    f0a24455 , f1a24455 ,                          &!
    f0a12555 , f1a12555 ,                          &!
    f0a22555 , f1a22555 ,                          &!
    f0a24555 , f1a24555 ,                          &!
    f0a15555 , f1a15555 ,                          &!
    f0a45555 , f1a45555 ,                          &!
    f0a111111, f1a111111,                          &!
    f0a111112, f1a111112,                          &!
    f0a111122, f1a111122,                          &!
    f0a122333, f1a122333,                          &!
    f0a123333, f1a123333,                          &!
    f0a122334, f1a122334,                          &!
    f0a123334, f1a123334,                          &!
    f0a223334, f1a223334,                          &!
    f0a133334, f1a133334,                          &!
    f0a233334, f1a233334,                          &!
    f0a333334, f1a333334,                          &!
    f0a111344, f1a111344,                          &!
    f0a112344, f1a112344,                          &!
    f0a223344, f1a223344,                          &!
    f0a133344, f1a133344,                          &!
    f0a233344, f1a233344,                          &!
    f0a333344, f1a333344,                          &!
    f0a122444, f1a122444,                          &!
    f0a222444, f1a222444,                          &!
    f0a113444, f1a113444,                          &!
    f0a123444, f1a123444,                          &!
    f0a124444, f1a124444,                          &!
    f0a444444, f1a444444,                          &!
    f0a233335, f1a233335,                          &!
    f0a233345, f1a233345,                          &!
    f0a113445, f1a113445,                          &!
    f0a133445, f1a133445,                          &!
    f0a233445, f1a233445,                          &!
    f0a235555, f1a235555,                          &!
    f0a335555, f1a335555,                          &!
    f0a145555, f1a145555,                          &!
    h0a11112 , h1a11112 ,                          &!
    h0a11122 , h1a11122 ,                          &!
    h0a11244 , h1a11244 ,                          &!
    h0a13444 , h1a13444 ,                          &!
    h0a33444 , h1a33444 ,                          &!
    h0a12335 , h1a12335 ,                          &!
    h0a22335 , h1a22335 ,                          &!
    h0a13335 , h1a13335 ,                          &!
    h0a33335 , h1a33335 ,                          &!
    h0a14445 , h1a14445 ,                          &!
    h0a34555 , h1a34555 ,                          &!
    h0a55555 , h1a55555 ,                          &!
    h0a111333, h1a111333,                          &!
    h0a122333, h1a122333,                          &!
    h0a223333, h1a223333,                          &!
    h0a233333, h1a233333,                          &!
    h0a111334, h1a111334,                          &!
    h0a113334, h1a113334,                          &!
    h0a223334, h1a223334,                          &!
    h0a233334, h1a233334,                          &!
    h0a113344, h1a113344,                          &!
    h0a123344, h1a123344,                          &!
    h0a133344, h1a133344,                          &!
    h0a233344, h1a233344,                          &!
    h0a133444, h1a133444,                          &!
    h0a233444, h1a233444,                          &!
    h0a333444, h1a333444,                          &!
    h0a122335, h1a122335,                          &!
    h0a123335, h1a123335,                          &!
    h0a333335, h1a333335,                          &!
    h0a111145, h1a111145,                          &!
    h0a123345, h1a123345,                          &!
    h0a223345, h1a223345,                          &!
    h0a234445, h1a234445,                          &!
    h0a334445, h1a334445,                          &!
    h0a144445, h1a144445,                          &!
    h0a344445, h1a344445,                          &!
    h0a444445, h1a444445,                          &!
    h0a333355, h1a333355,                          &!
    h0a333455, h1a333455,                          &!
    h0a334455, h1a334455,                          &!
    h0a344455, h1a344455,                          &!
    h0a222555, h1a222555,                          &!
    h0a223555, h1a223555,                          &!
    h0a114555, h1a114555,                          &!
    h0a234555, h1a234555,                          &!
    h0a334555, h1a334555,                          &!
    h0a135555, h1a135555,                          &!
    h0a345555, h1a345555,                          &!
    h0a355555, h1a355555


    rhoe      = param(1) * pi / 180.0_ark
    re14      =  param(  2)
    beta      =  param(  3)**2
    f0a1      =  param(  4)
    f1a1      =  param(  5)
    f2a1      =  param(  6)
    f3a1      =  param(  7)
    f4a1      =  param(  8)
    f0a4      =  param(  9)
    f1a4      =  param( 10)
    f2a4      =  param( 11)
    f3a4      =  param( 12)
    f4a4      =  param( 13)
    f0a11     =  param( 14)
    f1a11     =  param( 15)
    f2a11     =  param( 16)
    f3a11     =  param( 17)
    f0a12     =  param( 18)
    f1a12     =  param( 19)
    f2a12     =  param( 20)
    f3a12     =  param( 21)
    f0a14     =  param( 22)
    f1a14     =  param( 23)
    f2a14     =  param( 24)
    f3a14     =  param( 25)
    f0a24     =  param( 26)
    f1a24     =  param( 27)
    f2a24     =  param( 28)
    f3a24     =  param( 29)
    f0a44     =  param( 30)
    f1a44     =  param( 31)
    f2a44     =  param( 32)
    f3a44     =  param( 33)
    f0a111    =  param( 34)
    f1a111    =  param( 35)
    f2a111    =  param( 36)
    f3a111    =  param( 37)
    f0a112    =  param( 38)
    f1a112    =  param( 39)
    f2a112    =  param( 40)
    f3a112    =  param( 41)
    f0a122    =  param( 42)
    f1a122    =  param( 43)
    f2a122    =  param( 44)
    f3a122    =  param( 45)
    f0a114    =  param( 46)
    f1a114    =  param( 47)
    f2a114    =  param( 48)
    f3a114    =  param( 49)
    f0a124    =  param( 50)
    f1a124    =  param( 51)
    f2a124    =  param( 52)
    f3a124    =  param( 53)
    f0a224    =  param( 54)
    f1a224    =  param( 55)
    f2a224    =  param( 56)
    f3a224    =  param( 57)
    f0a144    =  param( 58)
    f1a144    =  param( 59)
    f2a144    =  param( 60)
    f3a144    =  param( 61)
    f0a244    =  param( 62)
    f1a244    =  param( 63)
    f2a244    =  param( 64)
    f3a244    =  param( 65)
    f0a444    =  param( 66)
    f1a444    =  param( 67)
    f2a444    =  param( 68)
    f3a444    =  param( 69)
    f0a135    =  param( 70)
    f1a135    =  param( 71)
    f2a135    =  param( 72)
    f3a135    =  param( 73)
    f0a155    =  param( 74)
    f1a155    =  param( 75)
    f2a155    =  param( 76)
    f3a155    =  param( 77)
    f0a1111   =  param( 78)
    f1a1111   =  param( 79)
    f2a1111   =  param( 80)
    f0a1112   =  param( 81)
    f1a1112   =  param( 82)
    f2a1112   =  param( 83)
    f0a1122   =  param( 84)
    f1a1122   =  param( 85)
    f2a1122   =  param( 86)
    f0a1222   =  param( 87)
    f1a1222   =  param( 88)
    f2a1222   =  param( 89)
    f0a1123   =  param( 90)
    f1a1123   =  param( 91)
    f2a1123   =  param( 92)
    f0a1114   =  param( 93)
    f1a1114   =  param( 94)
    f2a1114   =  param( 95)
    f0a1124   =  param( 96)
    f1a1124   =  param( 97)
    f2a1124   =  param( 98)
    f0a1224   =  param( 99)
    f1a1224   =  param(100)
    f2a1224   =  param(101)
    f0a2224   =  param(102)
    f1a2224   =  param(103)
    f2a2224   =  param(104)
    f0a1234   =  param(105)
    f1a1234   =  param(106)
    f2a1234   =  param(107)
    f0a1144   =  param(108)
    f1a1144   =  param(109)
    f2a1144   =  param(110)
    f0a1244   =  param(111)
    f1a1244   =  param(112)
    f2a1244   =  param(113)
    f0a1444   =  param(114)
    f1a1444   =  param(115)
    f2a1444   =  param(116)
    f0a2444   =  param(117)
    f1a2444   =  param(118)
    f2a2444   =  param(119)
    f0a4444   =  param(120)
    f1a4444   =  param(121)
    f2a4444   =  param(122)
    f0a1125   =  param(123)
    f1a1125   =  param(124)
    f2a1125   =  param(125)
    f0a1225   =  param(126)
    f1a1225   =  param(127)
    f2a1225   =  param(128)
    f0a1245   =  param(129)
    f1a1245   =  param(130)
    f2a1245   =  param(131)
    f0a2445   =  param(132)
    f1a2445   =  param(133)
    f2a2445   =  param(134)
    f0a1155   =  param(135)
    f1a1155   =  param(136)
    f2a1155   =  param(137)
    f0a1255   =  param(138)
    f1a1255   =  param(139)
    f2a1255   =  param(140)
    f0a3355   =  param(141)
    f1a3355   =  param(142)
    f2a3355   =  param(143)
    f0a1455   =  param(144)
    f1a1455   =  param(145)
    f2a1455   =  param(146)
    f0a4455   =  param(147)
    f1a4455   =  param(148)
    f2a4455   =  param(149)

    if (parmax>149) then
      !
      f0a11111  =  param(150)
      f1a11111  =  param(151)
      f0a11112  =  param(152)
      f1a11112  =  param(153)
      f0a11122  =  param(154)
      f1a11122  =  param(155)
      f0a11123  =  param(156)
      f1a11123  =  param(157)
      f0a11223  =  param(158)
      f1a11223  =  param(159)
      f0a11114  =  param(160)
      f1a11114  =  param(161)
      f0a11124  =  param(162)
      f1a11124  =  param(163)
      f0a11224  =  param(164)
      f1a11224  =  param(165)
      f0a11234  =  param(166)
      f1a11234  =  param(167)
      f0a22344  =  param(168)
      f1a22344  =  param(169)
      f0a11135  =  param(170)
      f1a11135  =  param(171)
      f0a22235  =  param(172)
      f1a22235  =  param(173)
      f0a11245  =  param(174)
      f1a11245  =  param(175)
      f0a12245  =  param(176)
      f1a12245  =  param(177)
      f0a22245  =  param(178)
      f1a22245  =  param(179)
      f0a11155  =  param(180)
      f1a11155  =  param(181)
      f0a11255  =  param(182)
      f1a11255  =  param(183)
      f0a12255  =  param(184)
      f1a12255  =  param(185)
      f0a22255  =  param(186)
      f1a22255  =  param(187)
      f0a12355  =  param(188)
      f1a12355  =  param(189)
      f0a11455  =  param(190)
      f1a11455  =  param(191)
      f0a12455  =  param(192)
      f1a12455  =  param(193)
      f0a22455  =  param(194)
      f1a22455  =  param(195)
      f0a23455  =  param(196)
      f1a23455  =  param(197)
      f0a24455  =  param(198)
      f1a24455  =  param(199)
      f0a12555  =  param(200)
      f1a12555  =  param(201)
      f0a22555  =  param(202)
      f1a22555  =  param(203)
      f0a24555  =  param(204)
      f1a24555  =  param(205)
      f0a15555  =  param(206)
      f1a15555  =  param(207)
      f0a45555  =  param(208)
      f1a45555  =  param(209)
      f0a111111 =  param(210)
      f1a111111 =  param(211)
      f0a111112 =  param(212)
      f1a111112 =  param(213)
      f0a111122 =  param(214)
      f1a111122 =  param(215)
      f0a122333 =  param(216)
      f1a122333 =  param(217)
      f0a123333 =  param(218)
      f1a123333 =  param(219)
      f0a122334 =  param(220)
      f1a122334 =  param(221)
      f0a123334 =  param(222)
      f1a123334 =  param(223)
      f0a223334 =  param(224)
      f1a223334 =  param(225)
      f0a133334 =  param(226)
      f1a133334 =  param(227)
      f0a233334 =  param(228)
      f1a233334 =  param(229)
      f0a333334 =  param(230)
      f1a333334 =  param(231)
      f0a111344 =  param(232)
      f1a111344 =  param(233)
      f0a112344 =  param(234)
      f1a112344 =  param(235)
      f0a223344 =  param(236)
      f1a223344 =  param(237)
      f0a133344 =  param(238)
      f1a133344 =  param(239)
      f0a233344 =  param(240)
      f1a233344 =  param(241)
      f0a333344 =  param(242)
      f1a333344 =  param(243)
      f0a122444 =  param(244)
      f1a122444 =  param(245)
      f0a222444 =  param(246)
      f1a222444 =  param(247)
      f0a113444 =  param(248)
      f1a113444 =  param(249)
      f0a123444 =  param(250)
      f1a123444 =  param(251)
      f0a124444 =  param(252)
      f1a124444 =  param(253)
      f0a444444 =  param(254)
      f1a444444 =  param(255)
      f0a233335 =  param(256)
      f1a233335 =  param(257)
      f0a233345 =  param(258)
      f1a233345 =  param(259)
      f0a113445 =  param(260)
      f1a113445 =  param(261)
      f0a133445 =  param(262)
      f1a133445 =  param(263)
      f0a233445 =  param(264)
      f1a233445 =  param(265)
      f0a235555 =  param(266)
      f1a235555 =  param(267)
      f0a335555 =  param(268)
      f1a335555 =  param(269)
      f0a145555 =  param(270)
      f1a145555 =  param(271)
      h0a11112  =  param(272)
      h1a11112  =  param(273)
      h0a11122  =  param(274)
      h1a11122  =  param(275)
      h0a11244  =  param(276)
      h1a11244  =  param(277)
      h0a13444  =  param(278)
      h1a13444  =  param(279)
      h0a33444  =  param(280)
      h1a33444  =  param(281)
      h0a12335  =  param(282)
      h1a12335  =  param(283)
      h0a22335  =  param(284)
      h1a22335  =  param(285)
      h0a13335  =  param(286)
      h1a13335  =  param(287)
      h0a33335  =  param(288)
      h1a33335  =  param(289)
      h0a14445  =  param(290)
      h1a14445  =  param(291)
      h0a34555  =  param(292)
      h1a34555  =  param(293)
      h0a55555  =  param(294)
      h1a55555  =  param(295)
      h0a111333 =  param(296)
      h1a111333 =  param(297)
      h0a122333 =  param(298)
      h1a122333 =  param(299)
      h0a223333 =  param(300)
      h1a223333 =  param(301)
      h0a233333 =  param(302)
      h1a233333 =  param(303)
      h0a111334 =  param(304)
      h1a111334 =  param(305)
      h0a113334 =  param(306)
      h1a113334 =  param(307)
      h0a223334 =  param(308)
      h1a223334 =  param(309)
      h0a233334 =  param(310)
      h1a233334 =  param(311)
      h0a113344 =  param(312)
      h1a113344 =  param(313)
      h0a123344 =  param(314)
      h1a123344 =  param(315)
      h0a133344 =  param(316)
      h1a133344 =  param(317)
      h0a233344 =  param(318)
      h1a233344 =  param(319)
      h0a133444 =  param(320)
      h1a133444 =  param(321)
      h0a233444 =  param(322)
      h1a233444 =  param(323)
      h0a333444 =  param(324)
      h1a333444 =  param(325)
      h0a122335 =  param(326)
      h1a122335 =  param(327)
      h0a123335 =  param(328)
      h1a123335 =  param(329)
      h0a333335 =  param(330)
      h1a333335 =  param(331)
      h0a111145 =  param(332)
      h1a111145 =  param(333)
      h0a123345 =  param(334)
      h1a123345 =  param(335)
      h0a223345 =  param(336)
      h1a223345 =  param(337)
      h0a234445 =  param(338)
      h1a234445 =  param(339)
      h0a334445 =  param(340)
      h1a334445 =  param(341)
      h0a144445 =  param(342)
      h1a144445 =  param(343)
      h0a344445 =  param(344)
      h1a344445 =  param(345)
      h0a444445 =  param(346)
      h1a444445 =  param(347)
      h0a333355 =  param(348)
      h1a333355 =  param(349)
      h0a333455 =  param(350)
      h1a333455 =  param(351)
      h0a334455 =  param(352)
      h1a334455 =  param(353)
      h0a344455 =  param(354)
      h1a344455 =  param(355)
      h0a222555 =  param(356)
      h1a222555 =  param(357)
      h0a223555 =  param(358)
      h1a223555 =  param(359)
      h0a114555 =  param(360)
      h1a114555 =  param(361)
      h0a234555 =  param(362)
      h1a234555 =  param(363)
      h0a334555 =  param(364)
      h1a334555 =  param(365)
      h0a135555 =  param(366)
      h1a135555 =  param(367)
      h0a345555 =  param(368)
      h1a345555 =  param(369)
      h0a355555 =  param(370)
      h1a355555 =  param(371)
      !
   endif

   r14    = local(1) ; r24     = local(2) ; r34    = local(3)
   s4a = local(4) ; s4b = local(5) ; rhobar = local(6)

   sinrho = sin(rhobar)
   drho= sin(rhoe)-sinrho

   y1=1.0_ark*(r14-re14) *exp(-beta*(r14-re14)**2)
   y2=1.0_ark*(r24-re14) *exp(-beta*(r24-re14)**2)
   y3=1.0_ark*(r34-re14) *exp(-beta*(r34-re14)**2)

 fea1      =   f0a1      +   f1a1     *drho   +   f2a1      *drho**2  +   f3a1     *drho**3+   f4a1 *drho**4
 fea4      =   f0a4      +   f1a4     *drho   +   f2a4      *drho**2  +   f3a4     *drho**3+   f4a4 *drho**4
 fea11     =   f0a11     +   f1a11    *drho   +   f2a11     *drho**2  +   f3a11    *drho**3
 fea12     =   f0a12     +   f1a12    *drho   +   f2a12     *drho**2  +   f3a12    *drho**3
 fea14     =   f0a14     +   f1a14    *drho   +   f2a14     *drho**2  +   f3a14    *drho**3
 fea24     =   f0a24     +   f1a24    *drho   +   f2a24     *drho**2  +   f3a24    *drho**3
 fea44     =   f0a44     +   f1a44    *drho   +   f2a44     *drho**2  +   f3a44    *drho**3
 fea111    =   f0a111    +   f1a111   *drho   +   f2a111    *drho**2  +   f3a111   *drho**3
 fea112    =   f0a112    +   f1a112   *drho   +   f2a112    *drho**2  +   f3a112   *drho**3
 fea122    =   f0a122    +   f1a122   *drho   +   f2a122    *drho**2  +   f3a122   *drho**3
 fea114    =   f0a114    +   f1a114   *drho   +   f2a114    *drho**2  +   f3a114   *drho**3
 fea124    =   f0a124    +   f1a124   *drho   +   f2a124    *drho**2  +   f3a124   *drho**3
 fea224    =   f0a224    +   f1a224   *drho   +   f2a224    *drho**2  +   f3a224   *drho**3
 fea144    =   f0a144    +   f1a144   *drho   +   f2a144    *drho**2  +   f3a144   *drho**3
 fea244    =   f0a244    +   f1a244   *drho   +   f2a244    *drho**2  +   f3a244   *drho**3
 fea444    =   f0a444    +   f1a444   *drho   +   f2a444    *drho**2  +   f3a444   *drho**3
 fea135    =   f0a135    +   f1a135   *drho   +   f2a135    *drho**2  +   f3a135   *drho**3
 fea155    =   f0a155    +   f1a155   *drho   +   f2a155    *drho**2  +   f3a155   *drho**3
 fea1111   =   f0a1111   +   f1a1111  *drho   +   f2a1111   *drho**2
 fea1112   =   f0a1112   +   f1a1112  *drho   +   f2a1112   *drho**2
 fea1122   =   f0a1122   +   f1a1122  *drho   +   f2a1122   *drho**2
 fea1222   =   f0a1222   +   f1a1222  *drho   +   f2a1222   *drho**2
 fea1123   =   f0a1123   +   f1a1123  *drho   +   f2a1123   *drho**2
 fea1114   =   f0a1114   +   f1a1114  *drho   +   f2a1114   *drho**2
 fea1124   =   f0a1124   +   f1a1124  *drho   +   f2a1124   *drho**2
 fea1224   =   f0a1224   +   f1a1224  *drho   +   f2a1224   *drho**2
 fea2224   =   f0a2224   +   f1a2224  *drho   +   f2a2224   *drho**2
 fea1234   =   f0a1234   +   f1a1234  *drho   +   f2a1234   *drho**2
 fea1144   =   f0a1144   +   f1a1144  *drho   +   f2a1144   *drho**2
 fea1244   =   f0a1244   +   f1a1244  *drho   +   f2a1244   *drho**2
 fea1444   =   f0a1444   +   f1a1444  *drho   +   f2a1444   *drho**2
 fea2444   =   f0a2444   +   f1a2444  *drho   +   f2a2444   *drho**2
 fea4444   =   f0a4444   +   f1a4444  *drho   +   f2a4444   *drho**2
 fea1125   =   f0a1125   +   f1a1125  *drho   +   f2a1125   *drho**2
 fea1225   =   f0a1225   +   f1a1225  *drho   +   f2a1225   *drho**2
 fea1245   =   f0a1245   +   f1a1245  *drho   +   f2a1245   *drho**2
 fea2445   =   f0a2445   +   f1a2445  *drho   +   f2a2445   *drho**2
 fea1155   =   f0a1155   +   f1a1155  *drho   +   f2a1155   *drho**2
 fea1255   =   f0a1255   +   f1a1255  *drho   +   f2a1255   *drho**2
 fea3355   =   f0a3355   +   f1a3355  *drho   +   f2a3355   *drho**2
 fea1455   =   f0a1455   +   f1a1455  *drho   +   f2a1455   *drho**2
 fea4455   =   f0a4455   +   f1a4455  *drho   +   f2a4455   *drho**2
 !
 if (parmax>149) then
   !
   fea11111  =   f0a11111  +   f1a11111 *drho
   fea11112  =   f0a11112  +   f1a11112 *drho
   fea11122  =   f0a11122  +   f1a11122 *drho
   fea11123  =   f0a11123  +   f1a11123 *drho
   fea11223  =   f0a11223  +   f1a11223 *drho
   fea11114  =   f0a11114  +   f1a11114 *drho
   fea11124  =   f0a11124  +   f1a11124 *drho
   fea11224  =   f0a11224  +   f1a11224 *drho
   fea11234  =   f0a11234  +   f1a11234 *drho
   fea22344  =   f0a22344  +   f1a22344 *drho
   fea11135  =   f0a11135  +   f1a11135 *drho
   fea22235  =   f0a22235  +   f1a22235 *drho
   fea11245  =   f0a11245  +   f1a11245 *drho
   fea12245  =   f0a12245  +   f1a12245 *drho
   fea22245  =   f0a22245  +   f1a22245 *drho
   fea11155  =   f0a11155  +   f1a11155 *drho
   fea11255  =   f0a11255  +   f1a11255 *drho
   fea12255  =   f0a12255  +   f1a12255 *drho
   fea22255  =   f0a22255  +   f1a22255 *drho
   fea12355  =   f0a12355  +   f1a12355 *drho
   fea11455  =   f0a11455  +   f1a11455 *drho
   fea12455  =   f0a12455  +   f1a12455 *drho
   fea22455  =   f0a22455  +   f1a22455 *drho
   fea23455  =   f0a23455  +   f1a23455 *drho
   fea24455  =   f0a24455  +   f1a24455 *drho
   fea12555  =   f0a12555  +   f1a12555 *drho
   fea22555  =   f0a22555  +   f1a22555 *drho
   fea24555  =   f0a24555  +   f1a24555 *drho
   fea15555  =   f0a15555  +   f1a15555 *drho
   fea45555  =   f0a45555  +   f1a45555 *drho
   fea111111 =   f0a111111 +   f1a111111*drho
   fea111112 =   f0a111112 +   f1a111112*drho
   fea111122 =   f0a111122 +   f1a111122*drho
   fea122333 =   f0a122333 +   f1a122333*drho
   fea123333 =   f0a123333 +   f1a123333*drho
   fea122334 =   f0a122334 +   f1a122334*drho
   fea123334 =   f0a123334 +   f1a123334*drho
   fea223334 =   f0a223334 +   f1a223334*drho
   fea133334 =   f0a133334 +   f1a133334*drho
   fea233334 =   f0a233334 +   f1a233334*drho
   fea333334 =   f0a333334 +   f1a333334*drho
   fea111344 =   f0a111344 +   f1a111344*drho
   fea112344 =   f0a112344 +   f1a112344*drho
   fea223344 =   f0a223344 +   f1a223344*drho
   fea133344 =   f0a133344 +   f1a133344*drho
   fea233344 =   f0a233344 +   f1a233344*drho
   fea333344 =   f0a333344 +   f1a333344*drho
   fea122444 =   f0a122444 +   f1a122444*drho
   fea222444 =   f0a222444 +   f1a222444*drho
   fea113444 =   f0a113444 +   f1a113444*drho
   fea123444 =   f0a123444 +   f1a123444*drho
   fea124444 =   f0a124444 +   f1a124444*drho
   fea444444 =   f0a444444 +   f1a444444*drho
   fea233335 =   f0a233335 +   f1a233335*drho
   fea233345 =   f0a233345 +   f1a233345*drho
   fea113445 =   f0a113445 +   f1a113445*drho
   fea133445 =   f0a133445 +   f1a133445*drho
   fea233445 =   f0a233445 +   f1a233445*drho
   fea235555 =   f0a235555 +   f1a235555*drho
   fea335555 =   f0a335555 +   f1a335555*drho
   fea145555 =   f0a145555 +   f1a145555*drho
   hea11112  =   h0a11112  +   h1a11112 *drho
   hea11122  =   h0a11122  +   h1a11122 *drho
   hea11244  =   h0a11244  +   h1a11244 *drho
   hea13444  =   h0a13444  +   h1a13444 *drho
   hea33444  =   h0a33444  +   h1a33444 *drho
   hea12335  =   h0a12335  +   h1a12335 *drho
   hea22335  =   h0a22335  +   h1a22335 *drho
   hea13335  =   h0a13335  +   h1a13335 *drho
   hea33335  =   h0a33335  +   h1a33335 *drho
   hea14445  =   h0a14445  +   h1a14445 *drho
   hea34555  =   h0a34555  +   h1a34555 *drho
   hea55555  =   h0a55555  +   h1a55555 *drho
   hea111333 =   h0a111333 +   h1a111333*drho
   hea122333 =   h0a122333 +   h1a122333*drho
   hea223333 =   h0a223333 +   h1a223333*drho
   hea233333 =   h0a233333 +   h1a233333*drho
   hea111334 =   h0a111334 +   h1a111334*drho
   hea113334 =   h0a113334 +   h1a113334*drho
   hea223334 =   h0a223334 +   h1a223334*drho
   hea233334 =   h0a233334 +   h1a233334*drho
   hea113344 =   h0a113344 +   h1a113344*drho
   hea123344 =   h0a123344 +   h1a123344*drho
   hea133344 =   h0a133344 +   h1a133344*drho
   hea233344 =   h0a233344 +   h1a233344*drho
   hea133444 =   h0a133444 +   h1a133444*drho
   hea233444 =   h0a233444 +   h1a233444*drho
   hea333444 =   h0a333444 +   h1a333444*drho
   hea122335 =   h0a122335 +   h1a122335*drho
   hea123335 =   h0a123335 +   h1a123335*drho
   hea333335 =   h0a333335 +   h1a333335*drho
   hea111145 =   h0a111145 +   h1a111145*drho
   hea123345 =   h0a123345 +   h1a123345*drho
   hea223345 =   h0a223345 +   h1a223345*drho
   hea234445 =   h0a234445 +   h1a234445*drho
   hea334445 =   h0a334445 +   h1a334445*drho
   hea144445 =   h0a144445 +   h1a144445*drho
   hea344445 =   h0a344445 +   h1a344445*drho
   hea444445 =   h0a444445 +   h1a444445*drho
   hea333355 =   h0a333355 +   h1a333355*drho
   hea333455 =   h0a333455 +   h1a333455*drho
   hea334455 =   h0a334455 +   h1a334455*drho
   hea344455 =   h0a344455 +   h1a344455*drho
   hea222555 =   h0a222555 +   h1a222555*drho
   hea223555 =   h0a223555 +   h1a223555*drho
   hea114555 =   h0a114555 +   h1a114555*drho
   hea234555 =   h0a234555 +   h1a234555*drho
   hea334555 =   h0a334555 +   h1a334555*drho
   hea135555 =   h0a135555 +   h1a135555*drho
   hea345555 =   h0a345555 +   h1a345555*drho
   hea355555 =   h0a355555 +   h1a355555*drho
   !
 endif
 !
 select case (ix)

 case default
    write(out, '(/a,1x,i3)') 'xy3_dms_E error: unknown value of ix =', ix
    stop 'STOP'

 case (1)

    s2 = fea122*y1*y3**2+fea111*y1**3+(-fea155/2.0_ark-fea144/2.0_ark-fea244)*y3*s4b**2+fea14*y1*s4a+&!
          fea444*s4a*s4b**2+fea1122*y1**2*y3**2+fea2444*y2*s4a**3+fea1114*y1**3*s4a+&!
          fea4455*s4a**2*s4b**2-fea111*y2**3/2.0_ark+fea2444*y3*s4a**3+fea244*y2*s4a**2+&!
          fea2224*y3**3*s4a+(-fea155/2.0_ark-fea144/2.0_ark-fea244)*y2*s4b**2+fea3355*y2**2*s4b**2+&!
          (-fea3355-fea1155/2.0_ark-fea1144/2.0_ark)*y2**2*s4a**2+fea24*y3*s4a+fea1112*y1**3*y2+&!
          fea1122*y1**2*y2**2+fea1144*y1**2*s4a**2+fea1155*y1**2*s4b**2+(-fea1222-fea1112)*y2**3*y3+&!
          fea224*y3**2*s4a+fea12*y1*y3+fea12*y1*y2+(5.0_ark/18.0_ark*fea1444*sqrt(3.0_ark)+&!
          fea1455*sqrt(3.0_ark)/6.0_ark-4.0_ark/9.0_ark*fea2444*sqrt(3.0_ark)+fea2445/3.0_ark)*y3*s4b**3+fea1444*y1*s4a**3

     s3 = fea114*y1**2*s4a+s2+sqrt(3.0_ark)*(fea224-fea114)*y2**2*s4b/3.0_ark+&!
          sqrt(3.0_ark)*(3.0_ark*fea1125+fea1224*sqrt(3.0_ark)+3.0_ark*fea1225+&!
          fea1124*sqrt(3.0_ark))*y2**2*y3*s4a/6.0_ark-fea1111*y3**4/2.0_ark-fea1123*y1*y2**2*y3/2.0_ark-&!
          fea1123*y1*y2*y3**2/2.0_ark-fea1225*y1*y3**2*s4b-fea135*y1*y2*s4b-fea2445*y3*s4a**2*s4b+&!
          (-fea4444-fea4455/3.0_ark)*s4b**4-fea1245*y1*y3*s4a*s4b+&!
          fea124*y1*y3*s4a-sqrt(3.0_ark)*(-fea1114+fea2224)*y3**3*s4b/3.0_ark
     s1 = s3+fea135*y1*y3*s4b+fea124*y1*y2*s4a+sqrt(3.0_ark)*(-fea1114+fea2224)*y2**3*s4b/3.0_ark+&!
          (-fea1455/2.0_ark+fea1444/2.0_ark+fea2444)*y3*s4a*s4b**2+fea4444*s4a**4-fea111*y3**3/2.0_ark-fea1111*y2**4/2.0_ark-&!
          sqrt(3.0_ark)*(4.0_ark*fea3355-fea1155+3.0_ark*fea1144)*y3**2*s4a*s4b/6.0_ark+&!
          (-fea1455/2.0_ark+fea1444/2.0_ark+fea2444)*y2*s4a*s4b**2+(-fea1255/2.0_ark-3.0_ark/2.0_ark*fea1244+&!
          fea1245*sqrt(3.0_ark)/2.0_ark)*y2*y3*s4b**2+fea2445*y2*s4a**2*s4b+&!
          sqrt(3.0_ark)*(fea24-fea14)*y2*s4b/3.0_ark-sqrt(3.0_ark)*(fea24-fea14)*y3*s4b/3.0_ark+&!
          (-fea1224*sqrt(3.0_ark)/2.0_ark-fea1125/2.0_ark+fea1225/2.0_ark+fea1124*sqrt(3.0_ark)/2.0_ark)*y2*y3**2*s4b

     s3 = s1+(-fea1244/2.0_ark-3.0_ark/2.0_ark*fea1255-fea1245*sqrt(3.0_ark)/2.0_ark)*y2*y3*s4a**2-&!
          sqrt(3.0_ark)*(fea224-fea114)*y3**2*s4b/3.0_ark+(fea1224*sqrt(3.0_ark)/2.0_ark+fea1125/2.0_ark-&!
          fea1225/2.0_ark-fea1124*sqrt(3.0_ark)/2.0_ark)*y2**2*y3*s4b+fea3355*y3**2*s4b**2+&!
          (-fea1222-fea1112)*y2*y3**3+fea1255*y1*y3*s4b**2+fea1455*y1*s4a*s4b**2+&!
          fea1244*y1*y3*s4a**2+fea224*y2**2*s4a+(-fea3355-fea1155/2.0_ark-fea1144/2.0_ark)*y3**2*s4a**2+&!
          fea244*y3*s4a**2+sqrt(3.0_ark)*(3.0_ark*fea1125+fea1224*sqrt(3.0_ark)+3.0_ark*fea1225+&!
          fea1124*sqrt(3.0_ark))*y2*y3**2*s4a/6.0_ark

     s2 = s3-fea44*s4b**2+fea1*y1+fea1245*y1*y2*s4a*s4b+fea1222*y1*y2**3+fea1234*y1*y2*y3*s4a+&!
          fea1222*y1*y3**3+sqrt(3.0_ark)*(4.0_ark*fea3355-fea1155+3.0_ark*fea1144)*y2**2*s4a*s4b/6.0_ark-&!
          fea1125*y1**2*y3*s4b+fea1224*y1*y3**2*s4a+fea1244*y1*y2*s4a**2+fea1255*y1*y2*s4b**2-&!
          fea1*y2/2.0_ark-fea1*y3/2.0_ark+fea1111*y1**4
     t4x= s2+fea1123*y1**2*y2*y3+fea1124*y1**2*y2*s4a+fea1224*y1*y2**2*s4a+fea1225*y1*y2**2*s4b+&!
          fea1124*y1**2*y3*s4a+fea1125*y1**2*y2*s4b+fea44*s4a**2+fea24*y2*s4a-fea11*y2**2/2.0_ark-&!
          fea11*y3**2/2.0_ark+fea4*s4a+sqrt(3.0_ark)*(3.0_ark*fea155-fea144+4.0_ark*fea244)*y3*s4a*s4b/6.0_ark+&!
          fea444*s4a**3+fea11*y1**2+(-fea112-fea122)*y2*y3**2+(-fea112-fea122)*y2**2*y3+&!
          sqrt(3.0_ark)*(fea124*sqrt(3.0_ark)-3.0_ark*fea135)*y2*y3*s4a/3.0_ark-sqrt(3.0_ark)*(3.0_ark*fea155-&!
          fea144+4.0_ark*fea244)*y2*s4a*s4b/6.0_ark+fea112*y1**2*y2+fea122*y1*y2**2+fea144*y1*s4a**2+&!
          fea155*y1*s4b**2+fea2224*y2**3*s4a-2.0_ark*fea1122*y2**2*y3**2-2.0_ark*fea12*y2*y3+&!
          fea1112*y1**3*y3+fea112*y1**2*y3+(-5.0_ark/18.0_ark*fea1444*sqrt(3.0_ark)-fea1455*sqrt(3.0_ark)/6.0_ark+&!
          4.0_ark/9.0_ark*fea2444*sqrt(3.0_ark)-fea2445/3.0_ark)*y2*s4b**3
          !
     t56x = 0.0_ark
     if (parmax>149) then
       !
       s6 = (-fea11124+fea11135*sqrt(3.0_ark)/3.0_ark+fea22235*sqrt(3.0_ark)/3.0_ark+2.0_ark*hea13335)*y2**3*y3*s4a+&!
            (hea233334/2.0_ark+fea233334*sqrt(3.0_ark)/3.0_ark-fea133334*sqrt(3.0_ark)/3.0_ark-&!
            fea233335/2.0_ark)*y1**4*y2*s4b-sqrt(3.0_ark)*(3.0_ark*fea122444*sqrt(3.0_ark)+2.0_ark*hea233444+&!
            2.0_ark*fea113444*sqrt(3.0_ark)-6.0_ark*sqrt(3.0_ark)*hea223555+2.0_ark*hea133444+2.0_ark*fea133445+&!
            2.0_ark*fea113445)*y1*y3**2*s4a*s4b**2/3.0_ark-sqrt(3.0_ark)*(21.0_ark*hea234555-9.0_ark*hea234445+&!
            72.0_ark*fea124444-8.0_ark*sqrt(3.0_ark)*hea135555+48.0_ark*fea235555)*y1*y3*s4a**3*s4b/72.0_ark

       s5 = s6+(6.0_ark*hea33444+2.0_ark/3.0_ark*fea22455*sqrt(3.0_ark)-2.0_ark/3.0_ark*fea11455*sqrt(3.0_ark)+&!
            3.0_ark*fea22555)*y3**2*s4a**2*s4b-sqrt(3.0_ark)*(2.0_ark*fea111344*sqrt(3.0_ark)+3.0_ark*hea133344+&!
            3.0_ark*fea233345-3.0_ark*hea233344+fea233344*sqrt(3.0_ark))*y1*y3**3*s4b**2/9.0_ark+&!
            (-47.0_ark/72.0_ark*fea24555*sqrt(3.0_ark)+fea15555/3.0_ark+3.0_ark/8.0_ark*hea14445-25.0_ark/24.0_ark*hea34555+&!
            5.0_ark/6.0_ark*fea24455)*y2*s4a**4-fea22555*y3**2*s4b**3+(9.0_ark/10.0_ark*hea344445-&!
            17.0_ark/30.0_ark*sqrt(3.0_ark)*hea344455+sqrt(3.0_ark)*hea345555/5.0_ark+3.0_ark/2.0_ark*hea355555+&!
            hea144445/5.0_ark-3.0_ark/5.0_ark*fea145555)*y3*s4a**3*s4b**2

       s4 = s5-sqrt(3.0_ark)*(6.0_ark*hea33444-fea22455*sqrt(3.0_ark)-2.0_ark*fea11455*sqrt(3.0_ark))*y2**2*s4a**3/9.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*hea233333+2.0_ark*sqrt(3.0_ark)*fea111112)*y2*y3**5/3.0_ark-sqrt(3.0_ark)*(3.0_ark*hea11122+&!
            sqrt(3.0_ark)*fea11122)*y2**2*y3**3/6.0_ark+sqrt(3.0_ark)*(fea222444*sqrt(3.0_ark)-hea333444+&!
            hea333455)*y3**3*s4a*s4b**2/3.0_ark-sqrt(3.0_ark)*(3.0_ark*fea113444*sqrt(3.0_ark)-4.0_ark*hea233444+&!
            4.0_ark*fea133445-2.0_ark*fea233445+2.0_ark*fea122444*sqrt(3.0_ark)-6.0_ark*sqrt(3.0_ark)*hea223555+&!
            2.0_ark*hea133444)*y1**2*y3*s4a*s4b**2/3.0_ark+(fea45555/3.0_ark+2.0_ark/3.0_ark*hea55555)*s4a**5+&!
            (-6.0_ark*hea33444-2.0_ark/3.0_ark*fea22455*sqrt(3.0_ark)+2.0_ark/3.0_ark*fea11455*sqrt(3.0_ark)-&!
            3.0_ark*fea22555)*y2**2*s4a**2*s4b+(-fea11124+fea11135*sqrt(3.0_ark)/3.0_ark+&!
            fea22235*sqrt(3.0_ark)/3.0_ark+2.0_ark*hea13335)*y2*y3**3*s4a+(-hea233334/2.0_ark-&!
            fea233334*sqrt(3.0_ark)/3.0_ark+fea133334*sqrt(3.0_ark)/3.0_ark+fea233335/2.0_ark)*y1**4*y3*s4b

       s5 = s4+(-193.0_ark/300.0_ark*hea344445-167.0_ark/300.0_ark*sqrt(3.0_ark)*hea344455-&!
            9.0_ark/50.0_ark*sqrt(3.0_ark)*hea345555+63.0_ark/20.0_ark*hea355555-127.0_ark/150.0_ark*hea144445-&!
            33.0_ark/50.0_ark*fea145555)*y3*s4a**5+(27.0_ark/4.0_ark*hea355555-17.0_ark/10.0_ark*fea145555-&!
            39.0_ark/20.0_ark*hea344445-21.0_ark/20.0_ark*sqrt(3.0_ark)*hea344455-sqrt(3.0_ark)*hea345555/10.0_ark-&!
            21.0_ark/10.0_ark*hea144445)*y3*s4a*s4b**4+(2.0_ark/3.0_ark*fea22235*sqrt(3.0_ark)+&!
            2.0_ark/3.0_ark*fea11135*sqrt(3.0_ark)+hea13335)*y1*y2**3*s4a-sqrt(3.0_ark)*(-fea11114+hea33335)*y3**4*s4b+&!
            (3.0_ark/2.0_ark*hea334445+2.0_ark*hea114555-hea334555/2.0_ark+4.0_ark/3.0_ark*sqrt(3.0_ark)*hea334455-&!
            6.0_ark*fea335555)*y3**2*s4a**2*s4b**2+(2.0_ark*fea12255+2.0_ark*fea11255+3.0_ark*fea22344+&!
            fea11245*sqrt(3.0_ark)+fea12245*sqrt(3.0_ark))*y2**2*y3*s4b**2+sqrt(3.0_ark)*(-fea11114+&!
            hea33335)*y2**4*s4b+(-fea133334*sqrt(3.0_ark)/3.0_ark-hea233334+&!
            fea233334*sqrt(3.0_ark)/3.0_ark)*y1*y2**4*s4b-sqrt(3.0_ark)*(hea114555-&!
            hea334555)*y3**2*s4a*s4b**3/3.0_ark

       s3 = s5+(-sqrt(3.0_ark)*hea135555/3.0_ark+hea234555/8.0_ark+3.0_ark/8.0_ark*hea234445)*y1*y2*s4b**4+&!
            sqrt(3.0_ark)*(fea233445+hea233444+2.0_ark*fea122444*sqrt(3.0_ark)+2.0_ark*fea113444*sqrt(3.0_ark)-&!
            3.0_ark*sqrt(3.0_ark)*hea223555-2.0_ark*hea133444+2.0_ark*fea113445)*y2*y3**2*s4a*s4b**2/3.0_ark+&!
            fea11112*y1**4*y3+fea11114*y1**4*s4a+fea333334*y2**5*s4a-sqrt(3.0_ark)*(fea22255*sqrt(3.0_ark)+&!
            2.0_ark*fea11155*sqrt(3.0_ark)+3.0_ark*fea22245)*y3**3*s4a**2/9.0_ark-sqrt(3.0_ark)*(6.0_ark*hea33444-&!
            fea22455*sqrt(3.0_ark)-2.0_ark*fea11455*sqrt(3.0_ark))*y3**2*s4a**3/9.0_ark+fea111122*y1**4*y3**2+&!
            (-fea222444*sqrt(3.0_ark)-2.0_ark*hea333444+sqrt(3.0_ark)*hea222555-hea333455)*y2**3*s4b**3-&!
            sqrt(3.0_ark)*(3.0_ark*hea11112+sqrt(3.0_ark)*fea11112)*y2**4*y3/6.0_ark

       s6 = -sqrt(3.0_ark)*(3.0_ark*fea12255*sqrt(3.0_ark)+6.0_ark*fea11245+6.0_ark*hea11244+8.0_ark*fea22344*sqrt(3.0_ark)+&!
            6.0_ark*fea11255*sqrt(3.0_ark)+9.0_ark*fea12245)*y1*y3**2*s4a**2/3.0_ark-sqrt(3.0_ark)*(3.0_ark*fea12255*sqrt(3.0_ark)+&!
            6.0_ark*fea11245+6.0_ark*hea11244+8.0_ark*fea22344*sqrt(3.0_ark)+6.0_ark*fea11255*sqrt(3.0_ark)+&!
            9.0_ark*fea12245)*y1*y2**2*s4a**2/3.0_ark+(-5.0_ark/8.0_ark*fea24555*sqrt(3.0_ark)+fea15555+&!
            9.0_ark/8.0_ark*hea14445-9.0_ark/8.0_ark*hea34555+3.0_ark/2.0_ark*fea24455)*y2*s4b**4+&!
            sqrt(3.0_ark)*(fea12455*sqrt(3.0_ark)+2.0_ark*fea23455*sqrt(3.0_ark)+6.0_ark*hea13444)*y1*y3*s4a**3/9.0_ark

       s5 = s6-sqrt(3.0_ark)*(3.0_ark*hea11112+sqrt(3.0_ark)*fea11112)*y2*y3**4/6.0_ark+(-4.0_ark/3.0_ark*fea235555-&!
            3.0_ark*fea124444-7.0_ark/9.0_ark*sqrt(3.0_ark)*hea135555-5.0_ark/24.0_ark*hea234555-&!
            5.0_ark/8.0_ark*hea234445)*y1*y2*s4a**2*s4b**2+sqrt(3.0_ark)*(fea12455*sqrt(3.0_ark)+&!
            2.0_ark*fea23455*sqrt(3.0_ark)+6.0_ark*hea13444)*y1*y2*s4a**3/9.0_ark-sqrt(3.0_ark)*(-3.0_ark*fea11124+&!
            3.0_ark*hea13335+fea11135*sqrt(3.0_ark))*y1*y3**3*s4b/3.0_ark+(2.0_ark*fea12255+2.0_ark*fea11255+&!
            3.0_ark*fea22344+fea11245*sqrt(3.0_ark)+fea12245*sqrt(3.0_ark))*y2*y3**2*s4b**2

       s6 = s5+fea444444*s4a**6-sqrt(3.0_ark)*(3.0_ark*hea233333+2.0_ark*sqrt(3.0_ark)*fea111112)*y2**5*y3/3.0_ark+&!
            sqrt(3.0_ark)*(2.0_ark*hea333355+sqrt(3.0_ark)*hea111145)*y2**4*s4b**2/6.0_ark+(2.0_ark/3.0_ark*fea22235*sqrt(3.0_ark)+&!
            2.0_ark/3.0_ark*fea11135*sqrt(3.0_ark)+hea13335)*y1*y3**3*s4a

       s4 = s6+sqrt(3.0_ark)*(-9.0_ark*hea14445-24.0_ark*fea15555+21.0_ark*hea34555+11.0_ark*fea24555*sqrt(3.0_ark)-&!
            24.0_ark*fea24455)*y2*s4a**3*s4b/18.0_ark+(-193.0_ark/300.0_ark*hea344445-&!
            167.0_ark/300.0_ark*sqrt(3.0_ark)*hea344455-9.0_ark/50.0_ark*sqrt(3.0_ark)*hea345555+63.0_ark/20.0_ark*hea355555-&!
            127.0_ark/150.0_ark*hea144445-33.0_ark/50.0_ark*fea145555)*y2*s4a**5+sqrt(3.0_ark)*(3.0_ark*hea233334+&!
            2.0_ark*fea133334*sqrt(3.0_ark)+3.0_ark*fea233335)*y1**4*y2*s4a/6.0_ark-&!
            sqrt(3.0_ark)*(-9.0_ark*hea14445-24.0_ark*fea15555+21.0_ark*hea34555+&!
            11.0_ark*fea24555*sqrt(3.0_ark)-24.0_ark*fea24455)*y3*s4a**3*s4b/18.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*sqrt(3.0_ark)*hea344455+3.0_ark*hea344445+6.0_ark*hea144445+&!
            2.0_ark*sqrt(3.0_ark)*hea345555-15.0_ark*hea355555+6.0_ark*fea145555)*y2*s4a**2*s4b**3/6.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*hea122333+sqrt(3.0_ark)*fea122333)*y1**3*y2*y3**2/6.0_ark

       s5 = s4-sqrt(3.0_ark)*(3.0_ark*hea122333+sqrt(3.0_ark)*fea122333)*y1**3*y2**2*y3/6.0_ark-&!
            sqrt(3.0_ark)*(34.0_ark*hea144445-7.0_ark*hea344445+27.0_ark*sqrt(3.0_ark)*hea344455-&!
            6.0_ark*sqrt(3.0_ark)*hea345555-45.0_ark*hea355555+18.0_ark*fea145555)*y2*s4a**4*s4b/60.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*hea122333-sqrt(3.0_ark)*fea122333)*y1**2*y2*y3**3/6.0_ark-&!
            sqrt(3.0_ark)*(-fea122334+hea122335)*y1**2*y2**2*y3*s4b/4.0_ark+&!
            sqrt(3.0_ark)*(-fea122334+hea122335)*y1**2*y2*y3**2*s4b/4.0_ark-&!
            sqrt(3.0_ark)*(4.0_ark*hea123344-fea112344*sqrt(3.0_ark)-&!
            2.0_ark*sqrt(3.0_ark)*hea123345)*y1**2*y2*y3*s4b**2/9.0_ark+&!
            2.0_ark/3.0_ark*sqrt(3.0_ark)*hea111333*y2**3*y3**3+&!
            sqrt(3.0_ark)*(3.0_ark*hea11112-sqrt(3.0_ark)*fea11112)*y1*y2**4/6.0_ark+&!
            sqrt(3.0_ark)*(hea114555-hea334555)*y2**2*s4a*s4b**3/3.0_ark
       s2 = s5+sqrt(3.0_ark)*(3.0_ark*hea11122-sqrt(3.0_ark)*fea11122)*y1**2*y3**3/6.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*hea11122-sqrt(3.0_ark)*fea11122)*y1**2*y2**3/6.0_ark-&!
            sqrt(3.0_ark)*(4.0_ark*fea22255*sqrt(3.0_ark)-6.0_ark*fea22245-fea11155*sqrt(3.0_ark))*y1**3*s4a**2/9.0_ark+&!
            sqrt(3.0_ark)*(-55.0_ark*hea355555+22.0_ark*fea145555+27.0_ark*hea344445+&!
            13.0_ark*sqrt(3.0_ark)*hea344455+6.0_ark*sqrt(3.0_ark)*hea345555+6.0_ark*hea144445)*y3*s4b**5/100.0_ark-&!
            sqrt(3.0_ark)*(fea22255*sqrt(3.0_ark)+2.0_ark*fea11155*sqrt(3.0_ark)+3.0_ark*fea22245)*y2**3*s4a**2/9.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*hea11112-sqrt(3.0_ark)*fea11112)*y1*y3**4/6.0_ark-sqrt(3.0_ark)*(-55.0_ark*hea355555+&!
            22.0_ark*fea145555+27.0_ark*hea344445+13.0_ark*sqrt(3.0_ark)*hea344455+6.0_ark*sqrt(3.0_ark)*hea345555+&!
            6.0_ark*hea144445)*y2*s4b**5/100.0_ark+sqrt(3.0_ark)*(hea333335-fea333334)*y3**5*s4b/2.0_ark+&!
            fea12455*y1*y2*s4a*s4b**2+sqrt(3.0_ark)*(12.0_ark*hea33444+4.0_ark*fea22455*sqrt(3.0_ark)-&!
            fea11455*sqrt(3.0_ark))*y1**2*s4a**3/9.0_ark+s3

       s6 = -sqrt(3.0_ark)*(-9.0_ark*hea334445-5.0_ark*hea114555+2.0_ark*hea334555+24.0_ark*fea335555-&!
            4.0_ark*sqrt(3.0_ark)*hea334455)*y3**2*s4a**3*s4b/9.0_ark+sqrt(3.0_ark)*(3.0_ark*fea11255*sqrt(3.0_ark)+&!
            4.0_ark*fea22344*sqrt(3.0_ark)+3.0_ark*fea11245+6.0_ark*fea12245+6.0_ark*hea11244)*y1**2*y3*s4a**2/3.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*fea11255*sqrt(3.0_ark)+4.0_ark*fea22344*sqrt(3.0_ark)+3.0_ark*fea11245+6.0_ark*fea12245+&!
            6.0_ark*hea11244)*y1**2*y2*s4a**2/3.0_ark+(-hea233344+fea133344*sqrt(3.0_ark)/3.0_ark+hea133344-&!
            fea233344*sqrt(3.0_ark)/3.0_ark)*y1**3*y2*s4a*s4b

       s5 = s6-sqrt(3.0_ark)*(hea333455-2.0_ark*fea222444*sqrt(3.0_ark)-7.0_ark*hea333444)*y1**3*s4a*s4b**2/6.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*hea122333-sqrt(3.0_ark)*fea122333)*y1**2*y2**3*y3/6.0_ark+sqrt(3.0_ark)*(hea111334+&!
            2.0_ark*hea113334+hea223334+fea223334*sqrt(3.0_ark))*y1**3*y3**2*s4a/3.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*fea113444*sqrt(3.0_ark)-4.0_ark*hea233444+4.0_ark*fea133445-2.0_ark*fea233445+&!
            2.0_ark*fea122444*sqrt(3.0_ark)-6.0_ark*sqrt(3.0_ark)*hea223555+2.0_ark*hea133444)*y1**2*y2*s4a*s4b**2/3.0_ark+&!
            sqrt(3.0_ark)*(2.0_ark*hea333355+sqrt(3.0_ark)*hea111145)*y3**4*s4b**2/6.0_ark

       s4 = s5+(-sqrt(3.0_ark)*hea111145/2.0_ark-fea333344*sqrt(3.0_ark)+hea333355)*y3**4*s4a*s4b-&!
            sqrt(3.0_ark)*(3.0_ark*hea11122+sqrt(3.0_ark)*fea11122)*y2**3*y3**2/6.0_ark+fea111111*y1**6-&!
            sqrt(3.0_ark)*(3.0_ark*fea122444*sqrt(3.0_ark)+2.0_ark*hea233444+2.0_ark*fea113444*sqrt(3.0_ark)-&!
            6.0_ark*sqrt(3.0_ark)*hea223555+2.0_ark*hea133444+2.0_ark*fea133445+2.0_ark*fea113445)*y1*y2**2*s4a*s4b**2/3.0_ark+&!
            sqrt(3.0_ark)*(hea123335-fea123334)*y1*y2*y3**3*s4b/2.0_ark+(-4.0_ark/3.0_ark*hea123344-&!
            sqrt(3.0_ark)*hea123345/3.0_ark-2.0_ark/3.0_ark*fea112344*sqrt(3.0_ark))*y1*y2*y3**2*s4a*s4b+&!
            sqrt(3.0_ark)*(34.0_ark*hea144445-7.0_ark*hea344445+27.0_ark*sqrt(3.0_ark)*hea344455-&!
            6.0_ark*sqrt(3.0_ark)*hea345555-45.0_ark*hea355555+18.0_ark*fea145555)*y3*s4a**4*s4b/60.0_ark+&!
            (-2.0_ark/3.0_ark*hea444445-7.0_ark/3.0_ark*fea444444)*s4a**2*s4b**4+(hea233344-fea133344*sqrt(3.0_ark)/3.0_ark-&!
            hea133344+fea233344*sqrt(3.0_ark)/3.0_ark)*y1**3*y3*s4a*s4b+fea133334*y1*y3**4*s4a

       s5 = s4+fea122444*y1*y2**2*s4a**3-fea24555*y3*s4a*s4b**3+sqrt(3.0_ark)*(-3.0_ark*fea11124+3.0_ark*hea13335+&!
            fea11135*sqrt(3.0_ark))*y1*y2**3*s4b/3.0_ark-sqrt(3.0_ark)*(-4.0_ark*fea12455*sqrt(3.0_ark)+fea23455*sqrt(3.0_ark)+&!
            12.0_ark*hea13444)*y2*y3*s4a**3/9.0_ark+(fea133344*sqrt(3.0_ark)/3.0_ark-2.0_ark/3.0_ark*fea111344*sqrt(3.0_ark)-&!
            hea133344-fea233345-hea233344+fea233344*sqrt(3.0_ark)/3.0_ark)*y1*y3**3*s4a*s4b+&!
            fea24455*y2*s4a**2*s4b**2+fea122333*y1*y2**2*y3**3+fea233344*y2**3*y3*s4a**2+&!
            (-2.0_ark/3.0_ark*hea113334+2.0_ark/3.0_ark*hea111334-hea223334/3.0_ark)*y2**3*y3**2*s4b

       s3 = s5+fea133344*y1*y2**3*s4a**2+fea133334*y1*y2**4*s4a+fea233334*y2**4*y3*s4a+&!
            fea223334*y2**3*y3**2*s4a+fea122333*y1*y2**3*y3**2+fea124444*y1*y3*s4a**4-fea233335*y2**4*y3*s4b+&!
            fea233334*y2*y3**4*s4a+(-sqrt(3.0_ark)*hea223345/6.0_ark+2.0_ark/3.0_ark*hea113344-&!
            fea223344*sqrt(3.0_ark)/3.0_ark)*y1**2*y2**2*s4a*s4b+sqrt(3.0_ark)*(3.0_ark*sqrt(3.0_ark)*hea344455+&!
            3.0_ark*hea344445+6.0_ark*hea144445+2.0_ark*sqrt(3.0_ark)*hea345555-15.0_ark*hea355555+&!
            6.0_ark*fea145555)*y3*s4a**2*s4b**3/6.0_ark

       s5 = (-hea444445-fea444444)*s4a**4*s4b**2+fea15555*y1*s4b**4+fea335555*y2**2*s4b**4+&!
            (-fea133344*sqrt(3.0_ark)/3.0_ark+2.0_ark/3.0_ark*fea111344*sqrt(3.0_ark)+hea133344+fea233345+&!
            hea233344-fea233344*sqrt(3.0_ark)/3.0_ark)*y1*y2**3*s4a*s4b+(4.0_ark/3.0_ark*hea123344+&!
            sqrt(3.0_ark)*hea123345/3.0_ark+2.0_ark/3.0_ark*fea112344*sqrt(3.0_ark))*y1*y2**2*y3*s4a*s4b+&!
            fea145555*y1*s4a*s4b**4+(27.0_ark/4.0_ark*hea355555-17.0_ark/10.0_ark*fea145555-&!
            39.0_ark/20.0_ark*hea344445-21.0_ark/20.0_ark*sqrt(3.0_ark)*hea344455-sqrt(3.0_ark)*hea345555/10.0_ark-&!
            21.0_ark/10.0_ark*hea144445)*y2*s4a*s4b**4+fea113444*y1**2*y3*s4a**3+&!
            (9.0_ark/10.0_ark*hea344445-17.0_ark/30.0_ark*sqrt(3.0_ark)*hea344455+sqrt(3.0_ark)*hea345555/5.0_ark+&!
            3.0_ark/2.0_ark*hea355555+hea144445/5.0_ark-3.0_ark/5.0_ark*fea145555)*y2*s4a**3*s4b**2

       s6 = s5+(3.0_ark*fea222444*sqrt(3.0_ark)+3.0_ark*hea333444-3.0_ark*sqrt(3.0_ark)*hea222555+&!
            2.0_ark*hea333455)*y2**3*s4a**2*s4b+fea122334*y1*y2**2*y3**2*s4a+(-3.0_ark*fea222444*sqrt(3.0_ark)-&!
            3.0_ark*hea333444+3.0_ark*sqrt(3.0_ark)*hea222555-2.0_ark*hea333455)*y3**3*s4a**2*s4b+&!
            (fea133334*sqrt(3.0_ark)/3.0_ark+hea233334-fea233334*sqrt(3.0_ark)/3.0_ark)*y1*y3**4*s4b

       s4 = s6+(-6.0_ark*hea13444+2.0_ark/3.0_ark*fea12455*sqrt(3.0_ark)-2.0_ark/3.0_ark*fea23455*sqrt(3.0_ark)-&!
            3.0_ark*fea12555)*y1*y2*s4a**2*s4b+fea133344*y1*y3**3*s4a**2+(-10.0_ark/3.0_ark*fea235555-&!
            12.0_ark*fea124444-4.0_ark/9.0_ark*sqrt(3.0_ark)*hea135555-10.0_ark/3.0_ark*hea234555-&!
            hea234445)*y2*y3*s4a**2*s4b**2+(-4.0_ark/3.0_ark*fea235555-3.0_ark*fea124444-&!
            7.0_ark/9.0_ark*sqrt(3.0_ark)*hea135555-5.0_ark/24.0_ark*hea234555-5.0_ark/8.0_ark*hea234445)*y1*y3*s4a**2*s4b**2+&!
            (-5.0_ark*fea11245-2.0_ark*fea12255*sqrt(3.0_ark)-4.0_ark*fea11255*sqrt(3.0_ark)-6.0_ark*fea22344*sqrt(3.0_ark)-&!
            7.0_ark*fea12245-6.0_ark*hea11244)*y2*y3**2*s4a*s4b+(6.0_ark*hea13444-2.0_ark/3.0_ark*fea12455*sqrt(3.0_ark)+&!
            2.0_ark/3.0_ark*fea23455*sqrt(3.0_ark)+3.0_ark*fea12555)*y1*y3*s4a**2*s4b

       s5 = s4-sqrt(3.0_ark)*(hea12335-fea11234)*y1*y2*y3**2*s4b+fea113444*y1**2*y2*s4a**3+&!
            sqrt(3.0_ark)*(2.0_ark*fea222444*sqrt(3.0_ark)+3.0_ark*hea333444+3.0_ark*hea333455)*y1**3*s4a**3/6.0_ark-&!
            sqrt(3.0_ark)*hea111333*y1**3*y3**3/3.0_ark+(4.0_ark*fea124444+4.0_ark/3.0_ark*sqrt(3.0_ark)*hea135555+hea234555+&!
            fea235555)*y2*y3*s4a**4-sqrt(3.0_ark)*hea111333*y1**3*y2**3/3.0_ark+(-4.0_ark/3.0_ark*fea133445+&!
            4.0_ark/3.0_ark*hea233444+2.0_ark/3.0_ark*fea233445-10.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)-&!
            8.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)+2.0_ark*sqrt(3.0_ark)*hea223555-2.0_ark*hea133444+&!
            fea113445/3.0_ark)*y1**2*y2*s4b**3+sqrt(3.0_ark)*(-sqrt(3.0_ark)*hea111145-3.0_ark*fea333344*sqrt(3.0_ark)+&!
            hea333355)*y1**4*s4b**2/3.0_ark+(5.0_ark*fea11245+2.0_ark*fea12255*sqrt(3.0_ark)+4.0_ark*fea11255*sqrt(3.0_ark)+&!
            6.0_ark*fea22344*sqrt(3.0_ark)+7.0_ark*fea12245+6.0_ark*hea11244)*y2**2*y3*s4a*s4b

       s1 = s5+sqrt(3.0_ark)*(21.0_ark*hea234555-9.0_ark*hea234445+72.0_ark*fea124444-8.0_ark*sqrt(3.0_ark)*hea135555+&!
            48.0_ark*fea235555)*y1*y2*s4a**3*s4b/72.0_ark+fea123444*y1*y2*y3*s4a*s4b**2-&!
            sqrt(3.0_ark)*(40.0_ark*sqrt(3.0_ark)*hea135555+75.0_ark*hea234555+9.0_ark*hea234445+216.0_ark*fea124444+&!
            48.0_ark*fea235555)*y1*y2*s4a*s4b**3/72.0_ark+fea22344*y2*y3**2*s4a**2+fea12255*y1*y3**2*s4b**2-&!
            fea11135*y1**3*y2*s4b+fea11223*y1**2*y2*y3**2+fea11224*y1**2*y2**2*s4a-fea111111*y2**6/2.0_ark+&!
            s2+s3

       s4 = s1+fea22235*y2**3*y3*s4b+fea22344*y2**2*y3*s4a**2+(-3.0_ark*hea22335+4.0_ark*fea11224)*y2**2*y3**2*s4a+&!
            fea22245*y2**3*s4a*s4b-fea22235*y2*y3**3*s4b+fea222444*y2**3*s4a**3+sqrt(3.0_ark)*(3.0_ark*hea233333+&!
            sqrt(3.0_ark)*fea111112)*y1*y2**5/3.0_ark+sqrt(3.0_ark)*(3.0_ark*hea233333+&!
            sqrt(3.0_ark)*fea111112)*y1*y3**5/3.0_ark-fea22245*y3**3*s4a*s4b+fea22455*y3**2*s4a*s4b**2+&!
            fea12255*y1*y2**2*s4b**2-fea12555*y1*y3*s4b**3+fea24555*y2*s4a*s4b**3+fea233345*y2*y3**3*s4a*s4b+&!
            fea123444*y1*y2*y3*s4a**3+fea23455*y2*y3*s4a*s4b**2+fea123334*y1*y2*y3**3*s4a-&!
            fea233445*y2**2*y3*s4a**2*s4b

       s3 = s4-fea111111*y3**6/2.0_ark+fea11223*y1**2*y2**2*y3+fea11224*y1**2*y3**2*s4a+&!
            fea11455*y1**2*s4a*s4b**2+fea11255*y1**2*y3*s4b**2-sqrt(3.0_ark)*(hea22335-fea11224)*y1**2*y2**2*s4b+&!
            sqrt(3.0_ark)*(hea22335-fea11224)*y1**2*y3**2*s4b+sqrt(3.0_ark)*(hea12335-fea11234)*y1*y2**2*y3*s4b-&!
            fea133445*y1*y2**2*s4a**2*s4b-fea11111*y2**5/2.0_ark+(sqrt(3.0_ark)*hea223345/6.0_ark-2.0_ark/3.0_ark*hea113344+&!
            fea223344*sqrt(3.0_ark)/3.0_ark)*y1**2*y3**2*s4a*s4b+fea113445*y1**2*y3*s4a**2*s4b-&!
            fea113445*y1**2*y2*s4a**2*s4b+sqrt(3.0_ark)*(5.0_ark*hea123344+2.0_ark*sqrt(3.0_ark)*hea123345+&!
            fea112344*sqrt(3.0_ark))*y1*y2*y3**2*s4b**2/9.0_ark-fea233345*y2**3*y3*s4a*s4b+fea12555*y1*y2*s4b**3+&!
            fea223344*y2**2*y3**2*s4a**2+fea11124*y1**3*y3*s4a+fea123333*y1*y2*y3**4

       s4 = fea124444*y1*y2*s4a**4+fea223334*y2**2*y3**3*s4a+fea11124*y1**3*y2*s4a+fea111122*y1**4*y2**2+&!
            fea133445*y1*y3**2*s4a**2*s4b+(-2.0_ark*fea11114+3.0_ark*hea33335)*y3**4*s4a+&!
            fea112344*y1**2*y2*y3*s4a**2+fea123334*y1*y2**3*y3*s4a+fea123333*y1*y2**4*y3+&!
            fea122444*y1*y3**2*s4a**3+(hea113334/3.0_ark+2.0_ark/3.0_ark*hea223334+&!
            2.0_ark/3.0_ark*hea111334)*y1**2*y3**3*s4b+fea233344*y2*y3**3*s4a**2-fea11111*y3**5/2.0_ark+&!
            (2.0_ark/3.0_ark*hea113334-2.0_ark/3.0_ark*hea111334+hea223334/3.0_ark)*y2**2*y3**3*s4b+fea11123*y1**3*y2*y3+&!
            fea233335*y2*y3**4*s4b+(fea122334/4.0_ark+3.0_ark/4.0_ark*hea122335)*y1**2*y2*y3**2*s4a+&!
            sqrt(3.0_ark)*(hea113334+2.0_ark*hea111334-hea223334+fea223334*sqrt(3.0_ark))*y1**2*y3**3*s4a/3.0_ark+&!
            (fea122334/4.0_ark+3.0_ark/4.0_ark*hea122335)*y1**2*y2**2*y3*s4a

       s5 = s4+fea235555*y2*y3*s4b**4+(-hea113334/3.0_ark-2.0_ark/3.0_ark*hea223334-&!
            2.0_ark/3.0_ark*hea111334)*y1**2*y2**3*s4b-2.0_ark*fea123333*y1**4*y2*y3+(-fea123334/2.0_ark+&!
            3.0_ark/2.0_ark*hea123335)*y1**3*y2*y3*s4a+fea111344*y1**3*y3*s4a**2+fea111112*y1**5*y3+&!
            (-2.0_ark/3.0_ark*hea113334-hea111334/3.0_ark+2.0_ark/3.0_ark*hea223334)*y1**3*y2**2*s4b+&!
            fea111344*y1**3*y2*s4a**2+(2.0_ark/3.0_ark*hea113334+hea111334/3.0_ark-&!
            2.0_ark/3.0_ark*hea223334)*y1**3*y3**2*s4b

       s2 = s5+(-2.0_ark*fea11234+3.0_ark*hea12335)*y1*y2*y3**2*s4a+(-hea334445-11.0_ark/18.0_ark*hea114555-&!
            hea334555/18.0_ark+7.0_ark/3.0_ark*fea335555-13.0_ark/18.0_ark*sqrt(3.0_ark)*hea334455)*y1**2*s4a**4+&!
            fea22455*y2**2*s4a*s4b**2+fea11255*y1**2*y2*s4b**2-2.0_ark*fea11223*y1*y2**2*y3**2-&!
            sqrt(3.0_ark)*(6.0_ark*fea133445-3.0_ark*fea233445-3.0_ark*hea233444+4.0_ark*fea122444*sqrt(3.0_ark)+&!
            4.0_ark*fea113444*sqrt(3.0_ark)-9.0_ark*sqrt(3.0_ark)*hea223555+6.0_ark*hea133444)*y2*y3**2*s4a**3/3.0_ark+&!
            (-32.0_ark/75.0_ark*hea144445-19.0_ark/75.0_ark*hea344445-11.0_ark/75.0_ark*sqrt(3.0_ark)*hea344455+&!
            6.0_ark/25.0_ark*sqrt(3.0_ark)*hea345555+9.0_ark/5.0_ark*hea355555-3.0_ark/25.0_ark*fea145555)*y1*s4a**5+&!
            fea111112*y1**5*y2+fea24455*y3*s4a**2*s4b**2+(10.0_ark/3.0_ark*hea55555-&!
            4.0_ark/3.0_ark*fea45555)*s4a**3*s4b**2+s3

       s5 = -sqrt(3.0_ark)*(3.0_ark*hea223333+2.0_ark*sqrt(3.0_ark)*fea111122)*y2**4*y3**2/3.0_ark+&!
            (hea444445/3.0_ark-fea444444/3.0_ark)*s4b**6+fea11111*y1**5+fea12455*y1*y3*s4a*s4b**2+&!
            sqrt(3.0_ark)*(5.0_ark*hea123344+2.0_ark*sqrt(3.0_ark)*hea123345+fea112344*sqrt(3.0_ark))*y1*y2**2*y3*s4b**2/9.0_ark-&!
            fea12245*y1*y3**2*s4a*s4b+fea11135*y1**3*y3*s4b+fea12355*y1*y2*y3*s4b**2+(-hea114555/2.0_ark+&!
            hea334555/2.0_ark+fea335555-sqrt(3.0_ark)*hea334455/2.0_ark)*y1**2*s4b**4

       s4 = s5-fea12355*y1*y2*y3*s4a**2+fea12245*y1*y2**2*s4a*s4b+sqrt(3.0_ark)*(-2.0_ark*fea133344*sqrt(3.0_ark)+&!
            6.0_ark*hea133344+3.0_ark*fea233345-fea233344*sqrt(3.0_ark))*y2*y3**3*s4b**2/9.0_ark+s2-&!
            fea11245*y1**2*y3*s4a*s4b+fea11234*y1**2*y2*y3*s4a-sqrt(3.0_ark)*(6.0_ark*fea133445-&!
            3.0_ark*fea233445-3.0_ark*hea233444+4.0_ark*fea122444*sqrt(3.0_ark)+4.0_ark*fea113444*sqrt(3.0_ark)-&!
            9.0_ark*sqrt(3.0_ark)*hea223555+6.0_ark*hea133444)*y2**2*y3*s4a**3/3.0_ark-sqrt(3.0_ark)*(hea123344+&!
            sqrt(3.0_ark)*hea123345+fea112344*sqrt(3.0_ark))*y1*y2**2*y3*s4a**2/3.0_ark-sqrt(3.0_ark)*(3.0_ark*hea223333+&!
            2.0_ark*sqrt(3.0_ark)*fea111122)*y2**2*y3**4/3.0_ark-sqrt(3.0_ark)*(sqrt(3.0_ark)*hea223345+&!
            2.0_ark*hea113344)*y1**2*y2**2*s4a**2/6.0_ark

       s5 = s4+(-47.0_ark/72.0_ark*fea24555*sqrt(3.0_ark)+fea15555/3.0_ark+3.0_ark/8.0_ark*hea14445-25.0_ark/24.0_ark*hea34555+&!
            5.0_ark/6.0_ark*fea24455)*y3*s4a**4-sqrt(3.0_ark)*(hea123335-fea123334)*y1*y2**3*y3*s4b/2.0_ark+&!
            sqrt(3.0_ark)*(fea222444*sqrt(3.0_ark)-hea333444+hea333455)*y2**3*s4a*s4b**2/3.0_ark+&!
            fea333334*y3**5*s4a+fea11245*y1**2*y2*s4a*s4b+(4.0_ark/3.0_ark*fea133445-4.0_ark/3.0_ark*hea233444-&!
            2.0_ark/3.0_ark*fea233445+10.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)+8.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)-&!
            2.0_ark*sqrt(3.0_ark)*hea223555+2.0_ark*hea133444-fea113445/3.0_ark)*y1**2*y3*s4b**3+(-6.0_ark*fea15555-&!
            8.0_ark*fea24455+7.0_ark/2.0_ark*fea24555*sqrt(3.0_ark)-9.0_ark/2.0_ark*hea14445+&!
            15.0_ark/2.0_ark*hea34555)*y1*s4a**2*s4b**2+fea233445*y2*y3**2*s4a**2*s4b+fea222444*y3**3*s4a**3

       s3 = s5+(-8.0_ark/3.0_ark*fea24455+25.0_ark/18.0_ark*fea24555*sqrt(3.0_ark)-5.0_ark/3.0_ark*fea15555-&!
            3.0_ark/2.0_ark*hea14445+11.0_ark/6.0_ark*hea34555)*y1*s4a**4+fea22255*y3**3*s4b**2+&!
            fea22555*y2**2*s4b**3+fea11122*y1**3*y3**2-sqrt(3.0_ark)*(-sqrt(3.0_ark)*hea223345+hea113344+&!
            fea223344*sqrt(3.0_ark))*y1**2*y2**2*s4b**2/9.0_ark+(-2.0_ark*fea11234+3.0_ark*hea12335)*y1*y2**2*y3*s4a+&!
            sqrt(3.0_ark)*(hea113334+2.0_ark*hea111334-hea223334+fea223334*sqrt(3.0_ark))*y1**2*y2**3*s4a/3.0_ark+&!
            sqrt(3.0_ark)*(3.0_ark*hea223333+sqrt(3.0_ark)*fea111122)*y1**2*y3**4/3.0_ark+sqrt(3.0_ark)*(3.0_ark*hea223333+&!
            sqrt(3.0_ark)*fea111122)*y1**2*y2**4/3.0_ark+sqrt(3.0_ark)*(-2.0_ark*fea133344*sqrt(3.0_ark)+6.0_ark*hea133344+&!
            3.0_ark*fea233345-fea233344*sqrt(3.0_ark))*y2**3*y3*s4b**2/9.0_ark

       s6 = s3+(hea334445/2.0_ark+2.0_ark/9.0_ark*hea114555-7.0_ark/18.0_ark*hea334555-5.0_ark/3.0_ark*fea335555+&!
            4.0_ark/9.0_ark*sqrt(3.0_ark)*hea334455)*y3**2*s4a**4+(-8.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)+&!
            2.0_ark/3.0_ark*hea233444-10.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)+2.0_ark*sqrt(3.0_ark)*hea223555-&!
            2.0_ark*hea133444-5.0_ark/3.0_ark*fea133445+2.0_ark/3.0_ark*fea113445+4.0_ark/3.0_ark*fea233445)*y1*y2**2*s4b**3-&!
            sqrt(3.0_ark)*(sqrt(3.0_ark)*hea223345+2.0_ark*hea113344)*y1**2*y3**2*s4a**2/6.0_ark-sqrt(3.0_ark)*(-&!
            sqrt(3.0_ark)*hea223345+hea113344+fea223344*sqrt(3.0_ark))*y1**2*y3**2*s4b**2/9.0_ark

       s5 = s6+(-fea233445/3.0_ark+2.0_ark/3.0_ark*fea133445-2.0_ark/3.0_ark*hea233444+2.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)-&!
            2.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)-2.0_ark/3.0_ark*fea113445)*y2**2*y3*s4b**3-fea11123*y1*y2*y3**3/2.0_ark+&!
            (8.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)-2.0_ark/3.0_ark*hea233444+10.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)-&!
            2.0_ark*sqrt(3.0_ark)*hea223555+2.0_ark*hea133444+5.0_ark/3.0_ark*fea133445-2.0_ark/3.0_ark*fea113445-&!
            4.0_ark/3.0_ark*fea233445)*y1*y3**2*s4b**3-fea11123*y1*y2**3*y3/2.0_ark+(fea222444*sqrt(3.0_ark)+&!
            2.0_ark*hea333444-sqrt(3.0_ark)*hea222555+hea333455)*y3**3*s4b**3

       s4 = s5-sqrt(3.0_ark)*(hea333335-fea333334)*y2**5*s4b/2.0_ark-sqrt(3.0_ark)*(hea123344+sqrt(3.0_ark)*hea123345+&!
            fea112344*sqrt(3.0_ark))*y1*y2*y3**2*s4a**2/3.0_ark+fea333344*y2**4*s4a**2+(hea334445/2.0_ark+&!
            2.0_ark/9.0_ark*hea114555-7.0_ark/18.0_ark*hea334555-5.0_ark/3.0_ark*fea335555+&!
            4.0_ark/9.0_ark*sqrt(3.0_ark)*hea334455)*y2**2*s4a**4+fea11155*y1**3*s4b**2+&!
            (3.0_ark/2.0_ark*hea333335-fea333334/2.0_ark)*y1**5*s4a+fea11112*y1**4*y2-&!
            sqrt(3.0_ark)*(2.0_ark*fea111344*sqrt(3.0_ark)+3.0_ark*hea133344+3.0_ark*fea233345-3.0_ark*hea233344+&!
            fea233344*sqrt(3.0_ark))*y1*y2**3*s4b**2/9.0_ark+(-5.0_ark/8.0_ark*fea24555*sqrt(3.0_ark)+&!
            fea15555+9.0_ark/8.0_ark*hea14445-9.0_ark/8.0_ark*hea34555+3.0_ark/2.0_ark*fea24455)*y3*s4b**4+&!
            fea45555*s4a*s4b**4

       s5 = s4+fea22255*y2**3*s4b**2+(-2.0_ark*fea11114+3.0_ark*hea33335)*y2**4*s4a+&!
            fea335555*y3**2*s4b**4+fea333344*y3**4*s4a**2+sqrt(3.0_ark)*(fea233445+&!
            hea233444+2.0_ark*fea122444*sqrt(3.0_ark)+2.0_ark*fea113444*sqrt(3.0_ark)-3.0_ark*sqrt(3.0_ark)*hea223555-&!
            2.0_ark*hea133444+2.0_ark*fea113445)*y2**2*y3*s4a*s4b**2/3.0_ark+sqrt(3.0_ark)*(8.0_ark*hea113344+&!
            sqrt(3.0_ark)*hea223345-fea223344*sqrt(3.0_ark))*y2**2*y3**2*s4b**2/9.0_ark+sqrt(3.0_ark)*(-&!
            9.0_ark*hea334445-5.0_ark*hea114555+2.0_ark*hea334555+24.0_ark*fea335555-&!
            4.0_ark*sqrt(3.0_ark)*hea334455)*y2**2*s4a**3*s4b/9.0_ark+fea11122*y1**3*y2**2-&!
            sqrt(3.0_ark)*(3.0_ark*hea233344+fea133344*sqrt(3.0_ark)+3.0_ark*hea133344+fea233344*sqrt(3.0_ark)+&!
            fea111344*sqrt(3.0_ark))*y1**3*y2*s4b**2/9.0_ark

       s6 = s5+(sqrt(3.0_ark)*hea111145/2.0_ark+fea333344*sqrt(3.0_ark)-hea333355)*y2**4*s4a*s4b+(-3.0_ark*hea334445-&!
            2.0_ark*hea114555+2.0_ark*hea334555+6.0_ark*fea335555-5.0_ark/3.0_ark*sqrt(3.0_ark)*hea334455)*y1**2*s4a**2*s4b**2+&!
            (fea233445/3.0_ark-2.0_ark/3.0_ark*fea133445+2.0_ark/3.0_ark*hea233444-2.0_ark/9.0_ark*fea122444*sqrt(3.0_ark)+&!
            2.0_ark/9.0_ark*fea113444*sqrt(3.0_ark)+2.0_ark/3.0_ark*fea113445)*y2*y3**2*s4b**3+(6.0_ark/5.0_ark*fea145555+&!
            8.0_ark/5.0_ark*hea144445+11.0_ark/5.0_ark*hea344445+17.0_ark/15.0_ark*sqrt(3.0_ark)*hea344455-&!
            2.0_ark/5.0_ark*sqrt(3.0_ark)*hea345555-3.0_ark*hea355555)*y1*s4a**3*s4b**2-&!
            sqrt(3.0_ark)*(3.0_ark*hea233344+fea133344*sqrt(3.0_ark)+3.0_ark*hea133344+&!
            fea233344*sqrt(3.0_ark)+fea111344*sqrt(3.0_ark))*y1**3*y3*s4b**2/9.0_ark

       t56x=s6+sqrt(3.0_ark)*(40.0_ark*sqrt(3.0_ark)*hea135555+75.0_ark*hea234555+9.0_ark*hea234445+216.0_ark*fea124444+&!
            48.0_ark*fea235555)*y1*y3*s4a*s4b**3/72.0_ark+sqrt(3.0_ark)*(hea111334+2.0_ark*hea113334+hea223334+&!
            fea223334*sqrt(3.0_ark))*y1**3*y2**2*s4a/3.0_ark+(-sqrt(3.0_ark)*hea135555/3.0_ark+hea234555/8.0_ark+&!
            3.0_ark/8.0_ark*hea234445)*y1*y3*s4b**4-sqrt(3.0_ark)*(3.0_ark*hea333355-&!
            fea333344*sqrt(3.0_ark))*y1**4*s4a**2/3.0_ark+(3.0_ark/2.0_ark*hea334445+2.0_ark*hea114555-&!
            hea334555/2.0_ark+4.0_ark/3.0_ark*sqrt(3.0_ark)*hea334455-6.0_ark*fea335555)*y2**2*s4a**2*s4b**2+&!
            sqrt(3.0_ark)*(3.0_ark*hea233334+2.0_ark*fea133334*sqrt(3.0_ark)+3.0_ark*fea233335)*y1**4*y3*s4a/6.0_ark
     endif
     !
     f =( t4x+t56x )
     !
case (2)

     s3 = sqrt(3.0_ark)*fea1*y2/2.0_ark+ (2.0_ark/3.0_ark*fea1114+ fea2224/3.0_ark)*y2**3*s4b+ (fea24/3.0_ark+ &!
          2.0_ark/3.0_ark*fea14)*y3*s4b+ (-fea1114/3.0_ark+ 4.0_ark/3.0_ark*fea2224)*y1**3*s4b+  (fea224/3.0_ark+&!
           2.0_ark/3.0_ark*fea114)*y2**2*s4b+ (2.0_ark/3.0_ark*fea1114+ fea2224/3.0_ark)*y3**3*s4b+ (-&!
          2.0_ark*fea4444-fea4455)*s4a**3*s4b-sqrt(3.0_ark)*(-fea112+ fea122)*y2**2*y3/3.0_ark-&!
          fea135*y1*y2*s4a+ sqrt(3.0_ark)*(fea122+ 2.0_ark*fea112)*y1*y2**2/3.0_ark-sqrt(3.0_ark)*(fea112+ &!
          2.0_ark*fea122)*y1**2*y3/3.0_ark+ sqrt(3.0_ark)*(-fea112+ fea122)*y2*y3**2/3.0_ark

     s2 = s3-sqrt(3.0_ark)*fea111*y3**3/2.0_ark-sqrt(3.0_ark)*fea1111*y3**4/2.0_ark+ &!
          sqrt(3.0_ark)*fea111*y2**3/2.0_ark+ sqrt(3.0_ark)*(fea3355-fea1155)*y3**2*s4a**2/3.0_ark-sqrt(3.0_ark)*(-&!
          fea1114+ fea2224)*y3**3*s4a/3.0_ark+ sqrt(3.0_ark)*(16.0_ark*fea2444*sqrt(3.0_ark)+ &!
          17.0_ark*fea1444*sqrt(3.0_ark)+ 3.0_ark*fea1455*sqrt(3.0_ark)-12.0_ark*fea2445)*y2*s4b**3/108.0_ark+ &!
          sqrt(3.0_ark)*(-fea1444-3.0_ark*fea1455+ 4.0_ark*fea2444)*y2*s4a**3/12.0_ark+ &!
          sqrt(3.0_ark)*(2.0_ark*fea3355+ fea1155+ 3.0_ark*fea1144)*y2**2*s4b**2/6.0_ark+ fea135*y1*y3*s4a-&!
          sqrt(3.0_ark)*(fea3355-fea1155)*y2**2*s4a**2/3.0_ark+ sqrt(3.0_ark)*(-fea1114+ &!
          fea2224)*y2**3*s4a/3.0_ark+ sqrt(3.0_ark)*(fea224-fea114)*y2**2*s4a/3.0_ark

     s3 = -sqrt(3.0_ark)*(fea1222-fea1112)*y2**3*y3/3.0_ark+ sqrt(3.0_ark)*(20.0_ark*fea2444*sqrt(3.0_ark)+ &!
          fea1444*sqrt(3.0_ark)-3.0_ark*fea1455*sqrt(3.0_ark)+ 12.0_ark*fea2445)*y1*s4b**3/54.0_ark+ &!
          sqrt(3.0_ark)*(fea1222+ 2.0_ark*fea1112)*y1*y2**3/3.0_ark-sqrt(3.0_ark)*fea1122*y1**2*y3**2+ &!
          sqrt(3.0_ark)*(2.0_ark*fea1222+ fea1112)*y1**3*y2/3.0_ark+ sqrt(3.0_ark)*(fea112+ &!
          2.0_ark*fea122)*y1**2*y2/3.0_ark+ (-fea155-fea144/3.0_ark-8.0_ark/3.0_ark*fea244)*y1*s4a*s4b-&!
          sqrt(3.0_ark)*(2.0_ark*fea3355+ fea1155+ 3.0_ark*fea1144)*y3**2*s4b**2/6.0_ark+ (2.0_ark/3.0_ark*fea3355+ &!
          5.0_ark/6.0_ark*fea1155-fea1144/2.0_ark)*y3**2*s4a*s4b-sqrt(3.0_ark)*(fea224-fea114)*y3**2*s4a/3.0_ark+ &!
          (2.0_ark/3.0_ark*fea3355+ 5.0_ark/6.0_ark*fea1155-fea1144/2.0_ark)*y2**2*s4a*s4b-sqrt(3.0_ark)*(-fea144+ &!
          fea244)*y2*s4b**2/3.0_ark

     s1 = s3+ fea444*s4a**2*s4b+ (-fea1224*sqrt(3.0_ark)/6.0_ark+ fea1125/2.0_ark-fea1225/2.0_ark+ &!
          fea1124*sqrt(3.0_ark)/6.0_ark)*y2*y3**2*s4a-sqrt(3.0_ark)*(fea24-fea14)*y3*s4a/3.0_ark+ &!
          sqrt(3.0_ark)*(fea24-fea14)*y2*s4a/3.0_ark+ sqrt(3.0_ark)*(16.0_ark*fea2444*sqrt(3.0_ark)+ &!
          17.0_ark*fea1444*sqrt(3.0_ark)+ 3.0_ark*fea1455*sqrt(3.0_ark)-12.0_ark*fea2445)*y3*s4b**3/108.0_ark-&!
          sqrt(3.0_ark)*(fea122+ 2.0_ark*fea112)*y1*y3**2/3.0_ark+ (-2.0_ark*fea4444+ &!
          fea4455/3.0_ark)*s4a*s4b**3+ (fea24/3.0_ark+ 2.0_ark/3.0_ark*fea14)*y2*s4b+ &!
          sqrt(3.0_ark)*fea1122*y1**2*y2**2+ (fea224/3.0_ark+ 2.0_ark/3.0_ark*fea114)*y3**2*s4b+ &!
          sqrt(3.0_ark)*(3.0_ark*fea155+ fea144+ 2.0_ark*fea244)*y2*s4a**2/6.0_ark+ (4.0_ark/3.0_ark*fea224-&!
          fea114/3.0_ark)*y1**2*s4b+ s2

     s3 = s1+ (fea1224*sqrt(3.0_ark)/6.0_ark-fea1125/2.0_ark+ fea1225/2.0_ark-&!
          fea1124*sqrt(3.0_ark)/6.0_ark)*y2**2*y3*s4a-sqrt(3.0_ark)*fea12*y1*y3-sqrt(3.0_ark)*(3.0_ark*fea155+ &!
          fea144+ 2.0_ark*fea244)*y3*s4a**2/6.0_ark+ (fea1124*sqrt(3.0_ark)/3.0_ark-fea1224*sqrt(3.0_ark)/3.0_ark+ &!
          fea1225)*y1**2*y2*s4a+ (-fea1124*sqrt(3.0_ark)/3.0_ark+ fea1224*sqrt(3.0_ark)/3.0_ark-&!
          fea1225)*y1**2*y3*s4a-sqrt(3.0_ark)*(2.0_ark*fea1222+ fea1112)*y1**3*y3/3.0_ark+ &!
          (8.0_ark/3.0_ark*fea3355+ fea1155/3.0_ark+ fea1144)*y1**2*s4a*s4b+ sqrt(3.0_ark)*(fea1222-&!
          fea1112)*y2*y3**3/3.0_ark-sqrt(3.0_ark)*(-fea1444-3.0_ark*fea1455+ 4.0_ark*fea2444)*y3*s4a**3/12.0_ark-&!
          sqrt(3.0_ark)*(fea1222+ 2.0_ark*fea1112)*y1*y3**3/3.0_ark+ (fea155/2.0_ark-5.0_ark/6.0_ark*fea144-&!
          2.0_ark/3.0_ark*fea244)*y2*s4a*s4b


     s2 = s3+ (fea155/2.0_ark-5.0_ark/6.0_ark*fea144-2.0_ark/3.0_ark*fea244)*y3*s4a*s4b+ sqrt(3.0_ark)*(-&!
          fea144+ fea244)*y3*s4b**2/3.0_ark+ sqrt(3.0_ark)*fea12*y1*y2+ sqrt(3.0_ark)*fea1111*y2**4/2.0_ark+ &!
          (4.0_ark/3.0_ark*fea24-fea14/3.0_ark)*y1*s4b-sqrt(3.0_ark)*fea1*y3/2.0_ark-2.0_ark*fea44*s4a*s4b-&!
          sqrt(3.0_ark)*fea11*y3**2/2.0_ark+ sqrt(3.0_ark)*fea11*y2**2/2.0_ark-sqrt(3.0_ark)*(-fea124*sqrt(3.0_ark)+ &!
          2.0_ark*fea135)*y1*y3*s4b/3.0_ark-sqrt(3.0_ark)*(3.0_ark*fea1244+ 3.0_ark*fea1255+ &!
          fea1245*sqrt(3.0_ark))*y1*y3*s4b**2/6.0_ark-sqrt(3.0_ark)*(-fea124*sqrt(3.0_ark)+ &!
          2.0_ark*fea135)*y1*y2*s4b/3.0_ark

     s4 = s2+ sqrt(3.0_ark)*(4.0_ark*fea2444*sqrt(3.0_ark)-fea1444*sqrt(3.0_ark)-fea1455*sqrt(3.0_ark)-&!
          4.0_ark*fea2445)*y1*s4a**2*s4b/6.0_ark+ sqrt(3.0_ark)*(-fea1125+ fea1124*sqrt(3.0_ark)-fea1225+ &!
          fea1224*sqrt(3.0_ark))*y2**2*y3*s4b/6.0_ark+ sqrt(3.0_ark)*(-3.0_ark*fea1244-3.0_ark*fea1255+ &!
          fea1245*sqrt(3.0_ark))*y1*y3*s4a**2/6.0_ark+ sqrt(3.0_ark)*(fea1125+ fea1124*sqrt(3.0_ark)+ &!
          fea1225)*y1*y3**2*s4b/3.0_ark-sqrt(3.0_ark)*(-fea1124+ fea1125*sqrt(3.0_ark)+ &!
          fea1224)*y1*y3**2*s4a/3.0_ark

     s3 = s4+ sqrt(3.0_ark)*(fea1125+ fea1124*sqrt(3.0_ark)+ fea1225)*y1*y2**2*s4b/3.0_ark+ (-&!
          fea1245*sqrt(3.0_ark)/3.0_ark-fea1244+ fea1255)*y2*y3*s4a*s4b+ sqrt(3.0_ark)*(-fea1124+ &!
          fea1125*sqrt(3.0_ark)+ fea1224)*y1*y2**2*s4a/3.0_ark-sqrt(3.0_ark)*(-3.0_ark*fea1244-3.0_ark*fea1255+ &!
          fea1245*sqrt(3.0_ark))*y1*y2*s4a**2/6.0_ark-sqrt(3.0_ark)*fea1123*y1*y2*y3**2/2.0_ark+ &!
          sqrt(3.0_ark)*(3.0_ark*fea1244+ 3.0_ark*fea1255+ fea1245*sqrt(3.0_ark))*y1*y2*s4b**2/6.0_ark-&!
          sqrt(3.0_ark)*(4.0_ark*fea2444-7.0_ark*fea1444+ 3.0_ark*fea1455)*y3*s4a*s4b**2/12.0_ark

     t4y= s3+ sqrt(3.0_ark)*fea1123*y1*y2**2*y3/2.0_ark+ sqrt(3.0_ark)*(3.0_ark*fea1444*sqrt(3.0_ark)+ &!
          fea1455*sqrt(3.0_ark)+ 4.0_ark*fea2445)*y3*s4a**2*s4b/12.0_ark+ sqrt(3.0_ark)*(-fea1125+ &!
          fea1124*sqrt(3.0_ark)-fea1225+ fea1224*sqrt(3.0_ark))*y2*y3**2*s4b/6.0_ark+ &!
          (2.0_ark/3.0_ark*fea1245*sqrt(3.0_ark)-fea1244+ fea1255)*y1*y3*s4a*s4b+ sqrt(3.0_ark)*(4.0_ark*fea2444-&!
          7.0_ark*fea1444+ 3.0_ark*fea1455)*y2*s4a*s4b**2/12.0_ark+ sqrt(3.0_ark)*(3.0_ark*fea1444*sqrt(3.0_ark)+ &!
          fea1455*sqrt(3.0_ark)+ 4.0_ark*fea2445)*y2*s4a**2*s4b/12.0_ark+ sqrt(3.0_ark)*(fea1125+ &!
          fea1224*sqrt(3.0_ark)+ fea1225)*y1**2*y2*s4b/3.0_ark+ sqrt(3.0_ark)*(fea1125+ fea1224*sqrt(3.0_ark)+ &!
          fea1225)*y1**2*y3*s4b/3.0_ark+ fea4*s4b+ (2.0_ark/3.0_ark*fea1245*sqrt(3.0_ark)-fea1244+ &!
          fea1255)*y1*y2*s4a*s4b+ fea1234*y1*y2*y3*s4b+ sqrt(3.0_ark)*(fea135+ &!
          fea124*sqrt(3.0_ark))*y2*y3*s4b/3.0_ark+ fea444*s4b**3

     t56y = 0.0_ark
     if (parmax>149) then
       !
       s5 = hea333335*y3**5*s4b+ (2.0_ark/9.0_ark*fea111344*sqrt(3.0_ark)-hea233344/3.0_ark-&!
            fea233344*sqrt(3.0_ark)/3.0_ark+ fea133344*sqrt(3.0_ark)/9.0_ark-fea233345/3.0_ark)*y1*y2**3*s4b**2+ &!
            (17.0_ark/5.0_ark*hea355555-24.0_ark/25.0_ark*fea145555-9.0_ark/25.0_ark*hea344445-&!
            21.0_ark/25.0_ark*sqrt(3.0_ark)*hea344455-2.0_ark/25.0_ark*sqrt(3.0_ark)*hea345555-&!
            27.0_ark/25.0_ark*hea144445)*y1*s4b**5+ (-hea233334/2.0_ark-fea133334*sqrt(3.0_ark)/3.0_ark+&!
            fea233334*sqrt(3.0_ark)/3.0_ark+ fea233335/2.0_ark)*y1**4*y2*s4a+ hea333355*y3**4*s4b**2+ &!
            (7.0_ark*hea355555-9.0_ark/5.0_ark*hea344445-12.0_ark/5.0_ark*hea144445-&!
            6.0_ark/5.0_ark*sqrt(3.0_ark)*hea344455+ 4.0_ark/15.0_ark*sqrt(3.0_ark)*hea345555-&!
            4.0_ark/5.0_ark*fea145555)*y2*s4a**2*s4b**3-hea333355*y2**4*s4b**2+ hea233333*y2*y3**5

       s4 = s5+ sqrt(3.0_ark)*fea11111*y2**5/2.0_ark+ (sqrt(3.0_ark)*fea11122/2.0_ark-&!
            hea11122/2.0_ark)*y2**3*y3**2-hea11122*y1**3*y3**2-hea111333*y1**3*y2**3+ &!
            hea355555*y3*s4b**5+ (2.0_ark/3.0_ark*hea444445-8.0_ark/3.0_ark*fea444444)*s4a**3*s4b**3+ &!
            sqrt(3.0_ark)*fea111111*y2**6/2.0_ark+ hea444445*s4a**5*s4b+ (-hea444445/3.0_ark-&!
            8.0_ark/3.0_ark*fea444444)*s4a*s4b**5

       s3 = s4+ hea333444*y3**3*s4a**3+ (-2.0_ark*hea334445-5.0_ark/3.0_ark*hea114555+ &!
            2.0_ark/3.0_ark*hea334555+ 8.0_ark*fea335555-4.0_ark/3.0_ark*sqrt(3.0_ark)*hea334455)*y1**2*s4a**3*s4b+ &!
            (-4.0_ark/9.0_ark*fea133344*sqrt(3.0_ark)+ 4.0_ark/9.0_ark*fea111344*sqrt(3.0_ark)+ fea233345/3.0_ark+ &!
            hea233344/3.0_ark)*y2**3*y3*s4b**2+ (-2.0_ark/3.0_ark*fea12255*sqrt(3.0_ark)-2.0_ark*fea11245-&!
            3.0_ark*hea11244-3.0_ark*fea22344*sqrt(3.0_ark)-7.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)-&!
            4.0_ark*fea12245)*y2*y3**2*s4b**2+ hea222555*y3**3*s4b**3-hea11112*y1**4*y3+ &!
            hea333335*y2**5*s4b+ hea11112*y1**4*y2+ (hea223333+ &!
            sqrt(3.0_ark)*fea111122)*y1**2*y2**4+ hea11122*y1**3*y2**2+ hea222555*y2**3*s4b**3-&!
            hea135555*y1*y2*s4b**4-hea133444*y1*y2**2*s4a**3+ hea334445*y3**2*s4a**3*s4b+ (-&!
            hea333335/2.0_ark+ 3.0_ark/2.0_ark*fea333334)*y1**5*s4b-hea233444*y2**2*y3*s4a**3+ &!
            (sqrt(3.0_ark)*fea11122/2.0_ark+ hea11122/2.0_ark)*y1**2*y2**3+ (hea11112/2.0_ark-&!
            sqrt(3.0_ark)*fea11112/2.0_ark)*y2*y3**4

       s4 = s3+ (-hea11112/2.0_ark+ sqrt(3.0_ark)*fea11112/2.0_ark)*y2**4*y3-&!
            hea344455*y2*s4a**3*s4b**2-hea345555*y2*s4a*s4b**4-sqrt(3.0_ark)*fea111111*y3**6/2.0_ark+ &!
            hea334445*y2**2*s4a**3*s4b-hea233333*y2**5*y3+ hea223334*y2**2*y3**3*s4a-&!
            hea113344*y1**2*y2**2*s4a**2+ hea233334*y2*y3**4*s4a-hea233334*y2**4*y3*s4a-&!
            hea333455*y2**3*s4a*s4b**2+ hea223555*y2**2*y3*s4b**3+ &!
            hea334555*y2**2*s4a*s4b**3-hea334455*y2**2*s4a**2*s4b**2-&!
            hea113334*y1**2*y2**3*s4a-hea223334*y2**3*y3**2*s4a+ hea233444*y2*y3**2*s4a**3

       s5 = s4+ (4.0_ark*hea33335-3.0_ark*fea11114)*y1**4*s4b+ (fea233445+ hea233444-hea133444-&!
            fea133445+ fea113445)*y1**2*y2*s4a**3-hea122333*y1*y2**3*y3**2-&!
            hea133344*y1*y2**3*s4a**2-hea111334*y1**3*y2**2*s4a+ &!
            2.0_ark*fea12355*y1*y2*y3*s4a*s4b+ hea334455*y3**2*s4a**2*s4b**2+ &!
            hea22335*y2**2*y3**2*s4b+ (-hea223333-sqrt(3.0_ark)*fea111122)*y1**2*y3**4

       s6 = s5-sqrt(3.0_ark)*(2.0_ark*fea133344*sqrt(3.0_ark)+ 4.0_ark*fea111344*sqrt(3.0_ark)+ &!
            6.0_ark*hea133344+ 3.0_ark*fea233345)*y2**3*y3*s4a*s4b/9.0_ark+ (fea22255*sqrt(3.0_ark)/3.0_ark-&!
            fea11155*sqrt(3.0_ark)/3.0_ark-fea22245)*y3**3*s4b**2+ sqrt(3.0_ark)*(3.0_ark*hea233334+ &!
            fea133334*sqrt(3.0_ark)+ 2.0_ark*fea233334*sqrt(3.0_ark)+ 3.0_ark*fea233335)*y1*y2**4*s4b/9.0_ark-&!
            hea33444*y2**2*s4a**3

       s2 = s6+ 2.0_ark/3.0_ark*sqrt(3.0_ark)*(2.0_ark*hea123344+ sqrt(3.0_ark)*hea123345+ &!
            fea112344*sqrt(3.0_ark))*y1**2*y2*y3*s4a*s4b-sqrt(3.0_ark)*(-3.0_ark*fea233445-3.0_ark*hea233444+ &!
            4.0_ark*fea122444*sqrt(3.0_ark)+ 5.0_ark*fea113444*sqrt(3.0_ark)-12.0_ark*sqrt(3.0_ark)*hea223555+ &!
            9.0_ark*hea133444+ 9.0_ark*fea133445)*y1**2*y2*s4b**3/9.0_ark+ (-hea11112/2.0_ark-&!
            sqrt(3.0_ark)*fea11112/2.0_ark)*y1*y3**4+ fea123444*y1*y2*y3*s4b**3+ &!
            sqrt(3.0_ark)*(4.0_ark*fea133334*sqrt(3.0_ark)+ 3.0_ark*hea233334-fea233334*sqrt(3.0_ark)+ &!
            3.0_ark*fea233335)*y2*y3**4*s4b/9.0_ark+ (-2.0_ark/3.0_ark*fea233445-hea233444-&!
            2.0_ark/3.0_ark*fea122444*sqrt(3.0_ark)+ 2.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)+ 2.0_ark/3.0_ark*fea133445-&!
            2.0_ark/3.0_ark*fea113445)*y2*y3**2*s4a*s4b**2

       s5 = s2+ (-sqrt(3.0_ark)*fea11122/2.0_ark-hea11122/2.0_ark)*y1**2*y3**3-&!
            sqrt(3.0_ark)*(56.0_ark*sqrt(3.0_ark)*hea135555+ 69.0_ark*hea234555-9.0_ark*hea234445+ &!
            216.0_ark*fea124444+ 96.0_ark*fea235555)*y1*y2*s4a**2*s4b**2/72.0_ark+ sqrt(3.0_ark)*(-&!
            2.0_ark*fea111344*sqrt(3.0_ark)-3.0_ark*hea233344-3.0_ark*fea233344*sqrt(3.0_ark)-fea133344*sqrt(3.0_ark)+ &!
            3.0_ark*hea133344+ 3.0_ark*fea233345)*y1*y2**3*s4a*s4b/9.0_ark+ sqrt(3.0_ark)*(-&!
            3.0_ark*sqrt(3.0_ark)*hea222555+ 5.0_ark*hea333444+ 4.0_ark*fea222444*sqrt(3.0_ark)+ &!
            3.0_ark*hea333455)*y2**3*s4a**2*s4b/3.0_ark+ (-fea22255*sqrt(3.0_ark)/3.0_ark+ &!
            fea11155*sqrt(3.0_ark)/3.0_ark+ fea22245)*y2**3*s4b**2+ sqrt(3.0_ark)*(-5.0_ark*fea133445+ &!
            4.0_ark*hea233444+ 2.0_ark*fea233445-3.0_ark*fea122444*sqrt(3.0_ark)-2.0_ark*fea113444*sqrt(3.0_ark)+ &!
            6.0_ark*sqrt(3.0_ark)*hea223555-5.0_ark*hea133444+ 2.0_ark*fea113445)*y1*y2**2*s4a**2*s4b/3.0_ark+ (-&!
            2.0_ark*hea233333-sqrt(3.0_ark)*fea111112)*y1**5*y3+ sqrt(3.0_ark)*(fea223334*sqrt(3.0_ark)+ &!
            2.0_ark*hea113334+ 2.0_ark*hea111334)*y2**2*y3**3*s4b/3.0_ark

       s4 = s5+ (5.0_ark/3.0_ark*hea55555-2.0_ark/3.0_ark*fea45555)*s4a**4*s4b+ (hea11112/2.0_ark+ &!
            sqrt(3.0_ark)*fea11112/2.0_ark)*y1*y2**4-sqrt(3.0_ark)*(-6.0_ark*fea233445-6.0_ark*hea233444+ &!
            5.0_ark*fea122444*sqrt(3.0_ark)+ 4.0_ark*fea113444*sqrt(3.0_ark)-12.0_ark*sqrt(3.0_ark)*hea223555+ &!
            9.0_ark*hea133444+ 9.0_ark*fea133445)*y1*y3**2*s4b**3/9.0_ark+ (-sqrt(3.0_ark)*fea111112-&!
            hea233333)*y1*y3**5+ hea33444*y3**2*s4a**3+ (2.0_ark*hea233333+ &!
            sqrt(3.0_ark)*fea111112)*y1**5*y2+ (-sqrt(3.0_ark)*fea11122/2.0_ark+ hea11122/2.0_ark)*y2**2*y3**3+ &!
            hea33335*y3**4*s4b+ hea111333*y1**3*y3**3

       s5 = s4+ hea33335*y2**4*s4b-hea123344*y1*y2**2*y3*s4a**2+ &!
            sqrt(3.0_ark)*(4.0_ark*fea133334*sqrt(3.0_ark)+ 3.0_ark*hea233334-fea233334*sqrt(3.0_ark)+ &!
            3.0_ark*fea233335)*y2**4*y3*s4b/9.0_ark-sqrt(3.0_ark)*(2.0_ark*fea233445+ 2.0_ark*hea233444-&!
            2.0_ark*fea122444*sqrt(3.0_ark)-2.0_ark*fea113444*sqrt(3.0_ark)+ 3.0_ark*sqrt(3.0_ark)*hea223555-&!
            4.0_ark*hea133444-2.0_ark*fea133445+ 2.0_ark*fea113445)*y2**2*y3*s4a**2*s4b/3.0_ark+ (-&!
            2.0_ark*hea223333-sqrt(3.0_ark)*fea111122)*y1**4*y3**2-sqrt(3.0_ark)*(2.0_ark*fea233445+ &!
            2.0_ark*hea233444-2.0_ark*fea122444*sqrt(3.0_ark)-2.0_ark*fea113444*sqrt(3.0_ark)+ &!
            3.0_ark*sqrt(3.0_ark)*hea223555-4.0_ark*hea133444-2.0_ark*fea133445+ &!
            2.0_ark*fea113445)*y2*y3**2*s4a**2*s4b/3.0_ark+ sqrt(3.0_ark)*(5.0_ark*hea34555-6.0_ark*fea24455+ &!
            3.0_ark*fea24555*sqrt(3.0_ark)-4.0_ark*fea15555-3.0_ark*hea14445)*y2*s4b**4/4.0_ark+ &!
            (sqrt(3.0_ark)*fea111112+ hea233333)*y1*y2**5

       s3 = s5+ sqrt(3.0_ark)*(fea223334*sqrt(3.0_ark)+ 2.0_ark*hea113334+ &!
            2.0_ark*hea111334)*y2**3*y3**2*s4b/3.0_ark+ hea123344*y1*y2*y3**2*s4a**2+ &!
            hea123345*y1*y2*y3**2*s4a*s4b+ sqrt(3.0_ark)*(sqrt(3.0_ark)*hea223345+ 4.0_ark*hea113344-&!
            2.0_ark*fea223344*sqrt(3.0_ark))*y1**2*y3**2*s4a*s4b/6.0_ark+ sqrt(3.0_ark)*(3.0_ark*hea33444+ &!
            fea22455*sqrt(3.0_ark)+ fea22555)*y3**2*s4b**3/3.0_ark+ hea355555*y2*s4b**5-&!
            hea223333*y2**4*y3**2+ hea223333*y2**2*y3**4-sqrt(3.0_ark)*fea11223*y1**2*y2*y3**2+ &!
            fea123444*y1*y2*y3*s4a**2*s4b

       s5 = s3+ hea123335*y1*y2*y3**3*s4b-sqrt(3.0_ark)*(hea123335-fea123334)*y1*y2**3*y3*s4a/2.0_ark+ sqrt(3.0_ark)*(9.0_ark*hea13444+ &!
            2.0_ark*fea23455*sqrt(3.0_ark)+ fea12455*sqrt(3.0_ark)+ 9.0_ark*fea12555)*y1*y2*s4a**2*s4b/9.0_ark+ sqrt(3.0_ark)*(-fea122334+ hea122335)*y1**2*y2*y3**2*s4a/4.0_ark-&!
            sqrt(3.0_ark)*(10.0_ark*fea222444*sqrt(3.0_ark)+ 17.0_ark*hea333444-12.0_ark*sqrt(3.0_ark)*hea222555+ 9.0_ark*hea333455)*y1**3*s4a**2*s4b/6.0_ark-sqrt(3.0_ark)*fea11111*y3**5/2.0_ark+ &!
            2.0_ark*fea45555*s4a**2*s4b**3+ (3.0_ark/4.0_ark*fea122334+ hea122335/4.0_ark)*y1**2*y2*y3**2*s4b

       s6 = s5-sqrt(3.0_ark)*(3.0_ark*fea133445-hea233444-3.0_ark*fea233445+ &!
            3.0_ark*fea113444*sqrt(3.0_ark)-6.0_ark*sqrt(3.0_ark)*hea223555+ 5.0_ark*hea133444+ &!
            2.0_ark*fea122444*sqrt(3.0_ark))*y1**2*y2*s4a**2*s4b/3.0_ark-hea333444*y2**3*s4a**3+ &!
            sqrt(3.0_ark)*(-5.0_ark*fea133445+ 4.0_ark*hea233444+ 2.0_ark*fea233445-3.0_ark*fea122444*sqrt(3.0_ark)-&!
            2.0_ark*fea113444*sqrt(3.0_ark)+ 6.0_ark*sqrt(3.0_ark)*hea223555-5.0_ark*hea133444+ &!
            2.0_ark*fea113445)*y1*y3**2*s4a**2*s4b/3.0_ark-sqrt(3.0_ark)*(-fea22455*sqrt(3.0_ark)+ &!
            9.0_ark*hea33444-2.0_ark*fea11455*sqrt(3.0_ark)+ 9.0_ark*fea22555)*y3**2*s4a**2*s4b/9.0_ark

       s4 = s6+ sqrt(3.0_ark)*(4.0_ark*fea22255*sqrt(3.0_ark)-3.0_ark*fea22245+ &!
            2.0_ark*fea11155*sqrt(3.0_ark))*y3**3*s4a*s4b/9.0_ark+ sqrt(3.0_ark)*(hea113334+ hea223334+ &!
            fea223334*sqrt(3.0_ark))*y1**2*y2**3*s4b/3.0_ark+ sqrt(3.0_ark)*(6.0_ark*hea333355-&!
            sqrt(3.0_ark)*hea111145-6.0_ark*fea333344*sqrt(3.0_ark))*y2**4*s4a*s4b/6.0_ark+ (-3.0_ark*hea13444-&!
            2.0_ark/3.0_ark*fea23455*sqrt(3.0_ark)+ 2.0_ark/3.0_ark*fea12455*sqrt(3.0_ark))*y1*y2*s4a*s4b**2+ &!
            sqrt(3.0_ark)*(-2.0_ark*fea111344*sqrt(3.0_ark)-3.0_ark*hea233344-3.0_ark*fea233344*sqrt(3.0_ark)-&!
            fea133344*sqrt(3.0_ark)+ 3.0_ark*hea133344+ 3.0_ark*fea233345)*y1*y3**3*s4a*s4b/9.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*fea133445-hea233444-3.0_ark*fea233445+ 3.0_ark*fea113444*sqrt(3.0_ark)-&!
            6.0_ark*sqrt(3.0_ark)*hea223555+ 5.0_ark*hea133444+ &!
            2.0_ark*fea122444*sqrt(3.0_ark))*y1**2*y3*s4a**2*s4b/3.0_ark

       s6 = s4+ sqrt(3.0_ark)*(56.0_ark*sqrt(3.0_ark)*hea135555+ 69.0_ark*hea234555-9.0_ark*hea234445+ &!
            216.0_ark*fea124444+ 96.0_ark*fea235555)*y1*y3*s4a**2*s4b**2/72.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*hea234555+ hea234445+ 8.0_ark*fea124444)*y1*y3*s4a**4/8.0_ark+ &!
            sqrt(3.0_ark)*(3.0_ark*fea24555*sqrt(3.0_ark)+ 3.0_ark*hea34555-4.0_ark*fea24455-&!
            3.0_ark*hea14445)*y2*s4a**2*s4b**2/4.0_ark+ sqrt(3.0_ark)*(fea23455*sqrt(3.0_ark)+ 2.0_ark*fea12555+ &!
            6.0_ark*hea13444)*y2*y3*s4b**3/3.0_ark

       s5 = s6-sqrt(3.0_ark)*(-4.0_ark*fea12455*sqrt(3.0_ark)+ fea23455*sqrt(3.0_ark)+ 18.0_ark*fea12555+ &!
            18.0_ark*hea13444)*y2*y3*s4a**2*s4b/9.0_ark-sqrt(3.0_ark)*(-2.0_ark*fea133334*sqrt(3.0_ark)+ &!
            3.0_ark*hea233334-4.0_ark*fea233334*sqrt(3.0_ark)+ 3.0_ark*fea233335)*y1**4*y2*s4b/18.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*fea24555*sqrt(3.0_ark)+ 3.0_ark*hea34555-4.0_ark*fea24455-&!
            3.0_ark*hea14445)*y3*s4a**2*s4b**2/4.0_ark+ sqrt(3.0_ark)*(sqrt(3.0_ark)*hea223345+ 4.0_ark*hea113344-&!
            2.0_ark*fea223344*sqrt(3.0_ark))*y1**2*y2**2*s4a*s4b/6.0_ark-sqrt(3.0_ark)*(-&!
            2.0_ark*fea133334*sqrt(3.0_ark)+ 3.0_ark*hea233334-4.0_ark*fea233334*sqrt(3.0_ark)+ &!
            3.0_ark*fea233335)*y1**4*y3*s4b/18.0_ark

       s1 = s5+ sqrt(3.0_ark)*(hea233344-fea133344*sqrt(3.0_ark)+ hea133344-&!
            fea233344*sqrt(3.0_ark))*y1**3*y2*s4a*s4b/3.0_ark+ sqrt(3.0_ark)*(-hea223334+ &!
            fea223334*sqrt(3.0_ark)+ hea111334)*y1**3*y2**2*s4b/3.0_ark+ &!
            hea223345*y2**2*y3**2*s4a*s4b+ sqrt(3.0_ark)*(hea233344-fea133344*sqrt(3.0_ark)+ hea133344-&!
            fea233344*sqrt(3.0_ark))*y1**3*y3*s4a*s4b/3.0_ark+ sqrt(3.0_ark)*(-hea223334+ &!
            fea223334*sqrt(3.0_ark)+ hea111334)*y1**3*y3**2*s4b/3.0_ark-&!
            sqrt(3.0_ark)*fea11123*y1*y2*y3**3/2.0_ark-sqrt(3.0_ark)*(5.0_ark*fea12245+ 2.0_ark*fea11245+ &!
            6.0_ark*hea11244+ 4.0_ark*fea22344*sqrt(3.0_ark)+ &!
            2.0_ark*fea11255*sqrt(3.0_ark))*y1*y2**2*s4a*s4b/3.0_ark-sqrt(3.0_ark)*(3.0_ark*hea13444-&!
            fea12455*sqrt(3.0_ark)+ fea12555)*y1*y2*s4b**3/3.0_ark+ hea123345*y1*y2**2*y3*s4a*s4b-&!
            sqrt(3.0_ark)*(5.0_ark*fea12245+ 2.0_ark*fea11245+ 6.0_ark*hea11244+ 4.0_ark*fea22344*sqrt(3.0_ark)+ &!
            2.0_ark*fea11255*sqrt(3.0_ark))*y1*y3**2*s4a*s4b/3.0_ark

       s5 = s1-sqrt(3.0_ark)*(2.0_ark*fea22235*sqrt(3.0_ark)-3.0_ark*fea11124-fea11135*sqrt(3.0_ark)+ &!
            3.0_ark*hea13335)*y1*y3**3*s4a/9.0_ark-sqrt(3.0_ark)*(-6.0_ark*fea233445-6.0_ark*hea233444+ &!
            5.0_ark*fea122444*sqrt(3.0_ark)+ 4.0_ark*fea113444*sqrt(3.0_ark)-12.0_ark*sqrt(3.0_ark)*hea223555+ &!
            9.0_ark*hea133444+ 9.0_ark*fea133445)*y1*y2**2*s4b**3/9.0_ark-sqrt(3.0_ark)*(-6.0_ark*fea11124+ &!
            fea11135*sqrt(3.0_ark)-2.0_ark*fea22235*sqrt(3.0_ark)+ 6.0_ark*hea13335)*y1**3*y3*s4a/9.0_ark+ &!
            sqrt(3.0_ark)*(-6.0_ark*fea11124+ fea11135*sqrt(3.0_ark)-2.0_ark*fea22235*sqrt(3.0_ark)+ &!
            6.0_ark*hea13335)*y1**3*y2*s4a/9.0_ark+ 2.0_ark/9.0_ark*sqrt(3.0_ark)*(3.0_ark*fea22245+ &!
            2.0_ark*fea22255*sqrt(3.0_ark)+ fea11155*sqrt(3.0_ark))*y1**3*s4a*s4b+ &!
            sqrt(3.0_ark)*(18.0_ark*hea33444+ 4.0_ark*fea22455*sqrt(3.0_ark)-fea11455*sqrt(3.0_ark)+ &!
            18.0_ark*fea22555)*y1**2*s4a**2*s4b/9.0_ark+ sqrt(3.0_ark)*fea11123*y1*y2**3*y3/2.0_ark+ &!
            sqrt(3.0_ark)*(2.0_ark*fea22235*sqrt(3.0_ark)-3.0_ark*fea11124-fea11135*sqrt(3.0_ark)+ &!
            3.0_ark*hea13335)*y1*y2**3*s4a/9.0_ark

       s4 = s5+ (2.0_ark/3.0_ark*fea233445+ hea233444+ 2.0_ark/3.0_ark*fea122444*sqrt(3.0_ark)-&!
            2.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)-2.0_ark/3.0_ark*fea133445+ &!
            2.0_ark/3.0_ark*fea113445)*y2**2*y3*s4a*s4b**2+ sqrt(3.0_ark)*fea123333*y1*y2*y3**4-&!
            sqrt(3.0_ark)*(2.0_ark*fea22555+ 6.0_ark*hea33444-fea11455*sqrt(3.0_ark))*y1**2*s4b**3/3.0_ark-&!
            sqrt(3.0_ark)*(3.0_ark*hea13444-fea12455*sqrt(3.0_ark)+ fea12555)*y1*y3*s4b**3/3.0_ark+ &!
            sqrt(3.0_ark)*(9.0_ark*hea13444+ 2.0_ark*fea23455*sqrt(3.0_ark)+ fea12455*sqrt(3.0_ark)+ &!
            9.0_ark*fea12555)*y1*y3*s4a**2*s4b/9.0_ark+ sqrt(3.0_ark)*(4.0_ark*fea22255*sqrt(3.0_ark)-&!
            3.0_ark*fea22245+ 2.0_ark*fea11155*sqrt(3.0_ark))*y2**3*s4a*s4b/9.0_ark+ &!
            sqrt(3.0_ark)*(fea22235*sqrt(3.0_ark)-6.0_ark*fea11124+ 4.0_ark*fea11135*sqrt(3.0_ark)+ &!
            6.0_ark*hea13335)*y2**3*y3*s4a/9.0_ark-sqrt(3.0_ark)*(fea22235*sqrt(3.0_ark)-6.0_ark*fea11124+ &!
            4.0_ark*fea11135*sqrt(3.0_ark)+ 6.0_ark*hea13335)*y2*y3**3*s4a/9.0_ark-sqrt(3.0_ark)*(-&!
            fea22455*sqrt(3.0_ark)+ 9.0_ark*hea33444-2.0_ark*fea11455*sqrt(3.0_ark)+ &!
            9.0_ark*fea22555)*y2**2*s4a**2*s4b/9.0_ark

       s5 = s4+ (-2.0_ark*hea22335+ 3.0_ark*fea11224)*y1**2*y2**2*s4b-sqrt(3.0_ark)*(hea334555-&!
            3.0_ark*hea334445+ 8.0_ark*fea335555-2.0_ark*sqrt(3.0_ark)*hea334455-&!
            2.0_ark*hea114555)*y3**2*s4b**4/4.0_ark+ hea122335*y1*y2**2*y3**2*s4b+ (-hea123344/3.0_ark-&!
            sqrt(3.0_ark)*hea123345/3.0_ark-2.0_ark/3.0_ark*fea112344*sqrt(3.0_ark))*y1*y2*y3**2*s4b**2-&!
            hea233344*y2**3*y3*s4a**2+ (-17.0_ark/8.0_ark*hea234555-3.0_ark/8.0_ark*hea234445-&!
            9.0_ark*fea124444-5.0_ark/3.0_ark*sqrt(3.0_ark)*hea135555-2.0_ark*fea235555)*y1*y2*s4a*s4b**3+ &!
            hea122333*y1*y2**2*y3**3+ (4.0_ark*hea12335-3.0_ark*fea11234)*y1**2*y2*y3*s4b

       s6 = s5+ (5.0_ark/8.0_ark*hea234445-sqrt(3.0_ark)*hea135555/3.0_ark+ 7.0_ark/8.0_ark*hea234555+ &!
            3.0_ark*fea124444+ 2.0_ark*fea235555)*y1*y2*s4a**3*s4b-sqrt(3.0_ark)*(-fea122334+ &!
            hea122335)*y1**2*y2**2*y3*s4a/4.0_ark+ (13.0_ark*hea355555-21.0_ark/5.0_ark*hea344445-&!
            18.0_ark/5.0_ark*hea144445-9.0_ark/5.0_ark*sqrt(3.0_ark)*hea344455-14.0_ark/15.0_ark*sqrt(3.0_ark)*hea345555-&!
            16.0_ark/5.0_ark*fea145555)*y1*s4a**2*s4b**3+ sqrt(3.0_ark)*(6.0_ark*hea333355-&!
            sqrt(3.0_ark)*hea111145-6.0_ark*fea333344*sqrt(3.0_ark))*y3**4*s4a*s4b/6.0_ark

       s3 = s6+ sqrt(3.0_ark)*(-3.0_ark*sqrt(3.0_ark)*hea222555+ 5.0_ark*hea333444+ &!
            4.0_ark*fea222444*sqrt(3.0_ark)+ 3.0_ark*hea333455)*y3**3*s4a**2*s4b/3.0_ark+ &!
            sqrt(3.0_ark)*(8.0_ark*fea22344*sqrt(3.0_ark)+ 6.0_ark*fea11255*sqrt(3.0_ark)+ 4.0_ark*fea12255*sqrt(3.0_ark)+ &!
            7.0_ark*fea11245+ 10.0_ark*fea12245+ 6.0_ark*hea11244)*y1**2*y3*s4a*s4b/3.0_ark+ &!
            hea223555*y2*y3**2*s4b**3+ sqrt(3.0_ark)*(8.0_ark*fea22344*sqrt(3.0_ark)+ &!
            6.0_ark*fea11255*sqrt(3.0_ark)+ 4.0_ark*fea12255*sqrt(3.0_ark)+ 7.0_ark*fea11245+ 10.0_ark*fea12245+ &!
            6.0_ark*hea11244)*y1**2*y2*s4a*s4b/3.0_ark+ hea123335*y1*y2**3*y3*s4b+ &!
            (fea11135*sqrt(3.0_ark)/3.0_ark+ fea11124+ fea22235*sqrt(3.0_ark)/3.0_ark)*y2**3*y3*s4b

       s5 = s3+ sqrt(3.0_ark)*(-fea11114+ hea33335)*y2**4*s4a+ (fea12255*sqrt(3.0_ark)+ &!
            fea11255*sqrt(3.0_ark)+ hea11244+ 2.0_ark*fea22344*sqrt(3.0_ark)+ 2.0_ark*fea11245+ &!
            2.0_ark*fea12245)*y1*y3**2*s4a**2+ (-fea133334*sqrt(3.0_ark)/3.0_ark+ fea233334*sqrt(3.0_ark)/3.0_ark-&!
            fea233335)*y1*y2**4*s4a+ (-5.0_ark/3.0_ark*fea233445-hea233444+ &!
            10.0_ark/3.0_ark*fea122444*sqrt(3.0_ark)+ 8.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)+ 11.0_ark/3.0_ark*fea133445+ &!
            fea113445/3.0_ark-6.0_ark*sqrt(3.0_ark)*hea223555+ 3.0_ark*hea133444)*y1**2*y3*s4a*s4b**2+ &!
            hea114555*y1**2*s4a*s4b**3+ (sqrt(3.0_ark)*hea223345/6.0_ark+ hea113344/3.0_ark+ &!
            fea223344*sqrt(3.0_ark)/3.0_ark)*y1**2*y3**2*s4b**2+ (2.0_ark*fea22344+ 2.0_ark*fea11255+ &!
            2.0_ark*fea12255+ fea11245*sqrt(3.0_ark)/3.0_ark+ fea12245*sqrt(3.0_ark)/3.0_ark)*y2**2*y3*s4a*s4b+ &!
            (5.0_ark/8.0_ark*hea234445-sqrt(3.0_ark)*hea135555/3.0_ark+ 7.0_ark/8.0_ark*hea234555+ 3.0_ark*fea124444+ &!
            2.0_ark*fea235555)*y1*y3*s4a**3*s4b

       s6 = s5+ (fea12255*sqrt(3.0_ark)/3.0_ark+ 5.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)+ &!
            2.0_ark*fea22344*sqrt(3.0_ark)+ 2.0_ark*fea11245+ 3.0_ark*fea12245+ &!
            3.0_ark*hea11244)*y1*y2**2*s4b**2+ (5.0_ark/2.0_ark*hea14445+ 4.0_ark*fea15555-&!
            7.0_ark/2.0_ark*hea34555-11.0_ark/6.0_ark*fea24555*sqrt(3.0_ark)+ 4.0_ark*fea24455)*y2*s4a**3*s4b+ &!
            (hea123344/3.0_ark+ sqrt(3.0_ark)*hea123345/3.0_ark+ &!
            2.0_ark/3.0_ark*fea112344*sqrt(3.0_ark))*y1*y2**2*y3*s4b**2+ (7.0_ark*hea355555-&!
            9.0_ark/5.0_ark*hea344445-12.0_ark/5.0_ark*hea144445-6.0_ark/5.0_ark*sqrt(3.0_ark)*hea344455+ &!
            4.0_ark/15.0_ark*sqrt(3.0_ark)*hea345555-4.0_ark/5.0_ark*fea145555)*y3*s4a**2*s4b**3

       s4 = s6+ (-2.0_ark/9.0_ark*fea111344*sqrt(3.0_ark)+ hea233344/3.0_ark+ fea233344*sqrt(3.0_ark)/3.0_ark-&!
            fea133344*sqrt(3.0_ark)/9.0_ark+ fea233345/3.0_ark)*y1*y3**3*s4b**2+ (-2.0_ark/3.0_ark*fea113445-&!
            10.0_ark/3.0_ark*fea133445+ 4.0_ark/3.0_ark*fea233445-3.0_ark*hea133444+ 2.0_ark*hea233444+ &!
            6.0_ark*sqrt(3.0_ark)*hea223555-10.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)-&!
            8.0_ark/3.0_ark*fea122444*sqrt(3.0_ark))*y1*y2**2*s4a*s4b**2+ (-sqrt(3.0_ark)*fea122333/2.0_ark-&!
            hea122333/2.0_ark)*y1**2*y2**3*y3+ (5.0_ark/2.0_ark*hea14445+ 4.0_ark*fea15555-&!
            7.0_ark/2.0_ark*hea34555-11.0_ark/6.0_ark*fea24555*sqrt(3.0_ark)+ 4.0_ark*fea24455)*y3*s4a**3*s4b+ &!
            (5.0_ark/3.0_ark*fea233445+ hea233444-10.0_ark/3.0_ark*fea122444*sqrt(3.0_ark)-&!
            8.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)-11.0_ark/3.0_ark*fea133445-fea113445/3.0_ark+ &!
            6.0_ark*sqrt(3.0_ark)*hea223555-3.0_ark*hea133444)*y1**2*y2*s4a*s4b**2

       s5 = s4-sqrt(3.0_ark)*fea123333*y1*y2**4*y3+ hea234445*y2*y3*s4a**3*s4b+ &!
            (hea233334/2.0_ark+ fea133334*sqrt(3.0_ark)/3.0_ark-fea233334*sqrt(3.0_ark)/3.0_ark-&!
            fea233335/2.0_ark)*y1**4*y3*s4a+ (2.0_ark/9.0_ark*fea133344*sqrt(3.0_ark)+ fea233345/3.0_ark-&!
            fea233344*sqrt(3.0_ark)/3.0_ark+ hea233344/3.0_ark+ fea111344*sqrt(3.0_ark)/9.0_ark)*y1**3*y2*s4b**2+ &!
            (fea11135*sqrt(3.0_ark)/3.0_ark+ fea11124+ fea22235*sqrt(3.0_ark)/3.0_ark)*y2*y3**3*s4b+ (-&!
            3.0_ark*hea33444-2.0_ark/3.0_ark*fea22455*sqrt(3.0_ark)+ &!
            2.0_ark/3.0_ark*fea11455*sqrt(3.0_ark))*y2**2*s4a*s4b**2+ (2.0_ark/3.0_ark*fea113445+ &!
            10.0_ark/3.0_ark*fea133445-4.0_ark/3.0_ark*fea233445+ 3.0_ark*hea133444-2.0_ark*hea233444-&!
            6.0_ark*sqrt(3.0_ark)*hea223555+ 10.0_ark/3.0_ark*fea113444*sqrt(3.0_ark)+ &!
            8.0_ark/3.0_ark*fea122444*sqrt(3.0_ark))*y1*y3**2*s4a*s4b**2+ sqrt(3.0_ark)*(3.0_ark*hea234555+ &!
            hea234445+ 8.0_ark*fea124444)*y1*y2*s4a**4/8.0_ark+ (hea11244+ fea22344*sqrt(3.0_ark)+ &!
            fea11255*sqrt(3.0_ark)+ fea11245+ fea12245)*y2**2*y3*s4a**2

       s2 = s5+ hea234555*y2*y3*s4a*s4b**3+ sqrt(3.0_ark)*(12.0_ark*hea34555-15.0_ark*fea24455+ &!
            5.0_ark*fea24555*sqrt(3.0_ark)-6.0_ark*fea15555-9.0_ark*hea14445)*y2*s4a**4/18.0_ark+ &!
            hea144445*y1*s4a**4*s4b+ (3.0_ark*hea33444+ 2.0_ark/3.0_ark*fea22455*sqrt(3.0_ark)-&!
            2.0_ark/3.0_ark*fea11455*sqrt(3.0_ark))*y3**2*s4a*s4b**2+ (fea24555*sqrt(3.0_ark)+ &!
            hea34555)*y1*s4a*s4b**3+ (hea122333/2.0_ark-sqrt(3.0_ark)*fea122333/2.0_ark)*y1**3*y2**2*y3-&!
            sqrt(3.0_ark)*(5.0_ark*hea34555-6.0_ark*fea24455+ 3.0_ark*fea24555*sqrt(3.0_ark)-4.0_ark*fea15555-&!
            3.0_ark*hea14445)*y3*s4b**4/4.0_ark+ hea55555*s4b**5-sqrt(3.0_ark)*(fea22255-&!
            fea11155)*y2**3*s4a**2/3.0_ark+ (2.0_ark/3.0_ark*fea11135*sqrt(3.0_ark)-fea11124+ &!
            2.0_ark/3.0_ark*fea22235*sqrt(3.0_ark)+ 2.0_ark*hea13335)*y1**3*y3*s4b

       s5 = s2+ (-17.0_ark/8.0_ark*hea234555-3.0_ark/8.0_ark*hea234445-9.0_ark*fea124444-&!
            5.0_ark/3.0_ark*sqrt(3.0_ark)*hea135555-2.0_ark*fea235555)*y1*y3*s4a*s4b**3+ sqrt(3.0_ark)*(-&!
            9.0_ark*hea334445-14.0_ark*hea114555+ 11.0_ark*hea334555+ 24.0_ark*fea335555-&!
            10.0_ark*sqrt(3.0_ark)*hea334455)*y2**2*s4a**4/36.0_ark-sqrt(3.0_ark)*(-3.0_ark*hea333444-&!
            6.0_ark*fea222444*sqrt(3.0_ark)+ 4.0_ark*sqrt(3.0_ark)*hea222555-&!
            3.0_ark*hea333455)*y1**3*s4b**3/6.0_ark+ hea113334*y1**2*y3**3*s4a+ (-fea233445-&!
            hea233444+ hea133444+ fea133445-fea113445)*y1**2*y3*s4a**3+ sqrt(3.0_ark)*(hea334555-&!
            3.0_ark*hea334445+ 8.0_ark*fea335555-2.0_ark*sqrt(3.0_ark)*hea334455-&!
            2.0_ark*hea114555)*y2**2*s4b**4/4.0_ark-sqrt(3.0_ark)*(-9.0_ark*hea334445-14.0_ark*hea114555+ &!
            11.0_ark*hea334555+ 24.0_ark*fea335555-10.0_ark*sqrt(3.0_ark)*hea334455)*y3**2*s4a**4/36.0_ark+ (-&!
            fea12255*sqrt(3.0_ark)/3.0_ark-5.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)-2.0_ark*fea22344*sqrt(3.0_ark)-&!
            2.0_ark*fea11245-3.0_ark*fea12245-3.0_ark*hea11244)*y1*y3**2*s4b**2

       s4 = s5+ (-fea12255*sqrt(3.0_ark)-fea11255*sqrt(3.0_ark)-hea11244-2.0_ark*fea22344*sqrt(3.0_ark)-&!
            2.0_ark*fea11245-2.0_ark*fea12245)*y1*y2**2*s4a**2+ (-3.0_ark*hea11244-&!
            4.0_ark*fea22344*sqrt(3.0_ark)-8.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)-4.0_ark/3.0_ark*fea12255*sqrt(3.0_ark)-&!
            3.0_ark*fea11245-4.0_ark*fea12245)*y1**2*y2*s4b**2+ (3.0_ark*hea11244+ &!
            4.0_ark*fea22344*sqrt(3.0_ark)+ 8.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)+ 4.0_ark/3.0_ark*fea12255*sqrt(3.0_ark)+ &!
            3.0_ark*fea11245+ 4.0_ark*fea12245)*y1**2*y3*s4b**2+ (-hea123335/2.0_ark+ &!
            3.0_ark/2.0_ark*fea123334)*y1**3*y2*y3*s4b+ (2.0_ark/3.0_ark*fea12255*sqrt(3.0_ark)+ 2.0_ark*fea11245+ &!
            3.0_ark*hea11244+ 3.0_ark*fea22344*sqrt(3.0_ark)+ 7.0_ark/3.0_ark*fea11255*sqrt(3.0_ark)+ &!
            4.0_ark*fea12245)*y2**2*y3*s4b**2+ sqrt(3.0_ark)*(fea22255-fea11155)*y3**3*s4a**2/3.0_ark+ &!
            hea135555*y1*y3*s4b**4+ sqrt(3.0_ark)*(3.0_ark*hea33444+ fea22455*sqrt(3.0_ark)+ &!
            fea22555)*y2**2*s4b**3/3.0_ark-sqrt(3.0_ark)*(hea333335-fea333334)*y2**5*s4a/2.0_ark

       s3 = s4-sqrt(3.0_ark)*(12.0_ark*hea34555-15.0_ark*fea24455+ 5.0_ark*fea24555*sqrt(3.0_ark)-&!
            6.0_ark*fea15555-9.0_ark*hea14445)*y3*s4a**4/18.0_ark+ hea133344*y1*y3**3*s4a**2+ &!
            hea13444*y1*y3*s4a**3-hea13444*y1*y2*s4a**3+ hea13335*y1*y3**3*s4b+ &!
            hea111334*y1**3*y3**2*s4a+ hea133444*y1*y3**2*s4a**3+ &!
            hea113344*y1**2*y3**2*s4a**2+ hea344455*y3*s4a**3*s4b**2+ &!
            hea233344*y2*y3**3*s4a**2+ (2.0_ark*hea223333+ sqrt(3.0_ark)*fea111122)*y1**4*y2**2+ &!
            hea111145*y1**4*s4a*s4b+ hea11244*y1**2*y2*s4a**2+ hea345555*y3*s4a*s4b**4+ hea344445*y2*s4a**4*s4b+ hea34555*y3*s4a*s4b**3+ (3.0_ark/4.0_ark*fea122334+ &!
            hea122335/4.0_ark)*y1**2*y2**2*y3*s4b+ hea333455*y3**3*s4a*s4b**2

       s5 = s3+ (fea133334*sqrt(3.0_ark)/3.0_ark-fea233334*sqrt(3.0_ark)/3.0_ark+ fea233335)*y1*y3**4*s4a-&!
            sqrt(3.0_ark)*(hea111145+ 2.0_ark*fea333344)*y2**4*s4a**2/2.0_ark+ hea344445*y3*s4a**4*s4b+ &!
            hea34555*y2*s4a*s4b**3+ hea12335*y1*y2**2*y3*s4b+ (4.0_ark/9.0_ark*fea133344*sqrt(3.0_ark)-&!
            4.0_ark/9.0_ark*fea111344*sqrt(3.0_ark)-fea233345/3.0_ark-hea233344/3.0_ark)*y2*y3**3*s4b**2+ &!
            hea13335*y1*y2**3*s4b+ (2.0_ark/3.0_ark*fea11135*sqrt(3.0_ark)-fea11124+ &!
            2.0_ark/3.0_ark*fea22235*sqrt(3.0_ark)+ 2.0_ark*hea13335)*y1**3*y2*s4b

       s4 = s5+ hea14445*y1*s4a**3*s4b+ (2.0_ark*fea22344+ 2.0_ark*fea11255+ 2.0_ark*fea12255+ &!
            fea11245*sqrt(3.0_ark)/3.0_ark+ fea12245*sqrt(3.0_ark)/3.0_ark)*y2*y3**2*s4a*s4b+ &!
            (fea133344*sqrt(3.0_ark)/3.0_ark-fea111344*sqrt(3.0_ark)/3.0_ark-hea133344-&!
            fea233345)*y1**3*y2*s4a**2+ (-hea11244-fea22344*sqrt(3.0_ark)-fea11255*sqrt(3.0_ark)-&!
            fea11245-fea12245)*y2*y3**2*s4a**2-sqrt(3.0_ark)*(-fea11114+ hea33335)*y3**4*s4a+ &!
            (sqrt(3.0_ark)*fea122333/2.0_ark+ hea122333/2.0_ark)*y1**2*y2*y3**3+ (-sqrt(3.0_ark)*hea223345/6.0_ark-&!
            hea113344/3.0_ark-fea223344*sqrt(3.0_ark)/3.0_ark)*y1**2*y2**2*s4b**2+ &!
            sqrt(3.0_ark)*(3.0_ark*hea233334+ fea133334*sqrt(3.0_ark)+ 2.0_ark*fea233334*sqrt(3.0_ark)+ &!
            3.0_ark*fea233335)*y1*y3**4*s4b/9.0_ark+ hea334555*y3**2*s4a*s4b**3-hea11244*y1**2*y3*s4a**2

       s5 = s4+ sqrt(3.0_ark)*(hea333335-fea333334)*y3**5*s4a/2.0_ark-sqrt(3.0_ark)*(28.0_ark*hea144445+ &!
            26.0_ark*hea344445+ 19.0_ark*sqrt(3.0_ark)*hea344455+ 3.0_ark*sqrt(3.0_ark)*hea345555-&!
            90.0_ark*hea355555+ 36.0_ark*fea145555)*y2*s4a**5/75.0_ark+ hea12335*y1*y2*y3**2*s4b-&!
            sqrt(3.0_ark)*(hea12335-fea11234)*y1*y2*y3**2*s4a+ sqrt(3.0_ark)*(hea22335-&!
            fea11224)*y1**2*y3**2*s4a+ sqrt(3.0_ark)*(hea12335-fea11234)*y1*y2**2*y3*s4a-&!
            sqrt(3.0_ark)*(2.0_ark*fea133344*sqrt(3.0_ark)+ 4.0_ark*fea111344*sqrt(3.0_ark)+ 6.0_ark*hea133344+ &!
            3.0_ark*fea233345)*y2*y3**3*s4a*s4b/9.0_ark+ sqrt(3.0_ark)*(hea111145+ &!
            2.0_ark*fea333344)*y3**4*s4a**2/2.0_ark+ (-2.0_ark*hea22335+ 3.0_ark*fea11224)*y1**2*y3**2*s4b

       s6 = s5+ (-hea122333/2.0_ark+ sqrt(3.0_ark)*fea122333/2.0_ark)*y1**3*y2*y3**2+ &!
            sqrt(3.0_ark)*(28.0_ark*hea144445+ 26.0_ark*hea344445+ 19.0_ark*sqrt(3.0_ark)*hea344455+ &!
            3.0_ark*sqrt(3.0_ark)*hea345555-90.0_ark*hea355555+ 36.0_ark*fea145555)*y3*s4a**5/75.0_ark+ &!
            (3.0_ark*hea13444+ 2.0_ark/3.0_ark*fea23455*sqrt(3.0_ark)-&!
            2.0_ark/3.0_ark*fea12455*sqrt(3.0_ark))*y1*y3*s4a*s4b**2+ (-2.0_ark/9.0_ark*fea133344*sqrt(3.0_ark)-&!
            fea233345/3.0_ark+ fea233344*sqrt(3.0_ark)/3.0_ark-hea233344/3.0_ark-&!
            fea111344*sqrt(3.0_ark)/9.0_ark)*y1**3*y3*s4b**2

       t56y= s6+ (-fea133344*sqrt(3.0_ark)/3.0_ark+ fea111344*sqrt(3.0_ark)/3.0_ark+ hea133344+ &!
            fea233345)*y1**3*y3*s4a**2-sqrt(3.0_ark)*(-3.0_ark*fea233445-3.0_ark*hea233444+ &!
            4.0_ark*fea122444*sqrt(3.0_ark)+ 5.0_ark*fea113444*sqrt(3.0_ark)-12.0_ark*sqrt(3.0_ark)*hea223555+ &!
            9.0_ark*hea133444+ 9.0_ark*fea133445)*y1**2*y3*s4b**3/9.0_ark+ sqrt(3.0_ark)*(hea113334+ &!
            hea223334+ fea223334*sqrt(3.0_ark))*y1**2*y3**3*s4b/3.0_ark-sqrt(3.0_ark)*(hea22335-&!
            fea11224)*y1**2*y2**2*s4a+ sqrt(3.0_ark)*fea11223*y1**2*y2**2*y3+ sqrt(3.0_ark)*(hea123335-&!
          fea123334)*y1*y2*y3**3*s4a/2.0_ark
     endif
     !
     f = ( t4y+t56y )
     !
   end select
   !
   f = f*sinrho
   !
end function xy3_dms_E



end module xy3_alpha
