diff --git a/bin/genetic b/bin/genetic index f3c4309..a64e8e0 100755 Binary files a/bin/genetic and b/bin/genetic differ diff --git a/obj/adia.o b/obj/adia.o index 6df3517..5e668a2 100644 Binary files a/obj/adia.o and b/obj/adia.o differ diff --git a/obj/adia_mod.mod b/obj/adia_mod.mod index 2797744..2c5f176 100644 Binary files a/obj/adia_mod.mod and b/obj/adia_mod.mod differ diff --git a/obj/data_matrix.mod b/obj/data_matrix.mod index 63679ba..d2b214d 100644 Binary files a/obj/data_matrix.mod and b/obj/data_matrix.mod differ diff --git a/obj/data_transform.o b/obj/data_transform.o index d9425c9..e144e00 100644 Binary files a/obj/data_transform.o and b/obj/data_transform.o differ diff --git a/obj/diab_mod.mod b/obj/diab_mod.mod index d740a78..d0d467c 100644 Binary files a/obj/diab_mod.mod and b/obj/diab_mod.mod differ diff --git a/obj/keys.o b/obj/keys.o index 1db6f2e..5c77517 100644 Binary files a/obj/keys.o and b/obj/keys.o differ diff --git a/obj/keys_mod.mod b/obj/keys_mod.mod index c671ccc..eec6610 100644 Binary files a/obj/keys_mod.mod and b/obj/keys_mod.mod differ diff --git a/obj/matrix_form.o b/obj/matrix_form.o index 789292f..2a428b9 100644 Binary files a/obj/matrix_form.o and b/obj/matrix_form.o differ diff --git a/obj/model.o b/obj/model.o index faa19d4..c2bb1bd 100644 Binary files a/obj/model.o and b/obj/model.o differ diff --git a/obj/write.o b/obj/write.o index 653f598..19ebce7 100644 Binary files a/obj/write.o and b/obj/write.o differ diff --git a/obj/write_mod.mod b/obj/write_mod.mod index 55055f1..52577ac 100644 Binary files a/obj/write_mod.mod and b/obj/write_mod.mod differ diff --git a/src/.dim_parameter.f.swp b/src/.dim_parameter.f.swp deleted file mode 100644 index 317a9f9..0000000 Binary files a/src/.dim_parameter.f.swp and /dev/null differ diff --git a/src/model/Maik_chngs/maik_ctrans.f90 b/src/model/Maik_chngs/maik_ctrans.f90 deleted file mode 100644 index 0115a07..0000000 --- a/src/model/Maik_chngs/maik_ctrans.f90 +++ /dev/null @@ -1,673 +0,0 @@ -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % SUBROUTINE CTRANS(...) -! % -! % M. Vossel 21.03.2023 -! % -! % Routine to transform symmetryinput coordinates to symmetrized -! % coordinates. Distances Are discribet by Morse coordinates or -! % TMC depending on Set Parameters in the Genetic Input. -! % -! % input variables -! % q: -! % q(1): H1x -! % q(2): y -! % q(3): z -! % q(4): H2x -! % q(5): y -! % q(6): z -! % q(7): H3x -! % q(8): y -! % q(9): z -! -! -! -! % Internal variables: -! % t: primitive coordinates (double[qn]) -! % t(1): -! % t(2): -! % t(3): -! % t(4): -! % t(5): -! % t(6): -! % t(7): -! % t(8): -! % t(9): -! % t: dummy (double[qn]) -! % p: parameter vector -! % npar: length of parameter vector -! % -! % Output variables -! % s: symmetrized coordinates (double[qn]) -! % s(1): CH-symetric streatch -! % s(2): CH-asymetric streatch-ex -! % s(3): CH-asymetric streatch-ey -! % s(4): CH-bend-ex -! % s(5): CH-bend-ey -! % s(6): CH-umbrella -! % s(7): CH-umbrella**2 -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -module ctrans_mod - use accuracy_constants, only: dp, idp - implicit none -! precalculate pi, 2*pi and angle to radian conversion - real(dp), parameter :: pi = 4.0_dp*datan(1.0_dp) - real(dp), parameter :: pi2 = 2.0_dp*pi - real(dp), parameter :: ang2rad = pi/180.0_dp -! precalculate roots - real(dp), parameter:: sq2 = 1.0_dp/dsqrt(2.0_dp) - real(dp), parameter:: sq3 = 1.0_dp/dsqrt(3.0_dp) - real(dp), parameter:: sq6 = 1.0_dp/dsqrt(6.0_dp) -! change distances for equilibrium - !real(dp), parameter :: dchequi = 1.02289024_dp - real(dp), parameter :: dchequi = 2.344419_dp ! NO3 - !real(dp), parameter :: dchequi = 2.34451900_dp - -! see changes - -contains - subroutine ctrans(q, x1,x2, invariants) - use dim_parameter, only: qn - integer(idp) k !running indices - real(dp), intent(in) :: q(qn) !given coordinates - real(dp), intent(out) :: x1(qn) !output coordinates symmetry adapted and scaled - real(dp), intent(out) :: x2(qn) !output coordinates symmetry adapted but not scaled -! ANN Variables - real(dp), optional, intent(out) :: invariants(:) - real(dp) :: s(qn),t(qn) - ! kartesian coordianates copy from MeF+ so substitute c by n and removed f - real(dp) ch1(3), ch2(3), ch3(3), c_atom(3) - real(dp) nh1(3), nh2(3), nh3(3) - real(dp) zaxis(3), xaxis(3), yaxis(3) - real(dp) ph1(3), ph2(3), ph3(3) -! primitive coordinates - real(dp) dch1, dch2, dch3 !nh-distances - real(dp) umb !Umbrella Angle from xy-plane - -! Symmetry coordinates - real(dp) aR !a1-modes H-Dist., - real(dp) exR, exAng !ex components H-Dist., H-Ang. - real(dp) eyR, eyAng !ey components H-Dist., H-Ang. -! debugging - logical, parameter :: dbg = .false. - -! initialize coordinate vectors - s = 0.0_dp - t = 0.0_dp - -! write kartesian coords for readability - c_atom(1:3) = q(1:3) - do k = 1, 3 - ch1(k) = q(k + 3) - ch2(k) = q(k + 6) - ch3(k) = q(k + 9) - end do - -! construct z-axis - nh1 = normalized(ch1) - nh2 = normalized(ch2) - nh3 = normalized(ch3) - zaxis = create_plane(nh1, nh2, nh3) - - ! calculate bonding distance - dch1 = norm(ch1) - dch2 = norm(ch2) - dch3 = norm(ch3) - - ! construct symmertic and antisymmetric strech - aR = symmetrize(dch1 - dchequi, dch2 - dchequi, dch3 - dchequi, 'a') - exR = symmetrize(dch1, dch2, dch3, 'x') - eyR = symmetrize(dch1, dch2, dch3, 'y') - - ! construc x-axis and y axis - ph1 = normalized(project_point_into_plane(nh1, zaxis, c_atom)) - xaxis = normalized(ph1) - yaxis = xproduct(zaxis, xaxis) ! right hand side koordinates - -! project H atoms into C plane - ph2 = normalized(project_point_into_plane(nh2, zaxis, c_atom)) - ph3 = normalized(project_point_into_plane(nh3, zaxis, c_atom)) - - call construct_HBend(exAng, eyAng, ph1, ph2, ph3, xaxis, yaxis) - umb = construct_umbrella(nh1, nh2, nh3, zaxis) - -! set symmetry coordinates and even powers of umbrella - !s(1) = dch1- dchequi !aR - !s(2) = dch2 - dchequi !exR - !s(3) = dch3 - dchequi !eyR - s(1) = aR - s(2) = exR - s(3) = eyR - s(4) = exAng - s(5) = eyAng - s(6) = umb - s(7) = umb**2 - s(8) = 0 - s(9) = 0 -! pairwise distances as second coordinate set - t = 0._dp - call pair_distance(q, t(1:6)) - - if (dbg) write (6, '("sym coords s=",9f16.8)') s(1:qn) - if (dbg) write (6, '("sym coords t=",9f16.8)') t(1:qn) - if (present(invariants)) then - call get_invariants(s, invariants) - end if - - ! transform s and t to x1 and x2 - x1(1:qn)=s(1:qn) - x1(5)=x1(5) - ! set other x coordinate to zero other than strech - !X1(4:qn)=0.0d0 - x2(1:qn)=t(1:qn) - end subroutine ctrans - - subroutine pair_distance(q, r) - real(dp), intent(in) :: q(9) - real(dp), intent(out) :: r(6) - real(dp) :: atom(3, 4) - integer :: n, k, count - - !atom order: H1 H2 H3 N - atom(:, 1:3) = reshape(q, [3, 3]) - atom(:, 4) = (/0.0_dp, 0.0_dp, 0.0_dp/) - - ! disntace order 12 13 14 23 24 34 - count = 0 - do n = 1, size(atom, 2) - do k = n + 1, size(atom, 2) - count = count + 1 - r(count) = sqrt(sum((atom(:, k) - atom(:, n))**2)) - end do - end do - end subroutine pair_distance - - function morse_and_symmetrize(x,p,pst) result(s) - real(dp), intent(in),dimension(3) :: x - real(dp), intent(in),dimension(11) :: p - integer, intent(in),dimension(2) :: pst - integer :: k - real(dp), dimension(3) :: s - real(dp), dimension(3) :: t - - ! Morse transform - do k=1,3 - t(k) = morse_transform(x(k), p, pst) - end do - s(1) = symmetrize(t(1), t(2), t(3), 'a') - s(2) = symmetrize(t(1), t(2), t(3), 'x') - s(3) = symmetrize(t(1), t(2), t(3), 'y') - end function morse_and_symmetrize - - subroutine get_invariants(s, inv_out) - use dim_parameter, only: qn - use select_monom_mod, only: v_e_monom, v_ee_monom - real(dp), intent(in) :: s(qn) - real(dp), intent(out) :: inv_out(:) - ! real(dp), parameter :: ck = 1.0_dp, dk = 1.0_dp/ck ! scaling for higher order invariants - real(dp) inv(24) - integer, parameter :: inv_order(12) = & ! the order in which the invariants are selected - & [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] - real(dp) Rch, umb, xR, yR, xAng, yAng -! for readability - Rch = s(1) - xR = s(2) - yR = s(3) - xAng = s(4) - yAng = s(5) - umb = s(6)**2 -! invarianten -! a moden - inv(1) = Rch - inv(2) = umb -! invariante e pairs - inv(3) = v_e_monom(xR, yR, 1) - inv(4) = v_e_monom(xAng, yAng, 1) -! third order e pairs - inv(5) = v_e_monom(xR, yR, 2) - inv(6) = v_e_monom(xAng, yAng, 2) - ! invariant ee coupling - inv(7) = v_ee_monom(xR, yR, xAng, yAng, 1) - ! mode combinations - inv(8) = Rch*umb - - inv(9) = Rch*v_e_monom(xR, yR, 1) - inv(10) = umb*v_e_monom(xR, yR, 1) - - inv(11) = Rch*v_e_monom(xAng, yAng, 1) - inv(12) = umb*v_e_monom(xAng, yAng, 1) - -! damp coordinates because of second order and higher invariants - inv(3) = sign(sqrt(abs(inv(3))), inv(3)) - inv(4) = sign(sqrt(abs(inv(4))), inv(4)) - inv(5) = sign((abs(inv(5))**(1./3.)), inv(5)) - inv(6) = sign((abs(inv(6))**(1./3.)), inv(6)) - inv(7) = sign((abs(inv(7))**(1./3.)), inv(7)) - inv(8) = sign(sqrt(abs(inv(8))), inv(8)) - inv(9) = sign((abs(inv(9))**(1./3.)), inv(9)) - inv(10) = sign((abs(inv(10))**(1./3.)), inv(10)) - inv(11) = sign((abs(inv(11))**(1./3.)), inv(11)) - inv(12) = sign((abs(inv(12))**(1./3.)), inv(12)) - - inv_out(:) = inv(inv_order(1:size(inv_out, 1))) - - end subroutine get_invariants - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % real part of spherical harmonics -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! Ylm shifted to 0 for theta=0 - real(dp) function ylm(theta, phi, l, m) - implicit none - real(dp) theta, phi - integer(idp) l, m - ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.0_dp, l, m) - end function ylm -!---------------------------------------------------------- - real(dp) function plm2(x, l, n) - implicit none - real(dp) x - integer(idp) l, m, n - - real(dp) pmm, p_mp1m, pllm - integer(idp) ll - -! negative m und bereich von x abfangen - if ((l .lt. 0)& - &.or. (abs(n) .gt. abs(l))& - &.or. (abs(x) .gt. 1.)) then - write (6, '(''bad arguments in legendre'')') - stop - end if - -! fix sign of m to compute the positiv m - m = abs(n) - - pmm = (-1)**m*dsqrt(fac(2*m))*1./((2**m)*fac(m))& !compute P(m,m) not P(l,l) - &*(dsqrt(1.-x**2))**m - - if (l .eq. m) then - plm2 = pmm !P(l,m)=P(m,m) - else - p_mp1m = x*dsqrt(dble(2*m + 1))*pmm !compute P(m+1,m) - if (l .eq. m + 1) then - plm2 = p_mp1m !P(l,m)=P(m+1,m) - else - do ll = m + 2, l - pllm = x*(2*l - 1)/dsqrt(dble(l**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively - &- dsqrt(dble((l - 1)**2 - m**2))& - &/dsqrt(dble(l**2 - m**2))*pmm -! schreibe m+2 und m+1 jeweils fuer die naechste iteration - pmm = p_mp1m !P(m,m) = P(m+1,m) - p_mp1m = pllm !P(m+1,m) = P(m+2,m) - end do - plm2 = pllm !P(l,m)=P(m+k,m), k element N - end if - end if - -! sets the phase of -m term right (ignored to gurantee Ylm=(Yl-m)* for JT terms -! if(n.lt.0) then -! plm2 = (-1)**m * plm2 !* fac(l-m)/fac(l+m) -! endif - - end function -!---------------------------------------------------------------------------------------------------- - real(dp) function fac(i) - integer(idp) i - select case (i) - case (0) - fac = 1.0_dp - case (1) - fac = 1.0_dp - case (2) - fac = 2.0_dp - case (3) - fac = 6.0_dp - case (4) - fac = 24.0_dp - case (5) - fac = 120.0_dp - case (6) - fac = 720.0_dp - case (7) - fac = 5040.0_dp - case (8) - fac = 40320.0_dp - case (9) - fac = 362880.0_dp - case (10) - fac = 3628800.0_dp - case (11) - fac = 39916800.0_dp - case (12) - fac = 479001600.0_dp - case default - write (*, *) 'ERROR: no case for given faculty, Max is 12!' - stop - end select - end function fac - - ! Does the simplest morse transform possible - ! one skaling factor + shift - function morse_transform(x, p, pst) result(t) - real(dp), intent(in) :: x - real(dp), intent(in) :: p(11) - integer, intent(in) :: pst(2) - real(dp) :: t - if (pst(2) == 11) then - t = 1.0_dp - exp(-abs(p(2))*(x - p(1))) - else - error stop 'in morse_transform key required or wrong number of parameters' - end if - end function morse_transform - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % FUNCTION F(...) ! MAIK DEPRICATING OVER THE TOP MORSE FUNCTION FOR MYSELF -! % -! % Returns exponent of tunable Morse coordinate -! % exponent is polynomial * gaussian (skewed) -! % ilabel = 1 or 2 selects the parameters a and sfac to be used -! % -! % Background: better representation of the prefector in the -! % exponend of the morse function. -! % Formular: f(r) = lest no3 paper -! % -! % Variables: -! % x: distance of atoms (double) -! % p: parameter vector (double[20]) -! % ii: 1 for CCl and 2 for CCH (int) -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function f(x, p, ii) - integer(idp), intent(in) :: ii !1 for CCL and 2 for CCH - real(dp), intent(in) :: x !coordinate - real(dp), intent(in) :: p(11) !parameter-vector - - integer(idp) i !running index - - real(dp) r !equilibrium distance - real(dp) gaus !gaus part of f - real(dp) poly !polynom part of f - real(dp) skew !tanh part of f - - real(dp) f !prefactor of exponent and returned value - - integer(idp) npoly(2) !order of polynom - -! Maximum polynom order - npoly(1) = 5 - npoly(2) = 5 - -! p(1): position of equilibrium -! p(2): constant of exponent -! p(3): constant for skewing the gaussian -! p(4): tuning for skewing the gaussian -! p(5): Gaussian exponent -! p(6): Shift of Gaussian maximum -! p(7)...: polynomial coefficients -! p(8+n)...: coefficients of Morse Power series - -! 1-exp{[p(2)+exp{-p(5)[x-p(6)]^2}[Taylor{p(7+n)}(x-p(6))]][x-p(1)]} - -! Tunable Morse function -! Power series in Tunable Morse coordinates of order m -! exponent is polynomial of order npoly * gaussian + switching function - -! set r r-r_e - r = x - r = r - p(1) - -! set up skewing function: - skew = 0.5_dp*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.0_dp) - -! set up gaussian function: - gaus = dexp(-dabs(p(5))*(r - p(6))**2) - -! set up power series: - poly = 0.0_dp - do i = 0, npoly(ii) - 1 - poly = poly + p(7 + i)*(r - p(6))**i - end do -! set up full exponent function: - f = dabs(p(2)) + skew + gaus*poly - - end function -!---------------------------------------------------------------------------------------------------- - pure function xproduct(a, b) result(axb) - real(dp), intent(in) :: a(3), b(3) - real(dp) :: axb(3) !crossproduct a x b - axb(1) = a(2)*b(3) - a(3)*b(2) - axb(2) = a(3)*b(1) - a(1)*b(3) - axb(3) = a(1)*b(2) - a(2)*b(1) - end function xproduct - - pure function normalized(v) result(r) - real(dp), intent(in) :: v(:) - real(dp) :: r(size(v)) - r = v/norm(v) - end function normalized - - pure function norm(v) result(n) - real(dp), intent(in) :: v(:) - real(dp) n - n = dsqrt(sum(v(:)**2)) - end function norm - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % FUNCTION Project_Point_Into_Plane(x,n,r0) result(p) -! % return the to n orthogonal part of a vector x-r0 -! % p: projected point in plane -! % x: point being projected -! % n: normalvector of plane -! % r0: Point in plane -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function project_point_into_plane(x, n, r0) result(p) - real(dp), intent(in) :: x(:), n(:), r0(:) - real(dp) :: p(size(x)), xs(size(x)) - xs = x - r0 - p = xs - plane_to_point(x, n, r0) - end function project_point_into_plane - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % Function Plane_To_Point(x,n,r0) result(p) -! % p: part of n in x -! % x: point being projected -! % n: normalvector of plane -! % r0: Point in plane -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function plane_to_point(x, n, r0) result(p) - real(dp), intent(in) :: x(:), n(:), r0(:) - real(dp) p(size(x)), xs(size(x)), nn(size(n)) - nn = normalized(n) - xs = x - r0 - p = dot_product(nn, xs)*nn - end function plane_to_point - - subroutine check_coordinates(q) -! check for faulty kartesain coordinates - real(dp), intent(in) :: q(:) - integer(idp) :: i - if (all(abs(q) <= epsilon(0.0_dp))) then - stop 'Error (ctrans): all kartesian coordinates are<=1d-8' - end if - do i = 1, 9, 3 - if (all(abs(q(i:i + 2)) <= epsilon(0.0_dp))) then - write (*, *) q - stop 'Error(ctrans):kartesian coordinates zero for one atom' - end if - end do - end subroutine - - pure function rotor_a_to_z(a, z) result(r) - real(dp), intent(in) :: a(3), z(3) - real(dp) :: r(3, 3) - real(dp) :: alpha - real(dp) :: s1(3), s(3, 3), rotor(3, 3) - s1 = xproduct(normalized(a), normalized(z)) - alpha = asin(norm(s1)) - s(:, 1) = normalized(s1) - s(:, 2) = normalized(z) - s(:, 3) = xproduct(s1, z) - rotor = init_rotor(alpha, 0.0_dp, 0.0_dp) - r = matmul(s, matmul(rotor, transpose(s))) - end function - -! function returning Rz(gamma) * Ry(beta) * Rx(alpha) for basis order xyz - pure function init_rotor(alpha, beta, gamma) result(rotor) - real(dp), intent(in) :: alpha, beta, gamma - real(dp) :: rotor(3, 3) - rotor = 0.0_dp - rotor(1, 1) = dcos(beta)*dcos(gamma) - rotor(1, 2) = dsin(alpha)*dsin(beta)*dcos(gamma)& - &- dcos(alpha)*dsin(gamma) - rotor(1, 3) = dcos(alpha)*dsin(beta)*dcos(gamma)& - &+ dsin(alpha)*dsin(gamma) - - rotor(2, 1) = dcos(beta)*dsin(gamma) - rotor(2, 2) = dsin(alpha)*dsin(beta)*dsin(gamma)& - &+ dcos(alpha)*dcos(gamma) - rotor(2, 3) = dcos(alpha)*dsin(beta)*dsin(gamma)& - &- dsin(alpha)*dcos(gamma) - - rotor(3, 1) = -dsin(beta) - rotor(3, 2) = dsin(alpha)*dcos(beta) - rotor(3, 3) = dcos(alpha)*dcos(beta) - end function init_rotor - - pure function create_plane(a, b, c) result(n) - real(dp), intent(in) :: a(3), b(3), c(3) - real(dp) :: n(3) - real(dp) :: axb(3), bxc(3), cxa(3) - axb = xproduct(a, b) - bxc = xproduct(b, c) - cxa = xproduct(c, a) - n = normalized(axb + bxc + cxa) - end function create_plane - - function symmetrize(q1, q2, q3, sym) result(s) - real(dp), intent(in) :: q1, q2, q3 - character, intent(in) :: sym - real(dp) :: s - select case (sym) - case ('a') - s = (q1 + q2 + q3)*sq3 - case ('x') - s = sq6*(2.0_dp*q1 - q2 - q3) - case ('y') - s = sq2*(q2 - q3) - case default - write (*, *) 'ERROR: no rule for symmetrize with sym=', sym - stop - end select - end function symmetrize - - subroutine construct_HBend(ex, ey, ph1, ph2, ph3, x_axis, y_axis) - real(dp), intent(in) :: ph1(3), ph2(3), ph3(3) - real(dp), intent(in) :: x_axis(3), y_axis(3) - real(dp), intent(out) :: ex, ey - real(dp) :: x1, y1, alpha1 - real(dp) :: x2, y2, alpha2 - real(dp) :: x3, y3, alpha3 -! get x and y components of projected points - x1 = dot_product(ph1, x_axis) - y1 = dot_product(ph1, y_axis) - x2 = dot_product(ph2, x_axis) - y2 = dot_product(ph2, y_axis) - x3 = dot_product(ph3, x_axis) - y3 = dot_product(ph3, y_axis) -! -> calculate H deformation angles - alpha3 = datan2(y2, x2) - alpha2 = -datan2(y3, x3) !-120*ang2rad -! write(*,*)' atan2' -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad - if (alpha2 .lt. 0) alpha2 = alpha2 + pi2 - if (alpha3 .lt. 0) alpha3 = alpha3 + pi2 - alpha1 = (pi2 - alpha2 - alpha3) -! write(*,*)' fixed break line' -! write(*,*) 'alpha1:' , alpha1/ang2rad -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad - alpha1 = alpha1 !- 120.0_dp*ang2rad - alpha2 = alpha2 !- 120.0_dp*ang2rad - alpha3 = alpha3 !- 120.0_dp*ang2rad -! write(*,*)' delta alpha' -! write(*,*) 'alpha1:' , alpha1/ang2rad -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad -! write(*,*) - -! construct symmetric and antisymmetric H angles - ex = symmetrize(alpha1, alpha2, alpha3, 'x') - ey = symmetrize(alpha1, alpha2, alpha3, 'y') - end subroutine construct_HBend - - pure function construct_umbrella(nh1, nh2, nh3, n)& - &result(umb) - real(dp), intent(in) :: nh1(3), nh2(3), nh3(3) - real(dp), intent(in) :: n(3) - real(dp) :: umb - real(dp) :: theta(3) -! calculate projections for umberella angle - theta(1) = dacos(dot_product(n, nh1)) - theta(2) = dacos(dot_product(n, nh2)) - theta(3) = dacos(dot_product(n, nh3)) -! construct umberella angle - umb = sum(theta(1:3))/3.0_dp - 90.0_dp*ang2rad - end function construct_umbrella - - pure subroutine construct_sphericals& - &(theta, phi, cf, xaxis, yaxis, zaxis) - real(dp), intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3) - real(dp), intent(out) :: theta, phi - real(dp) :: x, y, z, v(3) - v = normalized(cf) - x = dot_product(v, normalized(xaxis)) - y = dot_product(v, normalized(yaxis)) - z = dot_product(v, normalized(zaxis)) - theta = dacos(z) - phi = -datan2(y, x) - end subroutine construct_sphericals - - subroutine int2kart(internal, kart) - real(dp), intent(in) :: internal(6) - real(dp), intent(out) :: kart(9) - real(dp) :: h1x, h1y, h1z - real(dp) :: h2x, h2y, h2z - real(dp) :: h3x, h3y, h3z - real(dp) :: dch0, dch1, dch2, dch3 - real(dp) :: a1, a2, a3, wci - - kart = 0.0_dp - dch1 = dchequi + sq3*internal(1) + 2*sq6*internal(2) - dch2 = dchequi + sq3*internal(1) - sq6*internal(2) + sq2*internal(3) - dch3 = dchequi + sq3*internal(1) - sq6*internal(2) - sq2*internal(3) - a1 = 2*sq6*internal(4) - a2 = -sq6*internal(4) + sq2*internal(5) - a3 = -sq6*internal(4) - sq2*internal(5) - wci = internal(6) - - ! Berechnung kartesische Koordinaten - ! ----------------------- - h1x = dch1*cos(wci*ang2rad) - h1y = 0.0 - h1z = -dch1*sin(wci*ang2rad) - - h3x = dch2*cos((a2 + 120)*ang2rad)*cos(wci*ang2rad) - h3y = dch2*sin((a2 + 120)*ang2rad)*cos(wci*ang2rad) - h3z = -dch2*sin(wci*ang2rad) - - h2x = dch3*cos((-a3 - 120)*ang2rad)*cos(wci*ang2rad) - h2y = dch3*sin((-a3 - 120)*ang2rad)*cos(wci*ang2rad) - h2z = -dch3*sin(wci*ang2rad) - - kart(1) = h1x - kart(2) = h1y - kart(3) = h1z - kart(4) = h2x - kart(5) = h2y - kart(6) = h2z - kart(7) = h3x - kart(8) = h3y - kart(9) = h3z - end subroutine int2kart - -end module ctrans_mod diff --git a/src/model/Maik_chngs/select_monom_mod.f90 b/src/model/Maik_chngs/select_monom_mod.f90 deleted file mode 100644 index 71ddb9e..0000000 --- a/src/model/Maik_chngs/select_monom_mod.f90 +++ /dev/null @@ -1,2293 +0,0 @@ -module select_monom_mod - use accuracy_constants, only: dp, idp - implicit none - ! !taylor expansion prefactors - ! real(dp),parameter,private :: pre1 = 1.0_dp - ! real(dp),parameter,private :: pre2 = 0.5_dp - ! real(dp),parameter,private :: pre3 = 1.0_dp/6.0_dp - ! real(dp),parameter,private :: pre4 = 1.0_dp/24.0_dp - ! real(dp),parameter,private :: pre5 = 1.0_dp/120.0_dp - ! real(dp),parameter,private :: pre6 = 1.0_dp/720.0_dp - ! real(dp),parameter,private :: pre7 = 1.0_dp/5040.0_dp - ! real(dp),parameter,private :: pre8 = 1.0_dp/40320.0_dp - ! real(dp),parameter,private :: pre9 = 1.0_dp/362880.0_dp - ! real(dp),parameter,private :: pre10 = 1.0_dp/3628800.0_dp - - ! ! primitive factors - ! real(dp),parameter,private :: fac1 = 1.0_dp - ! real(dp),parameter,private :: fac2 = (1._dp/2._dp) - ! real(dp),parameter,private :: fac3 = (1._dp/3._dp) - ! real(dp),parameter,private :: fac4 = (1._dp/4._dp) - ! real(dp),parameter,private :: fac5 = (1._dp/5._dp) - ! real(dp),parameter,private :: fac6 = (1._dp/6._dp) - ! real(dp),parameter,private :: fac7 = (1._dp/7._dp) - ! real(dp),parameter,private :: fac8 = (1._dp/8._dp) - - ! !missing factors for mixed monomials - ! real(dp),parameter,private :: fac31 = fac3 * fac2 - ! real(dp),parameter,private :: fac32 = fac3 - - ! real(dp),parameter,private :: fac41 = fac4 * fac3 * fac2 - ! real(dp),parameter,private :: fac42 = fac4 * fac3 - ! real(dp),parameter,private :: fac43 = fac4 - - ! real(dp),parameter,private :: fac51 = fac5 * fac4 * fac3 * fac2 - ! real(dp),parameter,private :: fac52 = fac5 * fac4 * fac3 - ! real(dp),parameter,private :: fac53 = fac5 * fac4 - ! real(dp),parameter,private :: fac54 = fac5 - - ! real(dp),parameter,private :: fac61 = fac6 * fac5 * fac4 * fac3 * fac2 - ! real(dp),parameter,private :: fac62 = fac6 * fac5 * fac4 * fac3 - ! real(dp),parameter,private :: fac63 = fac6 * fac5 * fac4 - ! real(dp),parameter,private :: fac64 = fac6 * fac5 - ! real(dp),parameter,private :: fac65 = fac6 - -contains -! V Monom - function v_e_monom(x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = (x**2 + y**2) - case (2) - monom = (2.0_dp*(x**3) - 6.0_dp*(x*y**2)) - case (3) - monom = (x**4 + 2.0_dp*(x**2*y**2) + y**4) - case (4) - monom = (2.0_dp*(x**5) - 4.0_dp*(x**3*y**2) - 6.0_dp*(x*y**4)) - case (5) - monom = (2.0_dp*(x**6) - 30.0_dp*(x**4*y**2) + 30.0_dp*(x**2*y**4) - 2.0_dp*(y**6)) - case (6) - monom = (x**6 + 3.0_dp*(x**4*y**2) + 3.0_dp*(x**2*y**4) + y**6) - case (7) - monom = (2.0_dp*x**7 - 2.0_dp*x**5*y**2 - 10.0_dp*x**3*y**4 - 6.0_dp*x*y**6) - case (8) - monom = (2.0_dp*x**8 - 28.0_dp*x**6*y**2 + 28.0_dp*x**2*y**6 - 2.0_dp*y**8) - case (9) - monom = (x**8 + 4.0_dp*x**6*y**2 + 6.0_dp*x**4*y**4 + 4.0_dp*x**2*y**6 + y**8) - case default - error stop 'called case of v_e_monom not implemented' - monom = 0.0_dp - end select - end function v_e_monom - - function v_aa_monom(a, b, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b - integer(kind=idp), intent(in) :: o - select case (o) - case (1) !order 2 - monom = a*b - - case (2) !order 3 (2) - monom = a**2*b - case (3) - monom = a*b**2 - - case (4) !order 4 (3) - monom = a**3*b - case (5) - monom = a**2*b**2 - case (6) - monom = a*b**3 - - case (7) !order 5 (4) - monom = a**4*b - case (8) - monom = a**3*b**2 - case (9) - monom = a**2*b**3 - case (10) - monom = a*b**4 - - case (11) !order 6 (5) - monom = a**5*b - case (12) - monom = a**4*b**2 - case (13) - monom = a**3*b**3 - case (14) - monom = a**2*b**4 - case (15) - monom = a*b**5 - case default - error stop 'called case of v_aa_monom not implemented' - monom = 0.0_dp - end select - end function v_aa_monom - - function v_ee_monom(x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = (2.0_dp*x1*x2 + 2.0_dp*y1*y2) !order 2 - case (2) - monom = (2.0_dp*x1*x2**2 - 2.0_dp*x1*y2**2& !order 3 - &- 4.0_dp*x2*y1*y2) - case (3) - monom = (2.0_dp*x1**2*x2 - 2.0_dp*x2*y1**2& !order 3 - &- 4.0_dp*x1*y1*y2) - case (4) - monom = (2.0_dp*x1**2*x2**2 - 2.0_dp*x1**2*y2**2& !order 4 - &- 2.0_dp*x2**2*y1**2 + 2.0_dp*y1**2*y2**2& - &+ 8.0_dp*x1*x2*y1*y2) - case (5) - monom = (2.0_dp*x1**3*x2 + 2.0_dp*x1**2*y1*y2& - &+ 2.0_dp*x1*x2*y1**2 + 2.0_dp*y1**3*y2) - case (6) - monom = (2.0_dp*x2**3*x1 + 2.0_dp*x1*x2*y2**2& - &+ 2.0_dp*x2**2*y1*y2 + 2.0_dp*y2**3*y1) - case (7) - monom = (x1**2*x2**2 + x1**2*y2**2& - &+ x2**2*y1**2 + y1**2*y2**2) - case (8) !order 5 - monom = (x1*y2**4 - 4*x2*y1*y2**3 - 6*x1*x2**2*y2**2 + 4*x2**3*y1*y2& - &+ x1*x2**4) - case (9) - monom = (y2**2 + x2**2)*(x1*y2**2 + 2*x2*y1*y2 - x1*x2**2) - case (10) - monom = (y1**2 + x1**2)*(x1*y2**2 + 2*x2*y1*y2 - x1*x2**2) - case (11) - monom = x2*(y1**2 + x1**2)*(3.0_dp*y2**2 - x2**2) - case (12) - monom = (2.0_dp*x1*y1*y2 + x2*y1**2 - x1**2*x2)*(y2**2 + x2**2) - case (13) - monom = x1*(3.0_dp*y1**2 - x1**2)*(y2**2 + x2**2) - case (14) - monom = (4.0_dp*x1*y1**3*y2 - 4.0_dp*x1**3*y1*y2 - x2*y1**4& - &+ 6.0_dp*x1**2*x2*y1**2 - x1**4*x2) - case (15) - monom = (y1**2 + x1**2)*(2.0_dp*x1*y1*y2 + x2*y1**2 - x1**2*x2) - case (16) !order 6 - monom = (y1*y2**5 - 5.0_dp*x1*x2*y2**4 - 10.0_dp*x2**2*y1*y2**3& - &+ 10.0_dp*x1*x2**3*y2**2 + 5.0_dp*x2**4*y1*y2 - x1*x2**5) - case (17) - monom = (y1*y2**2 - x1*y2**2 - 2.0_dp*x2*y1*y2 - 2.0_dp*x1*x2*y2& - &- x2**2*y1 + x1*x2**2)*(y1*y2**2 + x1*y2**2 + 2.0_dp*x2*y1*y2& - &- 2.0_dp*x1*x2*y2 - x2**2*y1 - x1*x2**2) - case (18) - monom = (y1*y2 - x1*x2)*(y1**2*y2**2 - 3.0_dp*x1**2*y2**2& - &- 8.0_dp*x1*x2*y1*y2 - 3.0_dp*x2**2*y1**2 + x1**2*x2**2) - case (19) - monom = (y1**2*y2 - 2.0_dp*x1*y1*y2 - x1**2*y2 - x2*y1**2& - &- 2.0_dp*x1*x2*y1 + x1**2*x2)*(y1**2*y2 + 2.0_dp*x1*y1*y2& - &- x1**2*y2 + x2*y1**2 - 2.0_dp*x1*x2*y1 - x1**2*x2) - case (20) - monom = (y1**5*y2 - 10.0_dp*x1**2*y1**3*y2 + 5.0_dp*x1**4*y1*y2& - &- 5.0_dp*x1*x2*y1**4 + 10.0_dp*x1**3*x2*y1**2 - x1**5*x2) - case (21) - - monom = (y1*y2 + x1*x2)*(y2**2 + x2**2)**2 - case (22) - monom = (y1*y2 - x1*y2 + x2*y1 + x1*x2)*(y1*y2 + x1*y2& - &- x2*y1 + x1*x2)*(y2**2 + x2**2) - case (23) - monom = (y1**2 + x1**2)*(y2**2 + x2**2)**2 - case (24) - monom = (y1*y2 + x1*x2)*(y1**2*y2**2 - 3.0_dp*x1**2*y2**2& - &+ 8.0_dp*x1*x2*y1*y2 - 3.0_dp*x2**2*y1**2 + x1**2*x2**2) - case (25) - monom = (y1**2 + x1**2)*(y1*y2 + x1*x2)*(y2**2 + x2**2) - case (26) - monom = (y1**2 + x1**2)*(y1*y2 - x1*y2 + x2*y1 + x1*x2)& - &*(y1*y2 + x1*y2 - x2*y1 + x1*x2) - case (27) - monom = (y1**2 + x1**2)**2*(y2**2 + x2**2) - case (28) - monom = (y1**2 + x1**2)**2*(y1*y2 + x1*x2) - case default - error stop 'called case of v_ee_monom not implemented' - monom = 0.0_dp - end select - end function v_ee_monom - - function v_ae_monom(a, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*v_e_monom(x, y, 1) ! order 3 - - case (2) - monom = a**2*v_e_monom(x, y, 1) ! order 4 - case (3) - monom = a*v_e_monom(x, y, 2) ! order 4 - - case (4) - monom = a**3*v_e_monom(x, y, 1) ! order 5 - case (5) - monom = a**2*v_e_monom(x, y, 2) ! order 5 - case (6) - monom = a*v_e_monom(x, y, 3) ! order 5 - - case (7) - monom = a**4*v_e_monom(x, y, 1) ! order 6 - case (8) - monom = a**3*v_e_monom(x, y, 2) ! order 6 - case (9) - monom = a**2*v_e_monom(x, y, 3) ! order 6 - case (10) - monom = a**1*v_e_monom(x, y, 4) ! order 6 - - case default - error stop 'called case of v_ae_monom not implemented' - monom = 0.0_dp - end select - end function v_ae_monom - - function v_aee_monom(a, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*v_ee_monom(x1, y1, x2, y2, 1) ! order 3 - - case (2) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 1) ! order 4 - case (3) - monom = a*v_ee_monom(x1, y1, x2, y2, 2) ! order 4 - case (4) - monom = a*v_ee_monom(x1, y1, x2, y2, 3) ! order 4 - - case (5) - monom = a**3*v_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (6) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 2) ! order 5 - case (7) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 3) ! order 5 - case (8) - monom = a*v_ee_monom(x1, y1, x2, y2, 4) ! order 5 - case (9) - monom = a*v_ee_monom(x1, y1, x2, y2, 5) ! order 5 - case (10) - monom = a*v_ee_monom(x1, y1, x2, y2, 6) ! order 5 - case (11) - monom = a*v_ee_monom(x1, y1, x2, y2, 7) ! order 5 - - case (12) - monom = a**4*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (13) - monom = a**3*v_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case (14) - monom = a**3*v_ee_monom(x1, y1, x2, y2, 3) ! order 6 - case (15) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 4) ! order 6 - case (16) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 5) ! order 6 - case (17) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 6) ! order 6 - case (18) - monom = a**2*v_ee_monom(x1, y1, x2, y2, 7) ! order 6 - case (19) - monom = a*v_ee_monom(x1, y1, x2, y2, 8) ! order 6 - case (20) - monom = a*v_ee_monom(x1, y1, x2, y2, 9) ! order 6 - case (21) - monom = a*v_ee_monom(x1, y1, x2, y2, 10) ! order 6 - case (22) - monom = a*v_ee_monom(x1, y1, x2, y2, 11) ! order 6 - case (23) - monom = a*v_ee_monom(x1, y1, x2, y2, 12) ! order 6 - case (24) - monom = a*v_ee_monom(x1, y1, x2, y2, 13) ! order 6 - case (25) - monom = a*v_ee_monom(x1, y1, x2, y2, 14) ! order 6 - case (26) - monom = a*v_ee_monom(x1, y1, x2, y2, 15) ! order 6 - case default - error stop 'called case of v_aee_monom not implemented' - monom = 0.0_dp - end select - end function v_aee_monom - - function v_eee_monom(x1, y1, x2, y2, x3, y3, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2, x3, y3 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = (x1*x2*x3 - y1*y2*x3 - y1*x2*y3 - x1*y2*y3) - case (2) ! order 4 - monom = ((x1**2)*x2*x3 - (y1**2)*x2*x3 + 2.0_dp*x1*y1*y2*x3& - &+ 2.0_dp*x1*y1*x2*y3 - (x1**2)*y2*y3 + (y1**2)*y2*y3) - case (3) ! order 4 - monom = ((x1**2)*x2*x3 + (y1**2)*x2*x3 + (x1**2)*y2*y3& - &+ (y1**2)*y2*y3) - case (4) !order 4 - monom = ((x2**2)*x1*x3 - (y2**2)*x1*x3 + 2.0_dp*x2*y2*y1*x3& - &+ 2.0_dp*x2*y2*x1*y3 - (x2**2)*y1*y3 + (y2**2)*y1*y3) - case (5) !order 4 - monom = ((x2**2)*x1*x3 + (y2**2)*x1*x3 + (x2**2)*y1*y3& - &+ (y2**2)*y1*y3) - case (6) !order 4 - monom = ((x3**2)*x2*x1 - (y3**2)*x2*x1 + 2.0_dp*x3*y3*y2*x1& - &+ 2.0_dp*x3*y3*x2*y1 - (x3**2)*y2*y1 + (y3**2)*y2*y1) - case (7) !order 4 - monom = ((x3**2)*x2*x1 + (y3**2)*x2*x1 + (x3**2)*y2*y1& - &+ (y3**2)*y2*y1) - case default - error stop 'called case of v_e_monom not implemented' - monom = 0.0_dp - end select - end function v_eee_monom - - function v_aaa_monom(a, b, c, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, c - integer(kind=idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = a*b*c - - case (2) !order 4 (3) - monom = a**2*b*c - case (3) - monom = a*b**2*c - case (4) - monom = a*b*c**2 - - case (5) !oreder 5 (6) - monom = a**3*b*c - case (6) - monom = a*b**3*c - case (7) - monom = a*b*c**3 - case (8) - monom = a**2*b**2*c - case (9) - monom = a**2*b*c**2 - case (10) - monom = a*b**2*c**2 - - case (11) !order 6 (9) - monom = a**4*b*c - case (12) - monom = a*b**4*c - case (13) - monom = a*b*c**4 - - case (14) - monom = a**3*b**2*c - case (15) - monom = a**3*b*c**2 - case (16) - monom = a**2*b**3*c - case (17) - monom = a**2*b*c**3 - - case (18) - monom = a*b**3*c**2 - case (19) - monom = a*b**2*c**3 - case default - error stop 'called case of v_aaa_monom not implemented' - monom = 0.0_dp - end select - end function v_aaa_monom - - function v_aae_monom(a, b, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*v_e_monom(x, y, 1) ! order 4 - - case (2) - monom = v_aa_monom(a, b, 3)*v_e_monom(x, y, 1) ! order 5 - case (3) - monom = v_aa_monom(a, b, 2)*v_e_monom(x, y, 1) ! order 5 - case (4) - monom = v_aa_monom(a, b, 1)*v_e_monom(x, y, 2) ! order 5 - - case (5) - monom = v_aa_monom(a, b, 6)*v_e_monom(x, y, 1) ! order 6 - case (6) - monom = v_aa_monom(a, b, 5)*v_e_monom(x, y, 1) ! order 6 - case (7) - monom = v_aa_monom(a, b, 4)*v_e_monom(x, y, 1) ! order 6 - - case (8) - monom = v_aa_monom(a, b, 3)*v_e_monom(x, y, 2) ! order 6 - case (9) - monom = v_aa_monom(a, b, 2)*v_e_monom(x, y, 2) ! order 6 - - case (10) - monom = v_aa_monom(a, b, 1)*v_e_monom(x, y, 3) ! order 6 - - case default - error stop 'called case of v_aae_monom not implemented' - monom = 0.0_dp - end select - end function v_aae_monom - - function v_aaee_monom(a, b, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 1) ! order 4 - - case (2) - monom = v_aa_monom(a, b, 3)*v_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (3) - monom = v_aa_monom(a, b, 2)*v_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (4) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 2) ! order 5 - case (5) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 3) ! order 5 - - case (6) - monom = v_aa_monom(a, b, 6)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (7) - monom = v_aa_monom(a, b, 5)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (8) - monom = v_aa_monom(a, b, 4)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - - case (9) - monom = v_aa_monom(a, b, 3)*v_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case (10) - monom = v_aa_monom(a, b, 2)*v_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case (11) - monom = v_aa_monom(a, b, 3)*v_ee_monom(x1, y1, x2, y2, 3) ! order 6 - case (12) - monom = v_aa_monom(a, b, 2)*v_ee_monom(x1, y1, x2, y2, 3) ! order 6 - - case (13) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 4) ! order 6 - case (14) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 5) ! order 6 - case (15) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 6) ! order 6 - case (16) - monom = v_aa_monom(a, b, 1)*v_ee_monom(x1, y1, x2, y2, 7) ! order 6 - - case default - error stop 'called case of v_aaaee_monom not implemented' - monom = 0.0_dp - end select - end function v_aaee_monom - - function v_aeee_monom(a, x1, y1, x2, y2, x3, y3, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x1, y1, x2, y2, x3, y3 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 1) ! order 4 - - case (2) - monom = a**2*v_eee_monom(x1, y1, x2, y2, x3, y3, 1) ! order 5 - case (3) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 7) ! order 5 - case (4) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 6) ! order 5 - case (5) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 5) ! order 5 - case (6) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 4) ! order 5 - case (7) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 3) ! order 5 - case (8) - monom = a*v_eee_monom(x1, y1, x2, y2, x3, y3, 2) ! order 5 - - case (9) - error stop 'v_aeee_monom only fully implemented up to 5 order' -! monom = a**3 * v_eee_monom(x1,y1,x2,y2,x3,y3,1) ! order 6 -! case(8) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,7) ! order 6 -! case(9) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,6) ! order 6 -! case(10) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,5) ! order 6 -! case(11) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,4) ! order 6 -! case(12) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,3) ! order 6 -! case(13) -! monom = a**2 * v_eee_monom(x1,y1,x2,y2,x3,y3,2) ! order 6 -! -! case(14) -! monom = a**1 * v_eee_monom(x1,y1,x2,y2,x3,y3,8..) ! order 6 -! insert missing order of v_eee here - case default - error stop 'called case of v_eee_monom not implemented' - monom = 0.0_dp - end select - end function v_aeee_monom - - function v_aaae_monom(a, b, c, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, c, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aaa_monom(a, b, c, 1)*v_e_monom(x, y, 1) ! order 5 - - case (2) - monom = v_aaa_monom(a, b, c, 4)*v_e_monom(x, y, 1) ! order 6 - case (3) - monom = v_aaa_monom(a, b, c, 3)*v_e_monom(x, y, 1) ! order 6 - case (4) - monom = v_aaa_monom(a, b, c, 2)*v_e_monom(x, y, 1) ! order 6 - case (5) - monom = v_aaa_monom(a, b, c, 1)*v_e_monom(x, y, 2) ! order 6 - case default - error stop 'called case of v_aaae_monom not implemented' - monom = 0.0_dp - end select - end function v_aaae_monom - - function v_aaaee_monom(a, b, c, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, c, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aaa_monom(a, b, c, 1)*v_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (2) - monom = v_aaa_monom(a, b, c, 1)*v_ee_monom(x1, y1, x2, y2, 1) ! order 5 - - case (3) - monom = v_aaa_monom(a, b, c, 4)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (4) - monom = v_aaa_monom(a, b, c, 3)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (5) - monom = v_aaa_monom(a, b, c, 2)*v_ee_monom(x1, y1, x2, y2, 1) ! order 6 - - case (6) - monom = v_aaa_monom(a, b, c, 1)*v_ee_monom(x1, y1, x2, y2, 3) ! order 6 - case (7) - monom = v_aaa_monom(a, b, c, 1)*v_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case default - error stop 'called case of v_aaaee_monom not implemented' - monom = 0.0_dp - end select - end function v_aaaee_monom - -! WZ Monom - function w_e_monom(x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = x - case (2) - monom = (x**2 - y**2) - case (3) - monom = (x**3 + x*y**2) - case (4) - monom = (x**4 - 6.0_dp*(x**2*y**2) + (y**4)) - case (5) - monom = (x**4 - y**4) - case (6) - monom = (x**5 - 10.0_dp*(x**3*y**2) + 5.0_dp*(x*y**4)) - case (7) - monom = (x**5 + 2.0_dp*(x**3*y**2) + (x*y**4)) - case (8) - monom = (x**6 - 5.0_dp*(x**4*y**2)& - &- 5.0_dp*(x**2*y**4) + y**6) - case (9) - monom = (x**6 + (x**4*y**2)& - &- (x**2*y**4) - y**6) - case default - error stop 'called case of w_e_monom not implemented' - monom = 0.0_dp - end select - end function w_e_monom - - function z_e_monom(x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = y - case (2) - monom = -2.0_dp*(x*y) - case (3) - monom = ((x**2*y) + y**3) - case (4) - monom = (4.0_dp*(x**3*y) - 4.0_dp*(x*y**3)) - case (5) - monom = (-2.0_dp*(x**3*y) - 2.0_dp*(x*y**3)) - case (6) - monom = (-5.0_dp*(x**4*y) + 10.0_dp*(x**2*y**3) - y**5) - case (7) - monom = ((x**4*y) + 2.0_dp*(x**2*y**3) + y**5) - case (8) - monom = (4.0_dp*(x**5*y) - 4.0_dp*(x*y**5)) - case (9) - monom = (-2.0_dp*(x**5*y) - 4.0_dp*(x**3*y**3)& - &- 2.0_dp*(x*y**5)) - case default - error stop 'called case of z_e_monom not implemented' - monom = 0.0_dp - end select - end function z_e_monom - - function w_ae_monom(a, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*w_e_monom(x, y, 1) ! order 2 - - case (2) - monom = a**2*w_e_monom(x, y, 1) ! order 3 - case (3) - monom = a*w_e_monom(x, y, 2) ! order 3 - - case (4) - monom = a**3*w_e_monom(x, y, 1) ! order 4 - case (5) - monom = a**2*w_e_monom(x, y, 2) ! order 4 - case (6) - monom = a*w_e_monom(x, y, 3) ! order 4 - - case (7) - monom = a**4*w_e_monom(x, y, 1) ! order 5 - case (8) - monom = a**3*w_e_monom(x, y, 2) ! order 5 - case (9) - monom = a**2*w_e_monom(x, y, 3) ! order 5 - case (10) - monom = a*w_e_monom(x, y, 4) ! order 5 - case (11) - monom = a*w_e_monom(x, y, 5) ! order 5 - - case (12) - monom = a**5*w_e_monom(x, y, 1) ! order 6 - case (13) - monom = a**4*w_e_monom(x, y, 2) ! order 6 - case (14) - monom = a**3*w_e_monom(x, y, 3) ! order 6 - case (15) - monom = a**2*w_e_monom(x, y, 4) ! order 6 - case (16) - monom = a**2*w_e_monom(x, y, 5) ! order 6 - case (17) - monom = a*w_e_monom(x, y, 6) ! order 6 - case (18) - monom = a*w_e_monom(x, y, 7) ! order 6 - - case default - error stop 'called case of w_ae_monom not implemented' - monom = 0.0_dp - end select - end function w_ae_monom - - function z_ae_monom(a, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*z_e_monom(x, y, 1) ! order 2 - - case (2) - monom = a**2*z_e_monom(x, y, 1) ! order 3 - case (3) - monom = a*z_e_monom(x, y, 2) ! order 3 - - case (4) - monom = a**3*z_e_monom(x, y, 1) ! order 4 - case (5) - monom = a**2*z_e_monom(x, y, 2) ! order 4 - case (6) - monom = a*z_e_monom(x, y, 3) ! order 4 - - case (7) - monom = a**4*z_e_monom(x, y, 1) ! order 5 - case (8) - monom = a**3*z_e_monom(x, y, 2) ! order 5 - case (9) - monom = a**2*z_e_monom(x, y, 3) ! order 5 - case (10) - monom = a*z_e_monom(x, y, 4) ! order 5 - case (11) - monom = a*z_e_monom(x, y, 5) ! order 5 - - case (12) - monom = a**5*z_e_monom(x, y, 1) ! order 6 - case (13) - monom = a**4*z_e_monom(x, y, 2) ! order 6 - case (14) - monom = a**3*z_e_monom(x, y, 3) ! order 6 - case (15) - monom = a**2*z_e_monom(x, y, 4) ! order 6 - case (16) - monom = a**2*z_e_monom(x, y, 5) ! order 6 - case (17) - monom = a*z_e_monom(x, y, 6) ! order 6 - case (18) - monom = a*z_e_monom(x, y, 7) ! order 6 - - case default - error stop 'called case of z_ae_monom not implemented' - monom = 0.0_dp - end select - end function z_ae_monom - - function w_aae_monom(a, b, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 1) ! order 2,1 => 3 - - case (2) - monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 1) ! order 3,1 => 4 - case (3) - monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 1) ! order 3,1 => 4 - case (4) - monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 2) ! order 2,2 => 4 - - case (5) - monom = v_aa_monom(a, b, 6)*w_e_monom(x, y, 1) ! order 4,1 => 5 - case (6) - monom = v_aa_monom(a, b, 5)*w_e_monom(x, y, 1) ! order 4,1 => 5 - case (7) - monom = v_aa_monom(a, b, 4)*w_e_monom(x, y, 1) ! order 4,1 => 5 - - case (8) - monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 2) ! order 3,2 => 5 - case (9) - monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 2) ! order 3,2 => 5 - case (10) - monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 2) ! order 2,3 => 5 - - case (11) - monom = v_aa_monom(a, b, 10)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (12) - monom = v_aa_monom(a, b, 9)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (13) - monom = v_aa_monom(a, b, 8)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (14) - monom = v_aa_monom(a, b, 7)*w_e_monom(x, y, 1) ! order 5,1 => 6 - - case (15) - monom = v_aa_monom(a, b, 6)*w_e_monom(x, y, 2) ! order 4,2 => 6 - case (16) - monom = v_aa_monom(a, b, 5)*w_e_monom(x, y, 2) ! order 4,2 => 6 - case (17) - monom = v_aa_monom(a, b, 4)*w_e_monom(x, y, 2) ! order 4,2 => 6 - - case (18) - monom = v_aa_monom(a, b, 3)*w_e_monom(x, y, 3) ! order 3,3 => 6 - case (19) - monom = v_aa_monom(a, b, 2)*w_e_monom(x, y, 3) ! order 3,3 => 6 - case (20) - monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 4) ! order 2,4 => 6 - case (21) - monom = v_aa_monom(a, b, 1)*w_e_monom(x, y, 5) ! order 2,4 => 6 - - case default - error stop 'called case of w_aae_monom not implemented' - monom = 0.0_dp - end select - end function w_aae_monom - - function z_aae_monom(a, b, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 1) ! order 2,1 => 3 - - case (2) - monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 1) ! order 3,1 => 4 - case (3) - monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 1) ! order 3,1 => 4 - case (4) - monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 2) ! order 2,2 => 4 - - case (5) - monom = v_aa_monom(a, b, 6)*z_e_monom(x, y, 1) ! order 4,1 => 5 - case (6) - monom = v_aa_monom(a, b, 5)*z_e_monom(x, y, 1) ! order 4,1 => 5 - case (7) - monom = v_aa_monom(a, b, 4)*z_e_monom(x, y, 1) ! order 4,1 => 5 - - case (8) - monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 2) ! order 3,2 => 5 - case (9) - monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 2) ! order 3,2 => 5 - case (10) - monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 2) ! order 2,3 => 5 - - case (11) - monom = v_aa_monom(a, b, 10)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (12) - monom = v_aa_monom(a, b, 9)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (13) - monom = v_aa_monom(a, b, 8)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (14) - monom = v_aa_monom(a, b, 7)*z_e_monom(x, y, 1) ! order 5,1 => 6 - - case (15) - monom = v_aa_monom(a, b, 6)*z_e_monom(x, y, 2) ! order 4,2 => 6 - case (16) - monom = v_aa_monom(a, b, 5)*z_e_monom(x, y, 2) ! order 4,2 => 6 - case (17) - monom = v_aa_monom(a, b, 4)*z_e_monom(x, y, 2) ! order 4,2 => 6 - - case (18) - monom = v_aa_monom(a, b, 3)*z_e_monom(x, y, 3) ! order 3,3 => 6 - case (19) - monom = v_aa_monom(a, b, 2)*z_e_monom(x, y, 3) ! order 3,3 => 6 - case (20) - monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 4) ! order 2,4 => 6 - case (21) - monom = v_aa_monom(a, b, 1)*z_e_monom(x, y, 5) ! order 2,4 => 6 - - case default - error stop 'called case of z_aae_monom not implemented' - monom = 0.0_dp - end select - end function z_aae_monom - - function w_ee_monom(x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) ! order 2 - monom = x1*x2 - y1*y2 - - case (2) ! order 3 - monom = x1**2*x2 + x2*y1**2 - case (3) - monom = x1*y2**2 + x1*x2**2 - case (4) - monom = 2.0_dp*x1*y1*y2 + x1**2*x2 - x2*y1**2 - case (5) - monom = x1*x2**2 + 2.0_dp*x2*y1*y2 - x1*y2**2 - - case (6) ! order 4 - monom = x1**3*x2 - 3.0_dp*x1**2*y1*y2& - &- 3.0_dp*x1*x2*y1**2 + y1**3*y2 - case (7) - monom = x1**2*x2**2 - x1**2*y2**2 - x2**2*y1**2 + y1**2*y2**2& - &- 4.0_dp*x1*x2*y1*y2 - case (8) - monom = x2**3*x1 - 3.0_dp*x1*x2*y2**2& - &- 3.0_dp*x2**2*y1*y2 + y2**3*y1 - case (9) - monom = x1**3*x2 + 3.0_dp*x1**2*y1*y2& - &- 3.0_dp*x1*x2*y1**2 - y1**3*y2 - case (10) - monom = x2**3*x1 - 3.0_dp*x1*x2*y2**2& - &+ 3.0_dp*x2**2*y1*y2 - y2**3*y1 - case (11) - monom = x1**3*x2 - x1**2*y1*y2 + x1*x2*y1**2 - y1**3*y2 - case (12) - monom = x1**2*x2**2 + x1**2*y2**2 - x2**2*y1**2& - &- y1**2*y2**2 - case (13) - monom = x1**2*x2**2 - x1**2*y2**2 + x2**2*y1**2& - &- y1**2*y2**2 - case (14) - monom = x1*x2**3 + x1*x2*y2**2 - x2**2*y1*y2 - y2**3*y1 - - case (15) ! order 5 - monom = (x1*y2**4 + 4.0_dp*x2*y1*y2**3 - 6.0_dp*x1*x2**2*y2**2& - &- 4.0_dp*x2**3*y1*y2 + x1*x2**4) - case (16) - monom = (2.0_dp*x1*y1*y2**3 + 3.0_dp*x2*y1**2*y2**2& - &- 3.0_dp*x1**2*x2*y2**2 - 6.0_dp*x1*x2**2*y1*y2& - &- x2**3*y1**2 + x1**2*x2**3) - case (17) - monom = (3.0_dp*x1*y1**2*y2**2 - x1**3*y2**2& - &+ 2.0_dp*x2*y1**3*y2& - &- 6.0_dp*x1**2*x2*y1*y2 - 3.0_dp*x1*x2**2*y1**2& - &+ x1**3*x2**2) - case (18) - monom = (4.0_dp*x1*y1**3*y2 - 4.0_dp*x1**3*y1*y2 + x2*y1**4& - &- 6.0_dp*x1**2*x2*y1**2 + x1**4*x2) - case (19) - monom = -(y2**2 + x2**2)*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2) - case (20) - monom = x1*(y2**2 + x2**2)**2 - case (21) - monom = -(2.0_dp*x1*y1*y2**3 - 3.0_dp*x2*y1**2*y2**2& - &+ 3.0_dp*x1**2& - &*x2*y2**2 - 6.0_dp*x1*x2**2*y1*y2& - &+ x2**3*y1**2 - x1**2*x2**3) - case (22) - monom = x2*(y1**2 + x1**2)*(y2**2 + x2**2) - case (23) - monom = -(y1**2 + x1**2)*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2) - case (24) - monom = x1*(y1**2 + x1**2)*(y2**2 + x2**2) - case (25) - monom = x2*(y1**2 + x1**2)**2 - case (26) - monom = (y1**2 + x1**2)*(2.0_dp*x1*y1*y2 - x2*y1**2 + x1**2*x2) - - case (27) ! order 6 - monom = (y1*y2**5 + 5.0_dp*x1*x2*y2**4 - 10.0_dp*x2**2*y1*y2**3& - &- 10.0_dp*x1*x2**3*y2**2 + 5.0_dp*x2**4*y1*y2 + x1*x2**5) - case (28) - monom = (y2**2 + x2**2)*(y1*y2**3 - 3.0_dp*x1*x2*y2**2& - &- 3.0_dp*x2**2*y1*y2 + x1*x2**3) - case (29) - monom = (y1**2 + x1**2)*(y2**2 - 2.0_dp*x2*y2 - x2**2)& - &*(y2**2 + 2.0_dp*x2*y2 - x2**2) - case (30) - monom = (y1*y2 - x1*y2 - x2*y1 - x1*x2)*(y1*y2 + x1*y2& - &+ x2*y1 - x1*x2)*(y2**2 + x2**2) - case (31) - monom = (y1**2 + x1**2)*(y1*y2**3 - 3.0_dp*x1*x2*y2**2& - &- 3.0_dp*x2**2*y1*y2 + x1*x2**3) - case (32) - monom = (y1**3*y2 - 3.0_dp*x1**2*y1*y2& - &- 3.0_dp*x1*x2*y1**2 + x1**3*x2)*(y2**2 + x2**2) - case (33) - monom = (y1**2 + x1**2)*(y1*y2 - x1*y2 - x2*y1 - x1*x2)& - &*(y1*y2 + x1*y2 + x2*y1 - x1*x2) - case (34) - monom = (y1**2 - 2.0_dp*x1*y1 - x1**2)& - &*(y1**2 + 2.0_dp*x1*y1 - x1**2)& - &*(y2**2 + x2**2) - case (35) - monom = (y1**2 + x1**2)*(y1**3*y2 - 3.0_dp*x1**2*y1*y2& - &- 3.0_dp*x1*x2*y1**2 + x1**3*x2) - case (36) - monom = (y1**5*y2 - 10.0_dp*x1**2*y1**3*y2 + 5.0_dp*x1**4*y1*y2& - &+ 5.0_dp*x1*x2*y1**4 - 10.0_dp*x1**3*x2*y1**2 + x1**5*x2) - case (37) - monom = (y1*y2 - x1*x2)*(y2**2 + x2**2)**2 - case (38) - monom = (y2**2 + x2**2)*(y1*y2**3 + 3.0_dp*x1*x2*y2**2& - &- 3.0_dp*x2**2*y1*y2 - x1*x2**3) - case (39) - monom = (y1**2 + x1**2)*(y2 - x2)*(y2 + x2)*(y2**2 + x2**2) - case (40) - monom = (y1 - x1)*(y1 + x1)*(y2**2 + x2**2)**2 - case (41) - monom = (y1*y2**2 - x1*y2**2 + 2.0_dp*x2*y1*y2 + 2.0_dp*x1*x2*y2& - &- x2**2*y1 + x1*x2**2)*(y1*y2**2 + x1*y2**2 - 2.0_dp*x2*y1*y2& - &+ 2.0_dp*x1*x2*y2 - x2**2*y1 - x1*x2**2) - case (42) - monom = (y1**3*y2 - 3.0_dp*x1**2*y1*y2& - &+ 3.0_dp*x1*x2*y1**2 - x1**3*x2)*(y2**2 + x2**2) - case (43) - monom = (y1**2 + x1**2)*(y1*y2 - x1*x2)*(y2**2 + x2**2) - case (44) - monom = (y1**2 + x1**2)*(y1*y2**3 + 3.0_dp*x1*x2*y2**2& - &- 3.0_dp*x2**2*y1*y2 - x1*x2**3) - case (45) - monom = (y1**2*y2 - 2.0_dp*x1*y1*y2 - x1**2*y2 + x2*y1**2& - &+ 2.0_dp*x1*x2*y1 - x1**2*x2)& - &*(y1**2*y2 + 2.0_dp*x1*y1*y2 - x1**2*y2& - &- x2*y1**2 + 2.0_dp*x1*x2*y1 + x1**2*x2) - case (46) - monom = (y1 - x1)*(y1 + x1)*(y1**2 + x1**2)*(y2**2 + x2**2) - case (47) - monom = (y1**2 + x1**2)**2*(y2 - x2)*(y2 + x2) - case (48) - monom = (y1**2 + x1**2)*(y1**3*y2 - 3.0_dp*x1**2*y1*y2& - &+ 3.0_dp*x1*x2*y1**2 - x1**3*x2) - case (49) - monom = (y1**2 + x1**2)**2*(y1*y2 - x1*x2) - case default - error stop 'called case of w_ee_monom not implemented' - monom = 0.0_dp - end select - end function w_ee_monom - - function z_ee_monom(x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) ! order 2 - monom = -x1*y2 - x2*y1 - - case (2) ! order 3 - monom = y1**2*y2 + x1**2*y2 - case (3) - monom = y1*y2**2 + x2**2*y1 - case (4) - monom = y1**2*y2 + 2.0_dp*x1*x2*y1 - x1**2*y2 - case (5) - monom = y2**2*y1 + 2.0_dp*x1*x2*y2 - x2**2*y1 - - case (6) ! order 4 - monom = x1**3*y2 + 3.0_dp*x1**2*x2*y1& - &- 3.0_dp*x1*y1**2*y2 - x2*y1**3 - case (7) - monom = 2.0_dp*x1**2*x2*y2 + 2.0_dp*x1*x2**2*y1& - &- 2.0_dp*x1*y1*y2**2 - 2.0_dp*x2*y1**2*y2 - case (8) - monom = x2**3*y1 + 3.0_dp*x1*x2**2*y2& - &- 3.0_dp*x2*y1*y2**2 - x1*y2**3 - case (9) - monom = x1**3*y2 - 3.0_dp*x1**2*x2*y1 - 3.0_dp*x1*y1**2*y2& - &+ x2*y1**3 - case (10) - monom = x2**3*y1 - 3.0_dp*x1*x2**2*y2 - 3.0_dp*x2*y1*y2**2& - &+ x1*y2**3 - case (11) - monom = -x1**3*y2 - x1**2*x2*y1 - x1*y1**2*y2 - x2*y1**3 - case (12) - monom = -2.0_dp*x1*x2**2*y1 - 2.0_dp*x1*y1*y2**2 - case (13) - monom = -2.0_dp*x1**2*x2*y2 - 2.0_dp*x2*y1**2*y2 - case (14) - monom = -y2**3*x1 - x1*x2**2*y2 - x2*y1*y2**2 - y1*x2**3 - - case (15) ! order 5 - monom = -(y1*y2**4 - 4.0_dp*x1*x2*y2**3 - 6.0_dp*x2**2*y1*y2**2& - &+ 4.0_dp*x1*x2**3*y2 + x2**4*y1) - case (16) - monom = -(y1**2*y2**3 - x1**2*y2**3 - 6.0_dp*x1*x2*y1*y2**2& - &- 3.0_dp*x2**2*y1**2*y2 + 3.0_dp*x1**2*x2**2*y2& - &+ 2.0_dp*x1*x2**3*y1) - case (17) - monom = -(y1**3*y2**2 - 3.0_dp*x1**2*y1*y2**2& - &- 6.0_dp*x1*x2*y1**2*y2 + 2.0_dp*x1**3*x2*y2 - x2**2*y1**3& - &+ 3.0_dp*x1**2*x2**2*y1) - case (18) - monom = -(y1**4*y2 - 6.0_dp*x1**2*y1**2*y2 + x1**4*y2& - &- 4.0_dp*x1*x2*y1**3 + 4.0_dp*x1**3*x2*y1) - case (19) - monom = (y2**2 + x2**2)*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) - case (20) - monom = y1*(y2**2 + x2**2)**2 - case (21) - monom = (y1**2*y2**3 - x1**2*y2**3 + 6.0_dp*x1*x2*y1*y2**2& - &- 3.0_dp*x2**2*y1**2*y2 + 3.0_dp*x1**2*x2**2*y2& - &- 2.0_dp*x1*x2**3*y1) - case (22) - monom = (y1**2 + x1**2)*y2*(y2**2 + x2**2) - case (23) - monom = (y1**2 + x1**2)*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) - case (24) - monom = y1*(y1**2 + x1**2)*(y2**2 + x2**2) - case (25) - monom = (y1**2 + x1**2)**2*y2 - case (26) - monom = (y1**2 + x1**2)*(y1**2*y2 - x1**2*y2 + 2.0_dp*x1*x2*y1) - - case (27) ! order 6 - monom = (x1*y2**5 - 5.0_dp*x2*y1*y2**4 - 10.0_dp*x1*x2**2*y2**3& - &+ 10.0_dp*x2**3*y1*y2**2 + 5.0_dp*x1*x2**4*y2 - x2**5*y1) - case (28) - monom = -(y2**2 + x2**2)*(x1*y2**3 + 3.0_dp*x2*y1*y2**2& - &- 3.0_dp*x1*x2**2*y2 - x2**3*y1) - case (29) - monom = -4.0_dp*x2*(y1**2 + x1**2)*y2*(y2 - x2)*(y2 + x2) - case (30) - monom = -2.0_dp*(x1*y2 + x2*y1)*(y1*y2 - x1*x2)*(y2**2 + x2**2) - case (31) - monom = -(y1**2 + x1**2)*(x1*y2**3 + 3.0_dp*x2*y1*y2**2& - &- 3.0_dp*x1*x2**2*y2 - x2**3*y1) - case (32) - monom = -(3.0_dp*x1*y1**2*y2 - x1**3*y2 + x2*y1**3& - &- 3.0_dp*x1**2*x2*y1)*(y2**2 + x2**2) - case (33) - monom = -2.0_dp*(y1**2 + x1**2)*(x1*y2 + x2*y1)*(y1*y2 - x1*x2) - case (34) - monom = -4.0_dp*x1*y1*(y1 - x1)*(y1 + x1)*(y2**2 + x2**2) - case (35) - monom = -(y1**2 + x1**2)*(3.0_dp*x1*y1**2*y2 - x1**3*y2 + x2*y1**3& - &- 3.0_dp*x1**2*x2*y1) - case (36) - monom = -(5.0_dp*x1*y1**4*y2& - &- 10.0_dp*x1**3*y1**2*y2 + x1**5*y2& - &- x2*y1**5 + 10.0_dp*x1**2*x2*y1**3 - 5.0_dp*x1**4*x2*y1) - case (37) - monom = (x1*y2 + x2*y1)*(y2**2 + x2**2)**2 - case (38) - monom = -(y2**2 + x2**2)*(x1*y2**3 - 3.0_dp*x2*y1*y2**2& - &- 3.0_dp*x1*x2**2*y2 + x2**3*y1) - case (39) - monom = 2.0_dp*x1*y1*(y2**2 + x2**2)**2 - case (40) - monom = 2.0_dp*x2*(y1**2 + x1**2)*y2*(y2**2 + x2**2) - case (41) - monom = -2.0_dp*(x1*y2**2 - 2.0_dp*x2*y1*y2 - x1*x2**2)& - &*(y1*y2**2 + 2.0_dp*x1*x2*y2 - x2**2*y1) - case (42) - monom = (3.0_dp*x1*y1**2*y2 - x1**3*y2 - x2*y1**3& - &+ 3.0_dp*x1**2*x2*y1)*(y2**2 + x2**2) - case (43) - monom = (y1**2 + x1**2)*(x1*y2 + x2*y1)*(y2**2 + x2**2) - case (44) - monom = -(y1**2 + x1**2)*(x1*y2**3 - 3.0_dp*x2*y1*y2**2& - &- 3.0_dp*x1*x2**2*y2 + x2**3*y1) - case (45) - monom = 2.0_dp*(2.0_dp*x1*y1*y2 - x2*y1**2 + x1**2*x2)& - &*(y1**2*y2 - x1**2*y2 + 2.0_dp*x1*x2*y1) - case (46) - monom = 2.0_dp*x1*y1*(y1**2 + x1**2)*(y2**2 + x2**2) - case (47) - monom = 2.0_dp*x2*(y1**2 + x1**2)**2*y2 - case (48) - monom = (y1**2 + x1**2)*(3.0_dp*x1*y1**2*y2 - x1**3*y2 - x2*y1**3& - &+ 3.0_dp*x1**2*x2*y1) - case (49) - monom = (y1**2 + x1**2)**2*(x1*y2 + x2*y1) - case default - error stop 'called case of z_ee_monom not implemented' - monom = 0.0_dp - end select - end function z_ee_monom - - function w_aee_monom(a, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*w_ee_monom(x1, y1, x2, y2, 1) ! order 3 - - case (2) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 1) ! order 4 - case (3) - monom = a*w_ee_monom(x1, y1, x2, y2, 2) ! order 4 - case (4) - monom = a*w_ee_monom(x1, y1, x2, y2, 3) ! order 4 - case (5) - monom = a*w_ee_monom(x1, y1, x2, y2, 4) ! order 4 - case (6) - monom = a*w_ee_monom(x1, y1, x2, y2, 5) ! order 4 - - case (7) - monom = a**3*w_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (8) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 2) ! order 5 - case (9) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 3) ! order 5 - case (10) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 4) ! order 5 - case (11) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 5) ! order 5 - case (12) - monom = a*w_ee_monom(x1, y1, x2, y2, 6) ! order 5 - case (13) - monom = a*w_ee_monom(x1, y1, x2, y2, 7) ! order 5 - case (14) - monom = a*w_ee_monom(x1, y1, x2, y2, 8) ! order 5 - case (15) - monom = a*w_ee_monom(x1, y1, x2, y2, 9) ! order 5 - case (16) - monom = a*w_ee_monom(x1, y1, x2, y2, 10) ! order 5 - case (17) - monom = a*w_ee_monom(x1, y1, x2, y2, 11) ! order 5 - case (18) - monom = a*w_ee_monom(x1, y1, x2, y2, 12) ! order 5 - case (19) - monom = a*w_ee_monom(x1, y1, x2, y2, 13) ! order 5 - case (20) - monom = a*w_ee_monom(x1, y1, x2, y2, 14) ! order 5 - - case (21) - monom = a**4*w_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (22) - monom = a**3*w_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case (23) - monom = a**3*w_ee_monom(x1, y1, x2, y2, 3) ! order 6 - case (24) - monom = a**3*w_ee_monom(x1, y1, x2, y2, 4) ! order 6 - case (25) - monom = a**3*w_ee_monom(x1, y1, x2, y2, 5) ! order 6 - case (26) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 6) ! order 6 - case (27) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 7) ! order 6 - case (28) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 8) ! order 6 - case (29) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 9) ! order 6 - case (30) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 10) ! order 6 - case (31) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 11) ! order 6 - case (32) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 12) ! order 6 - case (33) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 13) ! order 6 - case (34) - monom = a**2*w_ee_monom(x1, y1, x2, y2, 14) ! order 6 - case (35) - monom = a*w_ee_monom(x1, y1, x2, y2, 15) ! order 6 - case (36) - monom = a*w_ee_monom(x1, y1, x2, y2, 16) ! order 6 - case (37) - monom = a*w_ee_monom(x1, y1, x2, y2, 17) ! order 6 - case (38) - monom = a*w_ee_monom(x1, y1, x2, y2, 18) ! order 6 - case (39) - monom = a*w_ee_monom(x1, y1, x2, y2, 19) ! order 6 - case (40) - monom = a*w_ee_monom(x1, y1, x2, y2, 20) ! order 6 - case (41) - monom = a*w_ee_monom(x1, y1, x2, y2, 21) ! order 6 - case (42) - monom = a*w_ee_monom(x1, y1, x2, y2, 22) ! order 6 - case (43) - monom = a*w_ee_monom(x1, y1, x2, y2, 23) ! order 6 - case (44) - monom = a*w_ee_monom(x1, y1, x2, y2, 24) ! order 6 - case (45) - monom = a*w_ee_monom(x1, y1, x2, y2, 25) ! order 6 - case (46) - monom = a*w_ee_monom(x1, y1, x2, y2, 26) ! order 6 - - case default - error stop 'called case of w_aee_monom not implemented' - monom = 0.0_dp - end select - end function w_aee_monom - - function z_aee_monom(a, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = a*z_ee_monom(x1, y1, x2, y2, 1) ! order 3 - - case (2) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 1) ! order 4 - case (3) - monom = a*z_ee_monom(x1, y1, x2, y2, 2) ! order 4 - case (4) - monom = a*z_ee_monom(x1, y1, x2, y2, 3) ! order 4 - case (5) - monom = a*z_ee_monom(x1, y1, x2, y2, 4) ! order 4 - case (6) - monom = a*z_ee_monom(x1, y1, x2, y2, 5) ! order 4 - - case (7) - monom = a**3*z_ee_monom(x1, y1, x2, y2, 1) ! order 5 - case (8) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 2) ! order 5 - case (9) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 3) ! order 5 - case (10) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 4) ! order 5 - case (11) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 5) ! order 5 - case (12) - monom = a*z_ee_monom(x1, y1, x2, y2, 6) ! order 5 - case (13) - monom = a*z_ee_monom(x1, y1, x2, y2, 7) ! order 5 - case (14) - monom = a*z_ee_monom(x1, y1, x2, y2, 8) ! order 5 - case (15) - monom = a*z_ee_monom(x1, y1, x2, y2, 9) ! order 5 - case (16) - monom = a*z_ee_monom(x1, y1, x2, y2, 10) ! order 5 - case (17) - monom = a*z_ee_monom(x1, y1, x2, y2, 11) ! order 5 - case (18) - monom = a*z_ee_monom(x1, y1, x2, y2, 12) ! order 5 - case (19) - monom = a*z_ee_monom(x1, y1, x2, y2, 13) ! order 5 - case (20) - monom = a*z_ee_monom(x1, y1, x2, y2, 14) ! order 5 - - case (21) - monom = a**4*z_ee_monom(x1, y1, x2, y2, 1) ! order 6 - case (22) - monom = a**3*z_ee_monom(x1, y1, x2, y2, 2) ! order 6 - case (23) - monom = a**3*z_ee_monom(x1, y1, x2, y2, 3) ! order 6 - case (24) - monom = a**3*z_ee_monom(x1, y1, x2, y2, 4) ! order 6 - case (25) - monom = a**3*z_ee_monom(x1, y1, x2, y2, 5) ! order 6 - case (26) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 6) ! order 6 - case (27) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 7) ! order 6 - case (28) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 8) ! order 6 - case (29) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 9) ! order 6 - case (30) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 10) ! order 6 - case (31) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 11) ! order 6 - case (32) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 12) ! order 6 - case (33) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 13) ! order 6 - case (34) - monom = a**2*z_ee_monom(x1, y1, x2, y2, 14) ! order 6 - case (35) - monom = a*z_ee_monom(x1, y1, x2, y2, 15) ! order 6 - case (36) - monom = a*z_ee_monom(x1, y1, x2, y2, 16) ! order 6 - case (37) - monom = a*z_ee_monom(x1, y1, x2, y2, 17) ! order 6 - case (38) - monom = a*z_ee_monom(x1, y1, x2, y2, 18) ! order 6 - case (39) - monom = a*z_ee_monom(x1, y1, x2, y2, 19) ! order 6 - case (40) - monom = a*z_ee_monom(x1, y1, x2, y2, 20) ! order 6 - case (41) - monom = a*z_ee_monom(x1, y1, x2, y2, 21) ! order 6 - case (42) - monom = a*z_ee_monom(x1, y1, x2, y2, 22) ! order 6 - case (43) - monom = a*z_ee_monom(x1, y1, x2, y2, 23) ! order 6 - case (44) - monom = a*z_ee_monom(x1, y1, x2, y2, 24) ! order 6 - case (45) - monom = a*z_ee_monom(x1, y1, x2, y2, 25) ! order 6 - case (46) - monom = a*z_ee_monom(x1, y1, x2, y2, 26) ! order 6 - - case default - error stop 'called case of z_aee_monom not implemented' - monom = 0.0_dp - end select - end function z_aee_monom - - function w_aaae_monom(a, b, c, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, c, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 1) ! order 3,1 => 4 - - case (2) - monom = v_aaa_monom(a, b, c, 4)*w_e_monom(x, y, 1) ! order 4,1 => 5 - case (3) - monom = v_aaa_monom(a, b, c, 3)*w_e_monom(x, y, 1) ! order 4,1 => 5 - case (4) - monom = v_aaa_monom(a, b, c, 2)*w_e_monom(x, y, 1) ! order 4,1 => 5 - case (5) - monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 2) ! order 3,2 => 5 - - case (6) - monom = v_aaa_monom(a, b, c, 10)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (7) - monom = v_aaa_monom(a, b, c, 9)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (8) - monom = v_aaa_monom(a, b, c, 8)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (9) - monom = v_aaa_monom(a, b, c, 7)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (10) - monom = v_aaa_monom(a, b, c, 6)*w_e_monom(x, y, 1) ! order 5,1 => 6 - case (11) - monom = v_aaa_monom(a, b, c, 5)*w_e_monom(x, y, 1) ! order 5,1 => 6 - - case (12) - monom = v_aaa_monom(a, b, c, 4)*w_e_monom(x, y, 2) ! order 4,2 => 6 - case (13) - monom = v_aaa_monom(a, b, c, 3)*w_e_monom(x, y, 2) ! order 4,2 => 6 - case (14) - monom = v_aaa_monom(a, b, c, 2)*w_e_monom(x, y, 2) ! order 4,2 => 6 - case (15) - monom = v_aaa_monom(a, b, c, 1)*w_e_monom(x, y, 3) ! order 3,3 => 6 - - case default - error stop 'called case of w_aaae_monom not implemented' - monom = 0.0_dp - end select - end function w_aaae_monom - - function z_aaae_monom(a, b, c, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, c, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 1) ! order 3,1 => 4 - - case (2) - monom = v_aaa_monom(a, b, c, 4)*z_e_monom(x, y, 1) ! order 4,1 => 5 - case (3) - monom = v_aaa_monom(a, b, c, 3)*z_e_monom(x, y, 1) ! order 4,1 => 5 - case (4) - monom = v_aaa_monom(a, b, c, 2)*z_e_monom(x, y, 1) ! order 4,1 => 5 - case (5) - monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 2) ! order 3,2 => 5 - - case (6) - monom = v_aaa_monom(a, b, c, 10)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (7) - monom = v_aaa_monom(a, b, c, 9)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (8) - monom = v_aaa_monom(a, b, c, 8)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (9) - monom = v_aaa_monom(a, b, c, 7)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (10) - monom = v_aaa_monom(a, b, c, 6)*z_e_monom(x, y, 1) ! order 5,1 => 6 - case (11) - monom = v_aaa_monom(a, b, c, 5)*z_e_monom(x, y, 1) ! order 5,1 => 6 - - case (12) - monom = v_aaa_monom(a, b, c, 4)*z_e_monom(x, y, 2) ! order 4,2 => 6 - case (13) - monom = v_aaa_monom(a, b, c, 3)*z_e_monom(x, y, 2) ! order 4,2 => 6 - case (14) - monom = v_aaa_monom(a, b, c, 2)*z_e_monom(x, y, 2) ! order 4,2 => 6 - case (15) - monom = v_aaa_monom(a, b, c, 1)*z_e_monom(x, y, 3) ! order 3,3 => 6 - - case default - error stop 'called case of z_aaae_monom not implemented' - monom = 0.0_dp - end select - end function z_aaae_monom - - function w_eee_monom(x1, y1, x2, y2, x3, y3, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2, x3, y3 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = x1*x2*x3 + y1*y2*x3 + y1*x2*y3 - x1*y2*y3 - case (2) - monom = x1*x2*x3 + y1*y2*x3 - y1*x2*y3 + x1*y2*y3 - case (3) - monom = x1*x2*x3 - y1*y2*x3 + y1*x2*y3 + x1*y2*y3 - - case (4) !order 4 - monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 - 2.d0*x1*y1*y2*x3& - &- 2.0_dp*x1*y1*x2*y3 - (x1**2)*y2*y3 + (y1**2)*y2*y3 - case (5) - monom = (x1**2)*x2*x3 + (y1**2)*x2*x3 - (x1**2)*y2*y3& - &- (y1**2)*y2*y3 - case (6) - monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 + 2.0_dp*x1*y1*y2*x3& - &- 2.0_dp*x1*y1*x2*y3 + (x1**2)*y2*y3 - (y1**2)*y2*y3 - case (7) - monom = (x1**2)*x2*x3 - (y1**2)*x2*x3 - 2.0_dp*x1*y1*y2*x3& - &+ 2.0_dp*x1*y1*x2*y3 + (x1**2)*y2*y3 - (y1**2)*y2*y3 - case (8) - - monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 - 2.0_dp*x2*y2*y1*x3& - &- 2.0_dp*x2*y2*x1*y3 - (x2**2)*y1*y3 + (y2**2)*y1*y3 - case (9) - monom = (x2**2)*x1*x3 + (y2**2)*x1*x3 - (x2**2)*y1*y3& - &- (y2**2)*y1*y3 - case (10) - monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 + 2.0_dp*x2*y2*y1*x3& - &- 2.0_dp*x2*y2*x1*y3 + (x2**2)*y1*y3 - (y2**2)*y1*y3 - case (11) - monom = (x2**2)*x1*x3 - (y2**2)*x1*x3 - 2.0_dp*x2*y2*y1*x3& - &+ 2.0_dp*x2*y2*x1*y3 + (x2**2)*y1*y3 - (y2**2)*y1*y3 - case (12) - monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 - 2.0_dp*x3*y3*y2*x1& - &- 2.0_dp*x3*y3*x2*y1 - (x3**2)*y2*y1 + (y3**2)*y2*y1 - case (13) - monom = (x3**2)*x2*x1 + (y3**2)*x2*x1 - (x3**2)*y2*y1& - &- (y3**2)*y2*y1 - case (14) - monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 + 2.0_dp*x3*y3*y2*x1& - &- 2.0_dp*x3*y3*x2*y1 + (x3**2)*y2*y1 - (y3**2)*y2*y1 - case (15) - monom = (x3**2)*x2*x1 - (y3**2)*x2*x1 - 2.0_dp*x3*y3*y2*x1& - &+ 2.0_dp*x3*y3*x2*y1 + (x3**2)*y2*y1 - (y3**2)*y2*y1 - case default - error stop 'called case of w_eee_monom not implemented' - monom = 0.0_dp - end select - end function w_eee_monom - - function z_eee_monom(x1, y1, x2, y2, x3, y3, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: x1, y1, x2, y2, x3, y3 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) ! order 3 - monom = -y1*x2*x3 + x1*y2*x3 + x1*x2*y3 + y1*y2*y3 - case (2) - monom = y1*x2*x3 - x1*y2*x3 + x1*x2*y3 + y1*y2*y3 - case (3) - monom = y1*x2*x3 + x1*y2*x3 - x1*x2*y3 + y1*y2*y3 - - case (4) ! order 4 - monom = (2.0*x1*y1*x2*x3 + (x1**2)*y2*x3 + (x1**2)*x2*y3& - &- 2.0*x1*y1*y2*y3 - (y1**2)*y2*x3 - (y1**2)*x2*y3) - case (5) - monom = (-(x1**2)*y2*x3 - (x1**2)*x2*y3 - (y1**2)*y2*x3& - &- (y1**2)*x2*y3) - case (6) - monom = -(2.0*x1*y1*x2*x3 - (x1**2)*y2*x3 + (x1**2)*x2*y3& - &+ 2.0*x1*y1*y2*y3 + (y1**2)*y2*x3 - (y1**2)*x2*y3) - case (7) - monom = -(2.0*x1*y1*x2*x3 + (x1**2)*y2*x3 - (x1**2)*x2*y3& - &+ 2.0*x1*y1*y2*y3 - (y1**2)*y2*x3 + (y1**2)*x2*y3) - - case (8) - monom = (2.d0*x2*y2*x1*x3 + (x2**2)*y1*x3 + (x2**2)*x1*y3& - &- 2.d0*x2*y2*y1*y3 - (y2**2)*y1*x3 - (y2**2)*x1*y3) - case (9) - monom = (-(x2**2)*y1*x3 - (x2**2)*x1*y3 - (y2**2)*y1*x3& - &- (y2**2)*x1*y3) - case (10) - monom = -(2.d0*x2*y2*x1*x3 - (x2**2)*y1*x3 + (x2**2)*x1*y3& - &+ 2.d0*x2*y2*y1*y3 + (y2**2)*y1*x3 - (y2**2)*x1*y3) - case (11) - monom = -(2.d0*x2*y2*x1*x3 + (x2**2)*y1*x3 - (x2**2)*x1*y3& - &+ 2.d0*x2*y2*y1*y3 - (y2**2)*y1*x3 + (y2**2)*x1*y3) - - case (12) - monom = (2.d0*x3*y3*x2*x1 + (x3**2)*y2*x1 + (x3**2)*x2*y1& - &- 2.d0*x3*y3*y2*y1 - (y3**2)*y2*x1 - (y3**2)*x2*y1) - case (13) - monom = (-(x3**2)*y2*x1 - (x3**2)*x2*y1 - (y3**2)*y2*x1& - &- (y3**2)*x2*y1) - case (14) - monom = -(2.d0*x3*y3*x2*x1 - (x3**2)*y2*x1 + (x3**2)*x2*y1& - &+ 2.d0*x3*y3*y2*y1 + (y3**2)*y2*x1 - (y3**2)*x2*y1) - case (15) - monom = -(2.d0*x3*y3*x2*x1 + (x3**2)*y2*x1 - (x3**2)*x2*y1& - &+ 2.d0*x3*y3*y2*y1 - (y3**2)*y2*x1 + (y3**2)*x2*y1) - case default - error stop 'called case of z_eee_monom not implemented' - monom = 0.0_dp - end select - end function z_eee_monom - - function w_aaee_monom(a, b, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 1) ! order 2,2 => 4 - - case (2) - monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 - case (3) - monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 - - case (4) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 2) ! order 2,3 => 5 - case (5) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 3) ! order 2,3 => 5 - case (6) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 4) ! order 2,3 => 5 - case (7) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 5) ! order 2,3 => 5 - - case (8) - monom = v_aa_monom(a, b, 6)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - case (9) - monom = v_aa_monom(a, b, 5)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - case (10) - monom = v_aa_monom(a, b, 4)*w_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - - case (11) - monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 - case (12) - monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 - case (13) - monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 - case (14) - monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 - case (15) - monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 - case (16) - monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 - case (17) - monom = v_aa_monom(a, b, 3)*w_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 - case (18) - monom = v_aa_monom(a, b, 2)*w_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 - - case (19) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 6) ! order 2,4 => 6 - case (20) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 7) ! order 2,4 => 6 - case (21) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 8) ! order 2,4 => 6 - case (22) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 9) ! order 2,4 => 6 - case (23) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 10) ! order 2,4 => 6 - case (24) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 11) ! order 2,4 => 6 - case (25) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 12) ! order 2,4 => 6 - case (26) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 13) ! order 2,4 => 6 - case (27) - monom = v_aa_monom(a, b, 1)*w_ee_monom(x1, y1, x2, y2, 14) ! order 2,4 => 6 - - case default - error stop 'called case of w_aaee_monom not implemented' - monom = 0.0_dp - end select - end function w_aaee_monom - - function z_aaee_monom(a, b, x1, y1, x2, y2, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: a, b, x1, y1, x2, y2 - integer(kind=idp), intent(in) :: o - select case (o) - case (1) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 1) ! order 2,2 => 4 - - case (2) - monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 - case (3) - monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 1) ! order 3,2 => 5 - - case (4) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 2) ! order 2,3 => 5 - case (5) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 3) ! order 2,3 => 5 - case (6) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 4) ! order 2,3 => 5 - case (7) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 5) ! order 2,3 => 5 - - case (8) - monom = v_aa_monom(a, b, 6)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - case (9) - monom = v_aa_monom(a, b, 5)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - case (10) - monom = v_aa_monom(a, b, 4)*z_ee_monom(x1, y1, x2, y2, 1) ! order 4,2 => 6 - - case (11) - monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 - case (12) - monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 2) ! order 3,3 => 6 - case (13) - monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 - case (14) - monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 3) ! order 3,3 => 6 - case (15) - monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 - case (16) - monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 4) ! order 3,3 => 6 - case (17) - monom = v_aa_monom(a, b, 3)*z_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 - case (18) - monom = v_aa_monom(a, b, 2)*z_ee_monom(x1, y1, x2, y2, 5) ! order 3,3 => 6 - - case (19) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 6) ! order 2,4 => 6 - case (20) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 7) ! order 2,4 => 6 - case (21) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 8) ! order 2,4 => 6 - case (22) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 9) ! order 2,4 => 6 - case (23) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 10) ! order 2,4 => 6 - case (24) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 11) ! order 2,4 => 6 - case (25) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 12) ! order 2,4 => 6 - case (26) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 13) ! order 2,4 => 6 - case (27) - monom = v_aa_monom(a, b, 1)*z_ee_monom(x1, y1, x2, y2, 14) ! order 2,4 => 6 - - case default - error stop 'called case of z_aaee_monom not implemented' - monom = 0.0_dp - end select - end function z_aaee_monom - - function q_u_monom(u, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: u - integer(kind=idp), intent(in) :: o - ! only odd orders of u are allowed - monom = u**(2*(o - 1) + 1) - end function q_u_monom - - function q_ua_monom(u, a, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: u, a - integer(kind=idp), intent(in) :: o - ! only odd orders of u are allowed - select case (o) - case (1) !order 2 - monom = u*a - case (2) !order 3 - monom = u*a**2 - case (3) !order 4 - monom = u**3*a - case (4) - monom = u*a**3 - case (5) ! order 5 - monom = u**3*a**2 - case (6) - monom = u*a**4 - case (7) ! order 6 - monom = u**5*a - case (8) - monom = u**3*a**3 - case (9) - monom = u*a**5 - case (10) !order 7 - monom = u**5*a**2 - case (11) - monom = u**3*a**4 - case (12) - monom = u*a**6 - case (13) ! order 8 - monom = u**7*a - case (14) - monom = u**5*a**3 - case (15) - monom = u**3*a**5 - case (16) - monom = u*a**7 - case default - error stop 'called case of q_ua_monom not implemented' - monom = 0.0_dp - end select - end function q_ua_monom - - function q_ue_monom(u, x, y, o) result(monom) - real(kind=dp) :: monom - real(kind=dp), intent(in) :: u, x, y - integer(kind=idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = u*v_e_monom(x, y, 1) - case (2) !order 4 - monom = u*v_e_monom(x, y, 2) - case (3) !order 5 - monom = u**3*v_e_monom(x, y, 1) - case (4) - monom = u*v_e_monom(x, y, 3) - case (5) !order 6 - monom = u**3*v_e_monom(x, y, 2) - case (6) - monom = u*v_e_monom(x, y, 4) - case (7) !order 7 - monom = u**5*v_e_monom(x, y, 1) - case (8) - monom = u**3*v_e_monom(x, y, 3) - case (9) - monom = u*v_e_monom(x, y, 5) - case (10) - monom = u*v_e_monom(x, y, 6) - case default - error stop 'called case of q_ue not implemented' - monom = 0.0_dp - end select - end function q_ue_monom - - function q_uae_monom(u, a, x, y, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x, y - integer(idp), intent(in) :: o - select case (o) - case (1) !order 4 - monom = u*v_ae_monom(a, x, y, 1) - case (2) !order 5 - monom = u*v_ae_monom(a, x, y, 2) - case (3) - monom = u*v_ae_monom(a, x, y, 3) - case (4) !order 6 - monom = u**3*v_ae_monom(a, x, y, 1) - case (5) - monom = u*v_ae_monom(a, x, y, 4) - case (6) - monom = u*v_ae_monom(a, x, y, 5) - case (7) - monom = u*v_ae_monom(a, x, y, 6) - case (8) !order 7 - monom = u**3*v_ae_monom(a, x, y, 2) - case (9) - monom = u**3*v_ae_monom(a, x, y, 3) - case (10) - monom = u*v_ae_monom(a, x, y, 7) - case (11) - monom = u*v_ae_monom(a, x, y, 8) - case (12) - monom = u*v_ae_monom(a, x, y, 9) - case (13) - monom = u*v_ae_monom(a, x, y, 10) - case default - error stop 'called case of q_uae not implemented' - monom = 0.0_dp - end select - end function q_uae_monom - - function q_uee_monom(u, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = u*v_ee_monom(x1, y1, x2, y2, 1) - case (2) !order 4 - monom = u*v_ee_monom(x1, y1, x2, y2, 2) - case (3) - monom = u*v_ee_monom(x1, y1, x2, y2, 3) - case (4) !order 5 - monom = u**3*v_ee_monom(x1, y1, x2, y2, 1) - case (5) - monom = u*v_ee_monom(x1, y1, x2, y2, 4) - case (6) - monom = u*v_ee_monom(x1, y1, x2, y2, 5) - case (7) - monom = u*v_ee_monom(x1, y1, x2, y2, 6) - case (8) - monom = u*v_ee_monom(x1, y1, x2, y2, 7) - case (9) ! order 6 - monom = u**3*v_ee_monom(x1, y1, x2, y2, 2) - case (10) - monom = u**3*v_ee_monom(x1, y1, x2, y2, 3) - case (11) - monom = u*v_ee_monom(x1, y1, x2, y2, 8) - case (12) - monom = u*v_ee_monom(x1, y1, x2, y2, 9) - case (13) - monom = u*v_ee_monom(x1, y1, x2, y2, 10) - case (14) - monom = u*v_ee_monom(x1, y1, x2, y2, 11) - case (15) - monom = u*v_ee_monom(x1, y1, x2, y2, 12) - case (16) - monom = u*v_ee_monom(x1, y1, x2, y2, 13) - case (17) - monom = u*v_ee_monom(x1, y1, x2, y2, 14) - case (18) - monom = u*v_ee_monom(x1, y1, x2, y2, 15) - case default - error stop 'called case of q_uee not implemented' - monom = 0.0_dp - end select - end function q_uee_monom - - function q_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 4 - monom = u*v_aee_monom(a, x1, y1, x2, y2, 1) - case (2) !order 5 - monom = u*v_aee_monom(a, x1, y1, x2, y2, 2) - case (3) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 3) - case (4) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 4) - case (5) !order 6 - monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 1) - case (6) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 5) - case (7) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 6) - case (8) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 7) - case (9) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 8) - case (10) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 9) - case (11) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 10) - case (12) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 11) - case (13) !order 7 - monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 2) - case (14) - monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 3) - case (15) - monom = u**3*v_aee_monom(a, x1, y1, x2, y2, 4) - case (16) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 12) - case (17) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 13) - case (18) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 14) - case (19) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 15) - case (20) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 16) - case (21) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 17) - case (22) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 18) - case (23) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 19) - case (24) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 20) - case (25) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 21) - case (26) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 22) - case (27) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 23) - case (28) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 24) - case (29) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 25) - case (30) - monom = u*v_aee_monom(a, x1, y1, x2, y2, 26) - case default - error stop 'called case of q_uaee not implemented' - monom = 0.0_dp - end select - end function q_uaee_monom - - function qw_ue_monom(u, x, y, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, x, y - integer(idp), intent(in) :: o - select case (o) - case (1) !order 2 - monom = u*w_e_monom(x, y, 1) - case (2) !order 3 - monom = u*w_e_monom(x, y, 2) - case (3) !order 4 - monom = u**3*w_e_monom(x, y, 1) - case (4) - monom = u*w_e_monom(x, y, 3) - case (5) !order 5 - monom = u**3*w_e_monom(x, y, 2) - case (6) - monom = u*w_e_monom(x, y, 4) - case (7) - monom = u*w_e_monom(x, y, 5) - case (8) !order 6 - monom = u**5*w_e_monom(x, y, 1) - case (9) - monom = u**3*w_e_monom(x, y, 3) - case (10) - monom = u*w_e_monom(x, y, 6) - case (11) - monom = u*w_e_monom(x, y, 7) - case (12) !order 7 - monom = u**5*w_e_monom(x, y, 2) - case (13) - monom = u**3*w_e_monom(x, y, 4) - case (14) - monom = u**3*w_e_monom(x, y, 5) - case (15) - monom = u*w_e_monom(x, y, 8) - case (16) - monom = u*w_e_monom(x, y, 9) - case default - error stop 'called case of qw_ue_monom not implemented' - monom = 0.0_dp - end select - end function qw_ue_monom - - function qz_ue_monom(u, x, y, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, x, y - integer(idp), intent(in) :: o - select case (o) - case (1) !order 2 - monom = u*z_e_monom(x, y, 1) - case (2) !order 3 - monom = u*z_e_monom(x, y, 2) - case (3) !order 4 - monom = u**3*z_e_monom(x, y, 1) - case (4) - monom = u*z_e_monom(x, y, 3) - case (5) !order 5 - monom = u**3*z_e_monom(x, y, 2) - case (6) - monom = u*z_e_monom(x, y, 4) - case (7) - monom = u*z_e_monom(x, y, 5) - case (8) !order 6 - monom = u**5*z_e_monom(x, y, 1) - case (9) - monom = u**3*z_e_monom(x, y, 3) - case (10) - monom = u*z_e_monom(x, y, 6) - case (11) - monom = u*z_e_monom(x, y, 7) - case (12) !order 7 - monom = u**5*z_e_monom(x, y, 2) - case (13) - monom = u**3*z_e_monom(x, y, 4) - case (14) - monom = u**3*z_e_monom(x, y, 5) - case (15) - monom = u*z_e_monom(x, y, 8) - case (16) - monom = u*z_e_monom(x, y, 9) - case default - error stop 'called case of qz_ue_monom not implemented' - monom = 0.0_dp - end select - end function qz_ue_monom - - function qw_uee_monom(u, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = u*w_ee_monom(x1, y1, x2, y2, 1) - case (2) !order 4 - monom = u*w_ee_monom(x1, y1, x2, y2, 2) - case (3) - monom = u*w_ee_monom(x1, y1, x2, y2, 3) - case (4) - monom = u*w_ee_monom(x1, y1, x2, y2, 4) - case (5) - monom = u*w_ee_monom(x1, y1, x2, y2, 5) - case (6) !order 5 - monom = u**3*w_ee_monom(x1, y1, x2, y2, 1) - case (7) - monom = u*w_ee_monom(x1, y1, x2, y2, 6) - case (8) - monom = u*w_ee_monom(x1, y1, x2, y2, 7) - case (9) - monom = u*w_ee_monom(x1, y1, x2, y2, 8) - case (10) - monom = u*w_ee_monom(x1, y1, x2, y2, 9) - case (11) - monom = u*w_ee_monom(x1, y1, x2, y2, 10) - case (12) - monom = u*w_ee_monom(x1, y1, x2, y2, 11) - case (13) - monom = u*w_ee_monom(x1, y1, x2, y2, 12) - case (14) - monom = u*w_ee_monom(x1, y1, x2, y2, 13) - case (15) - monom = u*w_ee_monom(x1, y1, x2, y2, 14) - case default - error stop 'called case of qw_uee_monom not implemented' - monom = 0.0_dp - end select - end function qw_uee_monom - - function qz_uee_monom(u, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 3 - monom = u*z_ee_monom(x1, y1, x2, y2, 1) - case (2) !order 4 - monom = u*z_ee_monom(x1, y1, x2, y2, 2) - case (3) - monom = u*z_ee_monom(x1, y1, x2, y2, 3) - case (4) - monom = u*z_ee_monom(x1, y1, x2, y2, 4) - case (5) - monom = u*z_ee_monom(x1, y1, x2, y2, 5) - case (6) !order 5 - monom = u**3*z_ee_monom(x1, y1, x2, y2, 1) - case (7) - monom = u*z_ee_monom(x1, y1, x2, y2, 6) - case (8) - monom = u*z_ee_monom(x1, y1, x2, y2, 7) - case (9) - monom = u*z_ee_monom(x1, y1, x2, y2, 8) - case (10) - monom = u*z_ee_monom(x1, y1, x2, y2, 9) - case (11) - monom = u*z_ee_monom(x1, y1, x2, y2, 10) - case (12) - monom = u*z_ee_monom(x1, y1, x2, y2, 11) - case (13) - monom = u*z_ee_monom(x1, y1, x2, y2, 12) - case (14) - monom = u*z_ee_monom(x1, y1, x2, y2, 13) - case (15) - monom = u*z_ee_monom(x1, y1, x2, y2, 14) - case default - error stop 'called case of qz_uee_monom not implemented' - monom = 0.0_dp - end select - end function qz_uee_monom - - function qw_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 4 - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 1) - case (2) !order 5 - monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 1) - case (3) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 2) - case (4) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 3) - case (5) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 4) - case (6) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 5) - case (7) !order 6 - monom = q_ua_monom(u, a, 3)*w_ee_monom(x1, y1, x2, y2, 1) - case (8) - monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 2) - case (9) - monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 3) - case (10) - monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 4) - case (11) - monom = q_ua_monom(u, a, 2)*w_ee_monom(x1, y1, x2, y2, 5) - case (12) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 6) - case (13) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 7) - case (14) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 8) - case (15) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 9) - case (16) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 10) - case (17) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 11) - case (18) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 12) - case (19) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 13) - case (20) - monom = q_ua_monom(u, a, 1)*w_ee_monom(x1, y1, x2, y2, 14) - case default - error stop 'called case of qw_uaee_monom not implemented' - monom = 0.0_dp - end select - end function qw_uaee_monom - - function qz_uaee_monom(u, a, x1, y1, x2, y2, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x1, y1, x2, y2 - integer(idp), intent(in) :: o - select case (o) - case (1) !order 4 - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 1) - case (2) !order 5 - monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 1) - case (3) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 2) - case (4) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 3) - case (5) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 4) - case (6) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 5) - case (7) !order 6 - monom = q_ua_monom(u, a, 3)*z_ee_monom(x1, y1, x2, y2, 1) - case (8) - monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 2) - case (9) - monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 3) - case (10) - monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 4) - case (11) - monom = q_ua_monom(u, a, 2)*z_ee_monom(x1, y1, x2, y2, 5) - case (12) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 6) - case (13) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 7) - case (14) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 8) - case (15) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 9) - case (16) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 10) - case (17) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 11) - case (18) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 12) - case (19) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 13) - case (20) - monom = q_ua_monom(u, a, 1)*z_ee_monom(x1, y1, x2, y2, 14) - case default - error stop 'called case of qz_uaee_monom not implemented' - monom = 0.0_dp - end select - end function qz_uaee_monom - - function qw_uae_monom(u, a, x, y, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x, y - integer(idp), intent(in) :: o - select case (o) - case (1) ! order 3 - monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 1) - case (2) ! order 4 - monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 1) - case (3) - monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 2) - case (4) !order 5 - monom = q_ua_monom(u, a, 3)*w_e_monom(x, y, 1) - case (5) - monom = q_ua_monom(u, a, 4)*w_e_monom(x, y, 1) - case (6) - monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 2) - case (7) - monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 3) - case (8) !order 6 - monom = q_ua_monom(u, a, 5)*w_e_monom(x, y, 1) - case (9) - monom = q_ua_monom(u, a, 6)*w_e_monom(x, y, 1) - case (10) - monom = q_ua_monom(u, a, 3)*w_e_monom(x, y, 2) - case (11) - monom = q_ua_monom(u, a, 4)*w_e_monom(x, y, 2) - case (12) - monom = q_ua_monom(u, a, 2)*w_e_monom(x, y, 3) - case (13) - monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 4) - case (14) - monom = q_ua_monom(u, a, 1)*w_e_monom(x, y, 5) - case default - error stop 'called case of qw_uae not implemented' - monom = 0.0_dp - end select - end function qw_uae_monom - - function qz_uae_monom(u, a, x, y, o) result(monom) - real(dp) :: monom - real(dp), intent(in) :: u, a, x, y - integer(idp), intent(in) :: o - select case (o) - case (1) ! order 3 - monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 1) - case (2) ! order 4 - monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 1) - case (3) - monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 2) - case (4) !order 5 - monom = q_ua_monom(u, a, 3)*z_e_monom(x, y, 1) - case (5) - monom = q_ua_monom(u, a, 4)*z_e_monom(x, y, 1) - case (6) - monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 2) - case (7) - monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 3) - case (8) !order 6 - monom = q_ua_monom(u, a, 5)*z_e_monom(x, y, 1) - case (9) - monom = q_ua_monom(u, a, 6)*z_e_monom(x, y, 1) - case (10) - monom = q_ua_monom(u, a, 3)*z_e_monom(x, y, 2) - case (11) - monom = q_ua_monom(u, a, 4)*z_e_monom(x, y, 2) - case (12) - monom = q_ua_monom(u, a, 2)*z_e_monom(x, y, 3) - case (13) - monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 4) - case (14) - monom = q_ua_monom(u, a, 1)*z_e_monom(x, y, 5) - case default - error stop 'called case of qz_uae not implemented' - monom = 0.0_dp - end select - end function qz_uae_monom - -end module select_monom_mod diff --git a/src/model/Maik_chngs/sphericalharmonics_mod.f90 b/src/model/Maik_chngs/sphericalharmonics_mod.f90 deleted file mode 100644 index 47af619..0000000 --- a/src/model/Maik_chngs/sphericalharmonics_mod.f90 +++ /dev/null @@ -1,575 +0,0 @@ -! Module contains the spherical harmonics up to l=5 m=-l,..,0,..,l listed on https://en.wikipedia.org/wiki/Table_of_spherical_harmonics from 19.07.2022 -! the functions are implementde by calling switch case function for given m or l value and return the corresdpondig value for given theta and phi -! the functions are split for diffrent l values and are named by P_lm. -! example for l=1 and m=-1 the realpart of the spherical harmonic for given theta and phi -! is returned by calling Re_Y_lm(1,-1,theta,phi) which itself calls the corresponding function P_1m(m,theta) and multilpies it by cos(m*phi) to account for the real part of exp(m*phi*i) -! Attention the legendre polynoms are shifted to account for the missing zero order term in spherical harmonic expansions -module sphericalharmonics_mod - use accuracy_constants, only: dp, idp - implicit none - real(kind=dp), parameter :: PI = 4.0_dp * atan( 1.0_dp ) -contains - !---------------------------------------------------------------------------------------------------- - function Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = Y_1m( m , theta , phi ) - case (2) - y = Y_2m( m , theta , phi ) - case (3) - y = Y_3m( m , theta , phi ) - case (4) - y = Y_4m( m , theta , phi ) - case (5) - y = Y_5m( m , theta , phi ) - case default - write(errmesg,'(A,i0)')& - &'order of spherical harmonics not implemented', l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Y_lm - !---------------------------------------------------------------------------------------------------- - function Re_Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = P_1m( m , theta ) * cos(m*phi) - case (2) - y = P_2m( m , theta ) * cos(m*phi) - case (3) - y = P_3m( m , theta ) * cos(m*phi) - case (4) - y = P_4m( m , theta ) * cos(m*phi) - case (5) - y = P_5m( m , theta ) * cos(m*phi) - case (6) - y = P_6m( m , theta ) * cos(m*phi) - case default - write(errmesg,'(A,i0)')& - &'order of spherical harmonics not implemented', l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Re_Y_lm - !---------------------------------------------------------------------------------------------------- - function Im_Y_lm( l , m , theta , phi ) result( y ) - integer(kind=idp), intent( in ) :: l , m - real(kind=dp), intent( in ) :: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( l ) - case (1) - y = P_1m( m , theta ) * sin(m*phi) - case (2) - y = P_2m( m , theta ) * sin(m*phi) - case (3) - y = P_3m( m , theta ) * sin(m*phi) - case (4) - y = P_4m( m , theta ) * sin(m*phi) - case (5) - y = P_5m( m , theta ) * sin(m*phi) - case (6) - y = P_6m( m , theta ) * sin(m*phi) - case default - write(errmesg,'(a,i0)')& - &'order of spherical harmonics not implemented',l - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Im_Y_lm - - !---------------------------------------------------------------------------------------------------- - !---------------------------------------------------------------------------------------------------- - function Y_1m( m , theta , phi ) result( y ) - integer(kind=idp),intent( in ):: m - real(kind=dp),intent( in ):: theta , phi - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-1) - y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi) - - case (0) - y = 0.5_dp*sqrt(3.0_dp/PI)*cos(theta) - - case (1) - y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta)*cos(phi) - - case default - write(errmesg,'(a,i0)') 'in y_1m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - end select - end function Y_1m - !---------------------------------------------------------------------------------------------------- - function Y_2m(m,theta,phi) result(y) - integer(kind=idp),intent(in):: m - real(kind=dp),intent(in):: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(2.0_dp*phi) - - case (-1) - y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta)*cos(phi) - - case (0) - y=0.25_dp*sqrt(5.0_dp/PI)& - &*(3.0_dp*cos(theta)**2-1.0_dp) - - case (1) - y=-0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta)*cos(phi) - - case (2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(2.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_2m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_2m - !---------------------------------------------------------------------------------------------------- - function Y_3m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-3) - y=0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(3.0_dp*phi) - - case (-2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta)*cos(2.0_dp*phi) - - case (-1) - y=0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi) - - case (0) - y=0.25_dp*sqrt(7.0_dp/PI)& - &*(5.0_dp*cos(theta)**3-3.0_dp*cos(theta)) - - case (1) - y=-0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp)*cos(phi) - - case (2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta)*cos(2.0_dp*phi) - - case (3) - y=-0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(3.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_3m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_3m - !---------------------------------------------------------------------------------------------------- - function Y_4m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(4.0_dp*phi) - - case (-3) - y=0.375_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta)*cos(3.0_dp*phi) - - case (-2) - y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2& - &*(7.0_dp*cos(theta)**2-1)*cos(2.0_dp*phi) - - case (-1) - y=0.375_dp*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7.0_dp*cos(theta)**3& - &-3.0_dp*cos(theta))*cos(phi) - - case (0) - y=(3.0_dp/16.0_dp)/sqrt(PI)& - &*(35.0_dp*cos(theta)**4& - &-30.0_dp*cos(theta)**2+3.0_dp) - - case (1) - y=-0.375_dp*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7.0_dp*cos(theta)**3& - &-3.0_dp*cos(theta))*cos(phi) - - case (2) - y=0.375_dp*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7.0_dp*cos(theta)**2-1.0_dp)& - &*cos(2*phi) - - case (3) - y=-0.375_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta)*cos(3.0_dp*phi) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(4.0_dp*phi) - - case default - write(errmesg,'(a,i0)')'in y_4m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_4m - !---------------------------------------------------------------------------------------------------- - function Y_5m(m,theta,phi) result(y) - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta,phi - real(kind=dp) y - character(len=70) :: errmesg - select case (m) - case (-5) - y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5*cos(5*phi) - - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta)*cos(4*phi) - - case (-3) - y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp)*cos(3*phi) - - case (-2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi) - - case (-1) - y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21.0_dp*cos(theta)**4& - &-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi) - - case (0) - y=(1.0_dp/16.0_dp)/sqrt(11.0_dp/PI)& - &*(63.0_dp*cos(theta)**5-70.0_dp*cos(theta)**3& - &+15.0_dp*cos(theta)) - - case (1) - y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21.0_dp*cos(theta)**4& - &-14.0_dp*cos(theta)**2+1.0_dp)*cos(phi) - - case (2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta))*cos(2*phi) - - case (3) - y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9.0_dp*cos(theta)**2-1.0_dp)& - &*cos(3.0_dp*phi) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta)*cos(4.0_dp*phi) - - case (5) - y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5*cos(5.0_dp*phi) - - case default - write(errmesg,'(A,i0)')'in y_5m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function Y_5m - !---------------------------------------------------------------------------------------------------- - function P_1m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=1 and given m and theta - integer(kind=idp),intent( in ):: m - real(kind=dp),intent( in ):: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-1) - y = 0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta) - - case (0) - y = 0.5_dp*sqrt(3.0_dp/PI)*(cos(theta)-1.0_dp) ! -1 is subtracted to shift so that for theta=0 y=0 - - case (1) - y = -0.5_dp*sqrt(3.0_dp/(PI*2.0_dp))*sin(theta) - - case default - write(errmesg,'(A,i0)')'in p_1m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_1m - !---------------------------------------------------------------------------------------------------- - function P_2m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=2 and given m and theta - integer(kind=idp),intent(in):: m - real(kind=dp),intent(in):: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-2) - y=0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2 - - case (-1) - y=0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta) - - case (0) - y = (3.0_dp*cos(theta)**2-1.0_dp) - y = y - 2.0_dp !2.0 is subtracted to shift so that for theta=0 y=0 - y = y * 0.25_dp*sqrt(5.0_dp/PI) ! normalize - - case (1) - y = -0.5_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)*cos(theta) - - case (2) - y = 0.25_dp*sqrt(15.0_dp/(PI*2.0_dp))& - &*sin(theta)**2 - - case default - write(errmesg,'(A,i0)')'in p_2m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_2m - !---------------------------------------------------------------------------------------------------- - function P_3m( m , theta ) result( y ) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=3 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-3) - y=0.125_dp*sqrt(35.0_dp/PI)& - &*sin(theta)**3 - - case (-2) - y=0.25_dp*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta) - - case (-1) - y=0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5*cos(theta)**2-1.0_dp) - - case (0) - y=(5.0_dp*cos(theta)**3-3*cos(theta)) - y=y-2.0_dp ! 2.0 is subtracted to shift so that for theta=0 y=0 - y=y*0.25_dp*sqrt(7.0_dp/PI) ! normalize - - case (1) - y=-0.125_dp*sqrt(21.0_dp/(PI))& - &*sin(theta)*(5.0_dp*cos(theta)**2-1.0_dp) - - case (2) - y=0.25*sqrt(105.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*cos(theta) - - case (3) - y=-0.125*sqrt(35.0_dp/PI)& - &*sin(theta)**3 - - case default - write(errmesg,'(A,i0)')'in p_3m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_3m - !---------------------------------------------------------------------------------------------------- - function P_4m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=4 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4 - - case (-3) - y=0.375*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta) - - case (-2) - y=0.375*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7*cos(theta)**2-1) - - case (-1) - y=0.375*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7*cos(theta)**3-3*cos(theta)) - - - case (0) - y=(35*cos(theta)**4-30*cos(theta)**2+3) - y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0 - y = y * (3.0_dp/16.0_dp)/sqrt(PI) - - case (1) - y=-0.375*sqrt(5.0_dp/(PI))& - &*sin(theta)*(7*cos(theta)**3-3*cos(theta)) - - case (2) - y=0.375*sqrt(5.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(7*cos(theta)**2-1) - - case (3) - y=-0.375*sqrt(35.0_dp/PI)& - &*sin(theta)**3*cos(theta) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(35.0_dp/2.0_dp*PI)& - &*sin(theta)**4 - - case default - write(errmesg,'(A,i0)')'in p_4m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_4m - !---------------------------------------------------------------------------------------------------- - function P_5m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=5 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp) y - character(len=70) :: errmesg - select case ( m ) - case (-5) - y=(3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5 - - case (-4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta) - - case (-3) - y=(1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp) - - case (-2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta)) - - case (-1) - y=(1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1) - - - case (0) - y = (63*cos(theta)**5-70*cos(theta)**3+15*cos(theta)) - y = y - 8.0_dp !8.0 is subtracted to shift so that for theta=0 y=0 - y = y * (1.0_dp/16.0_dp)/sqrt(11.0_dp/PI) - - case (1) - y=(-1.0_dp/16.0_dp)*sqrt(165.0_dp/(2.0_dp*PI))& - &*sin(theta)*(21*cos(theta)**4-14*cos(theta)**2+1) - - case (2) - y=0.125*sqrt(1155.0_dp/(PI*2.0_dp))& - &*sin(theta)**2*(3*cos(theta)**3-cos(theta)) - - case (3) - y=(-1.0_dp/32.0_dp)*sqrt(385.0_dp/PI)& - &*sin(theta)**3*(9*cos(theta)**2-1.0_dp) - - case (4) - y=(3.0_dp/16.0_dp)*sqrt(385.0_dp/2.0_dp*PI)& - &*sin(theta)**4*cos(theta) - - case (5) - y=(-3.0_dp/32.0_dp)*sqrt(77.0_dp/PI)& - &*sin(theta)**5 - - case default - write(errmesg,'(A,i0)')'in p_5m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_5m - !---------------------------------------------------------------------------------------------------- - function P_6m(m,theta) result(y) - ! >Function returns the value of the corresponding normalized Associated legendre polynom for l=6 and given m and theta - integer(kind=idp), intent(in) :: m - real(kind=dp), intent(in) :: theta - real(kind=dp):: y - character(len=70) :: errmesg - select case ( m ) - case (-6) - y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)& - &* sin(theta)**6 - case (-5) - y = (3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)& - &* sin(theta)**5& - &* cos(theta) - case (-4) - y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))& - &* sin(theta)**4& - &* (11*cos(theta)**2 - 1 ) - case (-3) - y= (1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**3& - &* (11*cos(theta)**3 - 3*cos(theta) ) - case (-2) - y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**2& - &* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 ) - case (-1) - y= (1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))& - &* sin(theta)& - &* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) ) - case (0) - y = 231*cos(theta)**6 - 315*cos(theta)**4 + 105*cos(theta)**2-5 - y = y - 16.0_dp !16.0 is subtracted to shift so that for theta=0 y=0 - y = y * (1.0_dp/32.0_dp)*sqrt(13.0_dp/PI) - case (1) - y= -(1.0_dp/16.0_dp)*sqrt(273.0_dp/(2.0_dp*PI))& - &* sin(theta)& - &* (33*cos(theta)**5 - 30*cos(theta)**3 + 5*cos(theta) ) - case (2) - y= (1.0_dp/64.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**2& - &* (33*cos(theta)**4 - 18*cos(theta)**2 + 1 ) - case (3) - y= -(1.0_dp/32.0_dp)*sqrt(1365.0_dp/PI)& - &* sin(theta)**3& - &* (11*cos(theta)**3 - 3*cos(theta) ) - case (4) - y= (3.0_dp/32.0_dp)*sqrt(91.0_dp/(2.0_dp*PI))& - &* sin(theta)**4& - &* (11*cos(theta)**2 - 1 ) - case (5) - y= -(3.0_dp/32.0_dp)*sqrt(1001.0_dp/PI)& - &* sin(theta)**5 * cos(theta) - case (6) - y = (1.0_dp/64.0_dp)*sqrt(3003.0_dp/PI)& - &* sin(theta)**6 - case default - write(errmesg,'(A,i0)')'in p_6m given m not logic, ', m - error stop 'error in spherical harmonics' !error stop errmesg - - end select - end function P_6m - !---------------------------------------------------------------------------------------------------- - -end module diff --git a/src/model/adia.f90 b/src/model/adia.f90 index e81f01b..78a0c91 100644 --- a/src/model/adia.f90 +++ b/src/model/adia.f90 @@ -66,9 +66,9 @@ ! init eigenvector matrix - TYPES = int(p(pst(1,32))) + TYPES = int(p(pst(1,33))) - BLK = int(p(pst(1,32)+1)) ! BLOCK IF TYPE IS 3 + BLK = int(p(pst(1,33)+1)) ! BLOCK IF TYPE IS 3 u = 0.d0 vx=0.0d0 skip=.false. diff --git a/src/model/data_transform.f90 b/src/model/data_transform.f90 index 67bf585..411a09a 100644 --- a/src/model/data_transform.f90 +++ b/src/model/data_transform.f90 @@ -32,14 +32,14 @@ - if (pst(2,32) .ne. 2) then + if (pst(2,33) .ne. 2) then - write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,32) + write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,33) stop end if - TYPES = int(p(pst(1,32)))! TYPE OF THE CALCULATION - BLK= int(p(pst(1,32)+1))! BLOCK IF TYPE IS 3 + TYPES = int(p(pst(1,33)))! TYPE OF THE CALCULATION + BLK= int(p(pst(1,33)+1))! BLOCK IF TYPE IS 3 write(*,*) "TYPE of calculation:",TYPES pt=1 @@ -109,7 +109,7 @@ call one_dia_upper(mat_z,y(1:ntot,pt)) else write(*,*) "Error in TYPE of calculationss",TYPES - write(*,*) "the value:,", p(pst(1,32)) + write(*,*) "the value:,", p(pst(1,33)) stop end if pt=pt+1 diff --git a/src/model/keys.f90 b/src/model/keys.f90 index ea9d2d1..4cffdd9 100644 --- a/src/model/keys.f90 +++ b/src/model/keys.f90 @@ -10,7 +10,7 @@ contains character(len=1) prefix(4) parameter (prefix=['N','P','A','S']) !character (len=20) key(4,25) - integer,parameter:: np=32 + integer,parameter:: np=33 character(len=16) parname(np) integer i,j ! Defining keys for potential @@ -63,11 +63,12 @@ contains parname(27)='LZPE1E2O1' parname(28)='LZPE1E2O2' parname(29)='LZPA2E1O1' - parname(30)='LZPA2E2O2' + parname(30)='LZPA2E1O2' parname(31)='LZPA2E2O1' + parname(32)='LZPA2E2O2' - parname(32)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE + parname(33)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE do i=1,np do j=1,4 diff --git a/src/model/matrix_form.f90 b/src/model/matrix_form.f90 index 71d335c..18fddf1 100644 --- a/src/model/matrix_form.f90 +++ b/src/model/matrix_form.f90 @@ -12,7 +12,7 @@ y=0.0d0 !y(1)=mx(4,4)+mx(5,5) - do i=2,3 + do i=1,ndiab y(1)=y(1)+mx(i,i) y(2)=y(2)+my(i,i) enddo @@ -298,8 +298,8 @@ double precision, intent(in) :: p(:) integer, intent(in) :: id_write integer :: type_calc, blk - type_calc = int(p(pst(1,32))) - blk = int(p(pst(1,32)+1)) + type_calc = int(p(pst(1,33))) + blk = int(p(pst(1,33)+1)) if (type_calc ==1) then write(id_write,*) "Type of calculation: TRACE" diff --git a/src/model/model.f90 b/src/model/model.f90 index 93e7d5e..c970e80 100644 --- a/src/model/model.f90 +++ b/src/model/model.f90 @@ -64,7 +64,8 @@ module diab_mod id=id+1 ! 4 ! order 2 - e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 5 p + + e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 3 p +p(pst(1,id)+2)*(xs*xb-ys*yb) id =id+1 ! 5 e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & @@ -388,13 +389,13 @@ module diab_mod ! w and z of E'' ! order 1 - id = id + id = id ! 22 e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb e(3,3) = e(3,3) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb ! order 2 - id = id +1 + id = id +1 ! 23 do i =1,3 e(2,2) = e(2,2) + p(pst(1,id)+(i-1))*v2(i+3) e(3,3) = e(3,3) - p(pst(1,id)+(i-1))*v2(i+3) @@ -404,13 +405,13 @@ module diab_mod ! W and Z of E' ! order 1 - id = id +1 + id = id +1 ! 24 e(4,4) = e(4,4) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb e(5,5) = e(5,5) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb e(4,5) = e(4,5) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb ! order 2 - id = id +1 + id = id +1 ! 25 do i =1,3 e(4,4) = e(4,4) + p(pst(1,id)+(i-1))*v2(i+3) e(5,5) = e(5,5) - p(pst(1,id)+(i-1))*v2(i+3) @@ -420,7 +421,7 @@ module diab_mod ! Pseudo of E' and E'' ! it must have odd power of b - id = id +1 + id = id +1 !26 ! order 0 e(2,4) = e(2,4) e(3,5) = e(3,5) @@ -428,13 +429,13 @@ module diab_mod e(3,4) = e(3,4) - b*(p(pst(1,id))) ! order 1 - id = id +1 + id = id +1 !27 e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) ! order 2 - id = id +1 + id = id +1 !28 do i=1,3 e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3) e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3) @@ -445,11 +446,12 @@ module diab_mod ! the coupling between A2' and E'' ! order 1 - id = id +1 + id = id +1 !29 e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb)) e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb)) - - id = id +1 + + + id = id +1 !30 ! order 2 do i=1,3 e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i) @@ -459,11 +461,11 @@ module diab_mod ! the coupling of A2' and E' ! order 1 - id = id +1 + id = id +1 !31 e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb)) e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb)) - id = id +1 + id = id +1 ! 32 ! order 2 do i=1,3 e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i) diff --git a/src/model/trans/JTmod.incl b/src/model/trans/JTmod.incl deleted file mode 100644 index 77e5f69..0000000 --- a/src/model/trans/JTmod.incl +++ /dev/null @@ -1,38 +0,0 @@ -!*** Relevant parameters for the analytic model -!*** offsets: -!*** offsets(1): morse equilibrium (N-H, Angström) -!*** offsets(2): reference angle (H-N-H) -!*** offsets(3): -- -!*** pat_index: vector giving the position of the -!*** various coordinates (see below) -!*** ppars: polynomial parameters for tmcs -!*** vcfs: coefficients for V expressions. -!*** wzcfs: coefficients for W & Z expressions. -!*** ifc: inverse factorials. - - integer matdim - parameter (matdim=5) ! matrix is (matdim)x(matdim) - - real*8 offsets(2) - integer pat_index(maxnin) - -! NH3 params - parameter (offsets=[1.0228710942d0,120.d0]) - -!########################################################################## -! coordinate order; the first #I number of coords are given to the -! ANN, where #I is the number of input neurons. The position i in -! pat_index corresponds to a coordinate, the value of pat_index(i) -! signifies its position. -! -! The vector is ordered as follows: -! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2, -! es*eb, es**3, eb**3,es**2*eb, es*eb**2 -! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b -! -! parts not supposed to be read by ANN are marked by ';' for your -! convenience. -!########################################################################## -! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D) - parameter (pat_index=[1,2,3,4,5,6,7,8,9,10,11,12,13,14]) -!************************************************************************** diff --git a/src/model/trans/cart2int.f b/src/model/trans/cart2int.f deleted file mode 100644 index f5d9e65..0000000 --- a/src/model/trans/cart2int.f +++ /dev/null @@ -1,260 +0,0 @@ - subroutine cart2int(cart,qint) - implicit none -! This version merges both coordinate transformation routines into -! one. JTmod's sscales(2:3) are ignored. -! This is the first version to be compatible with one of my proper 6D fits -! Time-stamp: <2024-10-22 13:52:59 dwilliams> - -! Input (cartesian, in Angström) -! cart(:,1): N -! cart(:,1+i): Hi -! Output -! qint(i): order defined in JTmod. -! Internal Variables -! no(1:3): NO distances 1-3 -! pat_in: temporary coordinates -! axis: main axis of NO3 - include 'nnparams.incl' - include 'JTmod.incl' - - - real*8 cart(3,4),qint(maxnin) - - real*8 no(3), r1, r2, r3 - real*8 v1(3), v2(3), v3(3) - real*8 n1(3), n2(3), n3(3), tr(3) - real*8 ortho(3) - real*8 pat_in(maxnin) - logical ignore_umbrella,dbg_umbrella - logical dbg_distances - -!.. Debugging parameters -!.. set umbrella to 0 - parameter (ignore_umbrella=.false.) -! parameter (ignore_umbrella=.true.) - -!.. break if umbrella is not 0 - parameter (dbg_umbrella=.false.) -! parameter (dbg_umbrella=.true.) - -!.. break for tiny distances - parameter (dbg_distances=.false.) -! parameter (dbg_distances=.true.) - - integer k - -!.. get N-O vectors and distances: - do k=1,3 - v1(k)=cart(k,2)-cart(k,1) - v2(k)=cart(k,3)-cart(k,1) - v3(k)=cart(k,4)-cart(k,1) - enddo - no(1)=norm(v1,3) - no(2)=norm(v2,3) - no(3)=norm(v3,3) - -!.. temporarily store displacements - do k=1,3 - pat_in(k)=no(k)-offsets(1) - enddo - - do k=1,3 - v1(k)=v1(k)/no(1) - v2(k)=v2(k)/no(2) - v3(k)=v3(k)/no(3) - enddo - -!.. compute three normal vectors for the ONO planes: - call xprod(n1,v1,v2) - call xprod(n2,v2,v3) - call xprod(n3,v3,v1) - - do k=1,3 - tr(k)=(n1(k)+n2(k)+n3(k))/3.d0 - enddo - r1=norm(tr,3) - do k=1,3 - tr(k)=tr(k)/r1 - enddo - -! rotate trisector - call rot_trisec(tr,v1,v2,v3) - -!.. determine trisector angle: - if (ignore_umbrella) then - pat_in(7)=0.0d0 - else - pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3)) - pat_in(7)=sign(pat_in(7),cart(1,2)) - endif - -!.. molecule now lies in yz plane, compute projected ONO angles: - v1(1)=0.d0 - v2(1)=0.d0 - v3(1)=0.d0 - r1=norm(v1,3) - r2=norm(v2,3) - r3=norm(v3,3) - do k=2,3 - v1(k)=v1(k)/r1 - v2(k)=v2(k)/r2 - v3(k)=v3(k)/r3 - enddo - -! make orthogonal vector to v3 - ortho(1)=0.0d0 - ortho(2)=v3(3) - ortho(3)=-v3(2) - -!.. projected ONO angles in radians - pat_in(4)=get_ang(v2,v3,ortho) - pat_in(5)=get_ang(v1,v3,ortho) - - pat_in(6)=dabs(pat_in(5)-pat_in(4)) - -!.. account for rotational order of atoms - if (pat_in(4).le.pat_in(5)) then - pat_in(5)=2*pi-pat_in(4)-pat_in(6) - else - pat_in(4)=2*pi-pat_in(5)-pat_in(6) - endif - - pat_in(4)=rad2deg*pat_in(4)-offsets(2) - pat_in(5)=rad2deg*pat_in(5)-offsets(2) - pat_in(6)=rad2deg*pat_in(6)-offsets(2) - pat_in(7)=rad2deg*pat_in(7) - - call genANN_ctrans(pat_in) - - qint(:)=pat_in(:) - - contains - -!------------------------------------------------------------------- -! compute vector product n1 of vectors v1 x v2 - subroutine xprod(n1,v1,v2) - implicit none - - real*8 n1(3), v1(3), v2(3) - - n1(1) = v1(2)*v2(3) - v1(3)*v2(2) - n1(2) = v1(3)*v2(1) - v1(1)*v2(3) - n1(3) = v1(1)*v2(2) - v1(2)*v2(1) - - end subroutine - - -!------------------------------------------------------------------- -! compute scalar product of vectors v1 and v2: - real*8 function scalar(v1,v2,n) - implicit none - integer i, n - real*8 v1(*), v2(*) - - scalar=0.d0 - do i=1,n - scalar=scalar+v1(i)*v2(i) - enddo - - end function - -!------------------------------------------------------------------- -! compute norm of vector: - real*8 function norm(x,n) - implicit none - integer i, n - real*8 x(*) - - norm=0.d0 - do i=1,n - norm=norm+x(i)**2 - enddo - norm=sqrt(norm) - - end function - -!------------------------------------------------------------------- - - subroutine rot_trisec(tr,v1,v2,v3) - implicit none - - real*8 tr(3),v1(3),v2(3),v3(3) - - real*8 vrot(3) - real*8 rot_ax(3) - real*8 cos_phi,sin_phi - -! evaluate cos(-phi) and sin(-phi), where phi is the angle between -! tr and (1,0,0) - cos_phi=tr(1) - sin_phi=dsqrt(tr(2)**2+tr(3)**2) - - if (sin_phi.lt.1.0d-12) then - return - endif - -! determine rotational axis - rot_ax(1) = 0.0d0 - rot_ax(2) = tr(3) - rot_ax(3) = -tr(2) -! normalize - rot_ax=rot_ax/sin_phi - -! now the rotation can be done using Rodrigues' rotation formula -! v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p)) -! for v=tr k*v vanishes by construction: - -! check that the rotation does what it should - call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi) - if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then - write(6,*) "ERROR: BROKEN TRISECTOR" - stop - endif - tr=vrot - - call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi) - v1=vrot - call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi) - v2=vrot - call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi) - v3=vrot - - - end subroutine - -!------------------------------------------------------------------- - - subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi) - implicit none - real*8 vrot(3),v(3),axis(3) - real*8 cos_phi,sin_phi - - real*8 ortho(3) - - call xprod(ortho,axis,v) - vrot = v*cos_phi + ortho*sin_phi - > + axis*scalar(axis,v,3)*(1-cos_phi) - - end subroutine - -!------------------------------------------------------------------- - - real*8 function get_ang(v,xaxis,yaxis) - implicit none -! get normalized [0:2pi) angle from vectors in the yz plane - real*8 v(3),xaxis(3),yaxis(3) - - real*8 phi - - real*8 pi - parameter (pi=3.141592653589793d0) - - phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3)) - if (phi.lt.0.0d0) then - phi=2*pi+phi - endif - get_ang=phi - - end function - - end subroutine cart2int diff --git a/src/model/trans/ctrans.f b/src/model/trans/ctrans.f deleted file mode 100644 index ca619dd..0000000 --- a/src/model/trans/ctrans.f +++ /dev/null @@ -1,88 +0,0 @@ -!------------------------------------------------------------------- -! Time-stamp: "2024-10-09 13:33:50 dwilliams" - - subroutine genANN_ctrans(pat_in) - implicit none - - include 'nnparams.incl' - include 'JTmod.incl' - - real*8 pat_in(maxnin) - - real*8 raw_in(maxnin),off_in(maxnin),ptrans_in(7) - real*8 r0 - real*8 a,b,xs,ys,xb,yb - - integer k - - off_in(1:7)=pat_in(1:7) - r0=offsets(1) - -! transform primitives -! recover raw distances from offset coords - do k=1,3 - raw_in(k)=off_in(k)+offsets(1) - enddo - - do k=1,3 - ptrans_in(k)=off_in(k) - enddo - -! rescale ONO angles - ptrans_in(4)=deg2rad*off_in(4) - ptrans_in(5)=deg2rad*off_in(5) - ptrans_in(6)=deg2rad*off_in(6) -! rescale umbrella - ptrans_in(7)=off_in(7)*deg2rad - -! compute symmetry coordinates - -! A (breathing) - a=(ptrans_in(1)+ptrans_in(2)+ptrans_in(3))/dsqrt(3.0d0) -! ES - call prim2emode(ptrans_in(1:3),xs,ys) -! EB - call prim2emode(ptrans_in(4:6),xb,yb) -! B (umbrella) - b=ptrans_in(7) - -! overwrite input with output - - pat_in(pat_index(1))=a ! 1 - pat_in(pat_index(2))=xs - pat_in(pat_index(3))=ys - pat_in(pat_index(4))=xb - pat_in(pat_index(5))=yb - pat_in(pat_index(6))=b -! totally symmetric monomials - pat_in(pat_index(7))=xs**2 + ys**2 ! 2 - pat_in(pat_index(8))=xb**2 + yb**2 ! 3 - pat_in(pat_index(9))=b**2 ! 9 - pat_in(pat_index(10))=xs*xb+ys*yb ! 4 -! S^3, B^3 - pat_in(pat_index(11))=xs*(xs**2-3*ys**2) ! 5 - pat_in(pat_index(12))=xb*(xb**2-3*yb**2) ! 6 -! S^2 B, S B^2 - pat_in(pat_index(13))=xb*(xs**2-ys**2) - 2*yb*xs*ys ! 7 - pat_in(pat_index(14))=xs*(xb**2-yb**2) - 2*ys*xb*yb ! 8 - - do k=11,14 - pat_in(pat_index(k))=tanh(0.1d0*pat_in(pat_index(k)))*10.0d0 - enddo - - contains - - subroutine prim2emode(prim,ex,ey) - implicit none -! Takes a 2D-vector prim and returns the degenerate modes x and y -! following our standard conventions. - - real*8 prim(3),ex,ey - - ex=(2.0d0*prim(1)-prim(2)-prim(3))/dsqrt(6.0d0) - ey=(prim(2)-prim(3))/dsqrt(2.0d0) - - end - - - end subroutine diff --git a/src/model/trans/nnparams.incl b/src/model/trans/nnparams.incl deleted file mode 100644 index 635c33c..0000000 --- a/src/model/trans/nnparams.incl +++ /dev/null @@ -1,43 +0,0 @@ -!**** Declarations - - real*8 pi - real*8 hart2eV, eV2hart - real*8 hart2icm, icm2hart - real*8 eV2icm, icm2eV - real*8 deg2rad, rad2deg - integer maxnin,maxnout - -!********************************************************** -!**** Parameters -!*** maxnin: max. number of neurons in input layer -!*** maxnout: max. number of neurons in output layer - - parameter (maxnin=14,maxnout=15) - -!********************************************************** -!**** Numerical Parameters -!*** infty: largest possible double precision real value. -!*** iinfty: largest possible integer value. - - ! 3.14159265358979323846264338327950... - parameter (pi=3.1415926536D0) - -!********************************************************** -!**** Unit Conversion Parameters -!*** X2Y: convert from X to Y. -!*** -!*** hart: hartree -!*** eV: electron volt -!*** icm: inverse centimeters (h*c/cm) -!**** -!*** deg: degree -!*** rad: radians - - parameter (hart2icm=219474.69d0) - parameter (hart2eV=27.211385d0) - parameter (eV2icm=hart2icm/hart2eV) - parameter (icm2hart=1.0d0/hart2icm) - parameter (eV2hart=1.0d0/hart2eV) - parameter (icm2eV=1.0d0/eV2icm) - parameter (deg2rad=pi/180.0d0) - parameter (rad2deg=1.0d0/deg2rad) diff --git a/src/model/write.f b/src/model/write.f index ec3b3f4..700fd30 100644 --- a/src/model/write.f +++ b/src/model/write.f @@ -400,10 +400,10 @@ call print_plotheader(id_plot,coord,npt,set) do i=1,npt - write(id_plot,"(ES16.8,*(3(ES16.8),:))") - > x(coord,i), ymod(:,i),y(:,i),(wt(:,i)) -! write(id_plot,"(2ES16.8,*(3(ES16.8),:))") -! > x(coord,i), x(coord+1,i),y(:,i) +! write(id_plot,"(ES16.8,*(3(ES16.8),:))") +! > x(coord,i), ymod(:,i),y(:,i),(wt(:,i)) + write(id_plot,"(2ES16.8,*(3(ES16.8),:))") + > x(coord,i), x(coord+1,i),y(:,i) enddo end subroutine