From d3cbc671099cb2da802ee4ad6097cd1bbf7064fe Mon Sep 17 00:00:00 2001
From: kamischi <karl-michael.schindler@web.de>
Date: Wed, 6 Mar 2024 13:55:20 +0100
Subject: [PATCH] New tests and summary

---
 .gitignore                                 |   1 +
 tests/phint optimization/README & RESULTS  |  53 ++++++++++++++++++
 tests/phint optimization/dooptimization.sh |  10 ----
 tests/phint optimization/phint10.f90       |  49 +++++++++++++++++
 tests/phint optimization/phint8.f90        |   8 +--
 tests/phint optimization/phint9.f90        |  61 +++++++++++++++++++++
 tests/phint optimization/phintdriver       | Bin 34616 -> 51288 bytes
 tests/phint optimization/phintdriver.f90   |  25 ++++++++-
 8 files changed, 192 insertions(+), 15 deletions(-)
 create mode 100644 tests/phint optimization/README & RESULTS
 delete mode 100755 tests/phint optimization/dooptimization.sh
 create mode 100644 tests/phint optimization/phint10.f90
 create mode 100644 tests/phint optimization/phint9.f90

diff --git a/.gitignore b/.gitignore
index ee9ddc1..7319c4b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,4 @@ tests/oxlog/bosou
 tests/oxlog/bosonf90
 tests/oxlog/eelsou
 tests/phintlog/eelsf90
+tests/phint optimization/phintdriver
diff --git a/tests/phint optimization/README & RESULTS b/tests/phint optimization/README & RESULTS
new file mode 100644
index 0000000..21ec55a
--- /dev/null
+++ b/tests/phint optimization/README & RESULTS	
@@ -0,0 +1,53 @@
+These programs compare different versions of phint.
+
+Phint is calculated for 
+  0 <= phi <= π
+  0.2 <= u <= 2.4
+  a = 1.73205
+
+Phint: original from Lambin
+Phint1: new formula 
+Phint2: new formula, use conjugate complex
+Phint3: new formula, partial formula ?
+Phint4: new formula, partial formula ?
+Phint5: new formula, 
+Phint6: new formula, consider conjugate complex property
+Phint7: new formula, consider conjugate complex property, extract common sub-expressions 
+Phint8: new formula, consider conjugate complex property, extract common sub-expressions, treat phi = π
+Phint9: new formula, consider conjugate complex property, extract common sub-expressions, treat phi = π, use rational approximation of original Lambin code
+Phint10: new formula, consider conjugate complex property, no common sub-expression, treat phi = π, leave everything to the compiler.
+
+Log:
+
+>% gfortran -O0 phintdriver.f90 phint.f90 phint1.f90 phint2.f90 phint3.f90 phint4.f90 phint5.f90 phint6.f90 phint7.f90 phint8.f90 phint9.f90 phint10.f90 -o phintdriver
+>% phintdriver  
+Time for phint:  9.277 s
+Time for phint1: 17.785 s
+Time for phint2: 14.707 s
+Time for phint3:  8.899 s
+Time for phint4: 17.914 s
+Time for phint5: 25.697 s
+Time for phint6: 12.598 s
+Time for phint7:  6.052 s
+Time for phint8:  6.455 s
+Time for phint9:  6.594 s
+Time for phint10:  7.036 s
+
+>% gfortran -O7 phintdriver.f90 phint.f90 phint1.f90 phint2.f90 phint3.f90 phint4.f90 phint5.f90 phint6.f90 phint7.f90 phint8.f90 phint9.f90 phint10.f90 -o phintdriver
+>% phintdriver
+Time for phint:  8.648 s
+Time for phint1: 13.632 s
+Time for phint2: 13.302 s
+Time for phint3:  8.099 s
+Time for phint4: 11.682 s
+Time for phint5: 14.071 s
+Time for phint6:  8.823 s
+Time for phint7:  5.859 s
+Time for phint8:  5.008 s
+Time for phint9:  5.023 s
+Time for phint10:  4.963 s
+
+Summary:
+Considering the conjugate complex property is a major plus. As expected.
+Manually extracting common sub-expression is a plus with -O0, but no advantage, maybe even a small disadvantage with -O7. Not really expected.
+The rational approximation of the original Lambin code does not give any advantage.
diff --git a/tests/phint optimization/dooptimization.sh b/tests/phint optimization/dooptimization.sh
deleted file mode 100755
index 65b4fb9..0000000
--- a/tests/phint optimization/dooptimization.sh	
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/sh -v
-for number in 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 ; do
-    echo '*** File number' $number '***'
-
-    cp ../inputFiles/'eelsin'$number eelsin
-    ./eelsf90 >/dev/null
-    cp phint.log logFiles/'phint'$number'.log'
-
-    rm -f eelsin eelsou
-done
diff --git a/tests/phint optimization/phint10.f90 b/tests/phint optimization/phint10.f90
new file mode 100644
index 0000000..93965b9
--- /dev/null
+++ b/tests/phint optimization/phint10.f90	
@@ -0,0 +1,49 @@
+double precision function phint10(phi, a, u)
+
+! ******************************************************************
+! *                                                                *
+! * evaluate the integral from zero to phi of                      *
+! *                                                                *
+! *                 u                 2                            *
+! *  ( ----------------------------- )  dphi                       *
+! *                          2     2                               *
+! *    (1 - a * u * cos(phi))  +  u                                *
+! *                                                                *
+! * for 0 <= phi <= pi , u >= 0 and a >= 0                         *
+! *                                                                *
+! ******************************************************************
+
+  implicit none
+
+  double precision, intent(in) :: phi
+  double precision, intent(in) :: a
+  double precision, intent(in) :: u
+
+  double precision :: e, pi, tanPhi2, sinPhi, cosPhi, Term1
+
+  double complex :: preFactor
+
+  if (phi == 0) then
+  	phint10 = 0
+  else
+    e = a*u
+    pi = 4 * atan (1.0d0)
+	preFactor = cmplx((e**2 - 1)/u + 2*u, 3) / sqrt(cmplx(e**2 - 1 + u**2, 2*u))**3
+	if (pi - phi < 1e-6) then
+	  phint10 = pi / 2 * real(preFactor)
+	else
+	  tanPhi2 = tan(phi/2)
+	  sinPhi = sin(phi)
+	  cosPhi = cos(phi)
+      Term1 = real(preFactor * atan(cmplx(u, 1 + e) / sqrt(cmplx(e**2 - 1 + u**2, 2*u)) * tanPhi2))
+      
+      phint10 = 																					&
+        Term1 +																					&
+          0.5 * e * (e**2 - 1 + 3*u**2 - e*(e**2 - 1 + u**2)*cosPhi) * sinPhi /					&
+            (((e**2 - 1)**2 + 2*(1 + e**2)*u**2 + u**4) * (1 + u**2 + e*cosPhi*(-2 + e*cosPhi)))
+     endif
+
+  endif
+
+  return
+end function phint10
diff --git a/tests/phint optimization/phint8.f90 b/tests/phint optimization/phint8.f90
index 7b2d140..647090c 100644
--- a/tests/phint optimization/phint8.f90	
+++ b/tests/phint optimization/phint8.f90	
@@ -40,11 +40,11 @@ double precision function phint8(phi, a, u)
 	  tanPhi2 = tan(phi/2)
 	  sinPhi = sin(phi)
 	  cosPhi = cos(phi)
-	  Term1 = real(preFactor	* atan(cmplx(u,1+e) / denominator * tanPhi2))
+	  Term1 = real(preFactor * atan(cmplx(u,1+e) / denominator * tanPhi2))
 	
-	  phint8 = 																	&
-	    Term1 +																	&
-		  0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /				&
+	  phint8 = 																		&
+	    Term1 +																		&
+		  0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /					&
 		    ((e21**2 + 2*(1 + e2)*u2 + u2**2) * (1 + u2 + e*cosPhi*(-2 + e*cosPhi)))
     endif
 
diff --git a/tests/phint optimization/phint9.f90 b/tests/phint optimization/phint9.f90
new file mode 100644
index 0000000..8e12f40
--- /dev/null
+++ b/tests/phint optimization/phint9.f90	
@@ -0,0 +1,61 @@
+double precision function phint9(phi, a, u)
+
+! ******************************************************************
+! *                                                                *
+! * evaluate the integral from zero to phi of                      *
+! *                                                                *
+! *                 u                 2                            *
+! *  ( ----------------------------- )  dphi                       *
+! *                          2     2                               *
+! *    (1 - a * u * cos(phi))  +  u                                *
+! *                                                                *
+! * for 0 <= phi <= pi , u >= 0 and a >= 0                         *
+! *                                                                *
+! ******************************************************************
+
+  implicit none
+
+  double precision, intent(in) :: phi
+  double precision, intent(in) :: a
+  double precision, intent(in) :: u
+
+  double precision :: e, e2, e21, pi, u2, zr, esr, tanPhi2, sinPhi, cosPhi, Term1
+
+  double complex :: nominator, denominator, preFactor
+
+  if (phi == 0) then
+  	phint9 = 0
+  else
+    e = a*u
+	e2 = e*e
+    u2 = u*u
+	e21 = e2 - 1
+    pi = 4 * atan (1.0d0)
+	nominator   =      cmplx(e21/u + 2*u ,3)
+	denominator = sqrt(cmplx(e21   +   u2,2*u))
+	preFactor   = nominator / denominator**3
+	if (pi - phi < 1e-6) then
+	  phint9 = pi / 2 * real(preFactor)
+	else
+	  tanPhi2 = tan(phi/2)
+	  sinPhi = sin(phi)
+	  cosPhi = cos(phi)
+      if (u < 1.0d0 .and. e < 1.0d-02 * (1.0d0 + u2)) then
+        zr = 1.0d0 + u2
+        esr = e / zr
+        phint9 = u2 / zr**2 * ((( (4.0d0 / 3.0d0) * (2.0d0 + cosPhi**2) * sinPhi * (5.0d0 - 3 * u2) *  &
+                esr + (phi + cosPhi * sinPhi) * (5.0d0 - u2)) * esr + 4 * sinPhi) * esr + phi)
+      else
+        Term1 = real(preFactor * atan(cmplx(u,1+e) / denominator * tanPhi2))
+        
+        phint9 = 																		&
+          Term1 +																		&
+            0.5 * e * (e21 + 3*u2 - e*(e21 + u2)*cosPhi) * sinPhi /						&
+              ((e21**2 + 2*(1 + e2)*u2 + u2**2) * (1 + u2 + e*cosPhi*(-2 + e*cosPhi)))
+       endif
+     endif
+
+  endif
+
+  return
+end function phint9
diff --git a/tests/phint optimization/phintdriver b/tests/phint optimization/phintdriver
index 5d6deeb066e287ef7622eaa0cfe8fb5f4694a901..bc3454f117149dc1a0762e451cb747f18e07d273 100755
GIT binary patch
delta 5883
zcmdnd$8=)?^Mp>p1_o5nFmaZd5Q7^710w?igNF`O(rw}+JE4Fgh$u`jplGrYqn%I<
zRQLc?ZB5bSCPq6Uhhm6o0jRJ;@#IU4cA^1L;R#F-aS)$*vLKuEWCtb%p$DZ9O$uBP
z*$1VQi<rzOJ2Hz;KEPzD(g2Z#ySagZff=F@OtLXBFbGVRWmb_q08z-`02XIpfVwS$
zV{#O;yw3p!h6Y9k24e<>29WC-4l*!oU|?VXxk!S6fg#`!1A{yR1H*wM5Ov%jF{nWx
zH6YB(z`(!{woZ^un$_V51B1fEK;g+hm^B0pepytk2{1|Q?QK(ftN!rEWCIoprhuZ&
zSuACYY6gcG7!(*77&Jk8LFOKU7%U6rI~-$RkYL)}=+JCGS%Iy5ast1iT#wUV1_p)?
zFVxv`IVV8apMx1!eqPKpkz<O}<oW!%Of1}!&++RqnKDiO$1hv|;Q_PEs`p$BQ<|9@
zf*2VXJ}5fyPGn+fxCAoh@nwdc%nS@49&!t=1j#YXaQf@Q%sWwYg45rh%*+$BdYt|~
zV3u6<AerIs@5xLPpQJPVlwfEGdd0;s#r#oA(1Z41U645snWa}XF*F2mOo6CpuuznH
zu$*D<GcEy!DIfzLaxqMK;O_AAEZ5{!0!EB$C*Kv2uHVVUFvW?PcOuAU#RI$(m7%VB
zl+N&z3FIn%gsYygGz57;T?29*$Q2-0fn4$AGQ-ZTAPWe%W)nyav*fC^Allu5`R7cq
zYnT$5C*Kp0s&7ShP68iZ=eV%pcTOY73Ifil1F2z_TvZ98LC#6#nj9r&%m{H!AhK&1
z`0=`C2|Iq*_=0R8;2JNG8fM8=t{@uZ8be6n#Ik?`M;6&R3-CLKg%iJXBtcdXaE=&A
z4YTAbK@jck!1VL)<dp)(OlKjkdC7^CavTKk22KhWe%CzZWSBxI9X|r8VU}EV7es?x
za}w;DNH%cTY(;j>1N^Qz!j0cG8$mV@2%B{vHO!K$R)T1dYo>x-!^sYIO(U{v8U*o%
zjS4S**VKY+AmExRkQ!#mRiz*rlyDL`C%+OlW`wxL7uhumLU>)%!iV2Ao*)|txW)~n
zhFNlzBZvltjV?HBjzU}`iR_vK_+9gaAHQovK{gO@jSxr;v*apX5Df~OpB$j_kI9!~
z@;w15v!@&kQ$WQhxByHLp6D#t5CkfMomd(!flAFNX5NXv(iwh&3OZ0x4JxLeb4<P`
zCdKuP16r_hOlA^~V!FmLIZoV-X&=YrW#V#D+c+4eFgqm3F+5?=Pv*#zTd0{V_t>1_
z?>3If_r$fu7C8MiP)LyT;bREt;bRD4U|_hwvBYVzyo4pwbdJeU5^8!gIT)t!IwZ(3
z3*^W(I3~zFG<W#BmNQwd5o$1_V**$_!y!SghGX(NiN#EFoF)fLK4enhm~1Ac#H7YC
zIZDcz%Sm+NCr^L*MXQ`9FPDmCde1(YO?sC6A$Eo-5)KJ+UoU&-nmHuMEpQT@sNk3&
z=jfOq2U5CUT9j!w`{ZlVbxc#)H;2iTF^XBSGfd$F`A8r~j)5~l?p-^3t_%C*N3zGc
znnWjla_7uV+T=8Ovz#i^FSg0I<dV#8voTB&c1Vyjbx4r=7Vq$v!7)KDp-FUN(k7?B
z8C+R%S{zw&E0n_IvsRVKGqX7SedzA+_X^wOB6%;S-E5op$@?*iWw9|#;fC7z(3;_|
z0%w9;3EN~P#g)nz7!Qa}d^(#wmuZL7UlyrMxtE#+a>^_Ye?Pi2{C#51@Yj@W@;gO8
zCTX_Kc1n7TOb=Nn=P9dconvK~!sL)3$0?8}2MTq4hXlD!%}hCl15SSrXeP+{IV8ww
zIwr}TWu3fFS&eBI>*Qz3X-v~uC;O?$G4-=f&QmFqD`aJuqU?|$_e9&_ucM)-Jf}l~
zoWsOi`7DP7xh&SnOsb!>4OkhbbTT>oWf91c%Vu%_H5e}RaHhyTXy?k+a7>cZ=FFB;
z)y$M*KH@Z4UQLxrn{~3AS``xy>*PafeN6XRCI_i|F&$=^JWsueX*J7aC5<?yE|$$r
z8hwmRYAln5w7Qs>SvJqp@?&BuVxIg?M~kV9d9sx5CME~w$>(&vm{gc2OX>A7aWYR{
zrKiXAg=zCGy+w>neN2<<47|h|nHZ+<zyj5t;qL*?1i4nG$@dJLn0%QgD;c&iIWTQr
zW$47n^qX<=J0mHk%Z!t`j7@nqFq{Clb~Z4cahmLHEGu=H5fYXR0y%P84heD%oC$J1
z4heF17$-LwCoyegoczpqDaVEv>W;Y#7n~+fHnFeYu@aOa{xUO6e5KEn%gx~MH$gE$
zPJy2xM415;T6uCW)ERR<7##k-F=x)@;!Ku%XU~$$q?sYd&zT|Tz?miI=#V7GYmp-N
zz?k9hI>`*V2b#HZ5A+%SGB_m3ozqN`JLi}rckZB)TskAe6wU)mata5O<PsPff*2hV
z<P_>b8A+SruRSBf6t)8ia!Lmh<Sr;A$+d7Zgs4d*$ngjy$|*3o$se-tkw3H1P2Nc&
zL9RzK0hCfUFkFH9{)W?E28(Pt0fj`lBU}t2AP*%vCdfT9XZY*Lh{z+L^qRn#Am`1<
zFhxKhPmWC>QLg@xxx-%p#R9o!oN00&Gz;XOa598E)OPsGDv&2<crZbZ!!b$DRWnsi
z@j!wc$3Z2zr4C7Q21<!?OAjc?E#YAZd93a5H&G%}E>p!(exrh;JnK|Pc|)m8xup`B
za*R(H^f@?E<sNA>{QY3g@b{TL!{4PGiE@fcj`9ps3gyKa8S1BaaOTN9(uY|1%--QI
zr$CO}2XlwNiX1s|&Ky~CkBlAuF4Iht^VUd}V--k{OPk>*@64Gg$7%^O*x_%sX0{w-
zOOpHpYlpuAj;`{c$j#==ms_IiD}Rh5S?-OwbMA9NfB9n^33Bi3U2+Au(&dgx6v!RW
zNR)e}U+<LrxBA421I+Q)*d!9c!Oh6xCV#-tP5#`D9C;>*M7cAPiE>N-{}*RqV1T5Y
zJ5GOJ={x0u40eDT3Nr(whVg;ZUq-G(IRk|(IUfdwkhTB+i#tme$T12e$Q@+$lYj7;
zVbg+t&p&Z6q|bX=Z|(4xX+n}516P3@i$Ic`n?;)3IvqcGCyRW!B^)Vo4<r-h6ge{F
zK)wQ5sK}8j=fsf*_K6}#nw%nsv%I2!v%IQ~v-|>%1i3T+|BHj7L{TzXPEgcQUYXZX
ze#L@Bc|)mWIabMJxjFy;i~p@gB>H;fRL{W25TXK2>AnmOe?cjI9cKzS2XF}F$t~c_
zm0PKqD|d=B51iT+9h2ldEHmU7G!x_=m_t(hdIp9mIt~hQ{F1rgRQ@2EVefm*Ou6@t
zNpdX=phT`9r+Pp^?t+4<+!HQ_`Vch<1vyRuML7iqHu<AAoboT1vdJ%%P>|D<Q~;&&
z37~|=@C2OFCosN1q;wW8h7eFH2c`C>_6&bJ7#ODT!wL~=hQFYcK8b;0O4~sNxrfmX
zdzl3C<eu6){5_<ZB&U8rL2jXAlAM!6eUjV<1vR-N91J0kqaF6PODfBC^Rmmkv9QZ4
zb+F4HlhBrP0^tKn;M9EVfRfx(eTKj890_tx9EoxVIuhg$$~elWF)&Oy&6zL9CXgq`
zB#<cg*xcbSyFh}Rmu8yWBF-ear-xm0L6vU0Wx5=XK#rW3W|rI{MJIW-daX3MBb+&M
zng^BS(j?vGy*P8_k~;k5K`B*}Gf{5kk_2$d28nZMddNF*q{#6q1;~RE-Z73`aO!4|
zNRwmG$ONYmP{M>IMm7mWId%a>IYt&X`KPjM^5^%c%P*BslzS<uD7OldpeLk(5*EW7
zM2ZF_B9P@ULqMUDE}<i*#=$P1W56!&*rG0fOhQLaNkT_%H^dNd2la3D{*|CoREU9L
z3L^u<1jP@+6CW@%1PL-QOnLfc@;j450s;&S8wBheeljr0Ujymcyumbwak7q?2&3iX
zkEZHEh6Yv&N;(Q@W_rdt3Q7vaN}7|U%-2ocU@pO4Ewax2VO7iI_2!mb3>W@`K@R8S
zPv((KES!^lEtHu~{GME7;l{M!#^hZV@0bd1O<rYrjmhHf<R+_QOdIY`_Od?5WbkmZ
zl+89Kj>nTP*_>mtcsY5I?J*{XH<PpMjxkldoor-(Y;uB?gy0(n1_l`h28IXkz}Ec#
zFj>(?G?0gZfq~VFht->#1vCV}>dDLM!N<Y@67}U|_2Ggr{n=UlI9TQ|GBB`tvU#w8
z1`pVL*?d?)y4n2M{8&U77#Qx|W%FY5p4@Ju7RaK(z`(!-Rtr`FQVTK@td<3?mXCo2
z#NlR%fhz~`QH=&Ui+eJky%Oh%4-l7opt~GoJ1d5BK^8*&Qvh-~$ekc}A&dw41*8h$
za(xnAKKYM*6|+DC|K>u6ZhKzP07n1=1H%N+pvGju1_M9%hz(?1hLMSZfl-fvfiZ%C
zfw75!fpHZB1LF;l93umx5hDY00V4x*9U}vC4<iHf3`Pd#C5#Nrn;03G4>2+@pJANb
z-=M*45}!HwK+C~;P*P)HV2o#I1C3aL#%5R;E<9mmU}0!@03tUmVPs%nGKgmoj%Sc$
zVYsjwBsgIch<vaOM1I%}CQpLMhKnHbz*P`A;TDMeFprUefy<CV-YA|y!#JKn&m^9~
z#5A74$}FD2!91RU$IT+1!N)S5At;_9V)B_rrTRFCix?Vu7(WOdU|qm|fb9ay1l9(=
z3G4@0F0g!Hy#R__28K*V1_nNmo1uekg&;mF1H+#ZP(VP2qxcy?X$-;_0}VMbF)+Y}
z;XsL)AEdtkBCiILXJn{nFbIMO=tBe;E`&n(7Er!GB!uq><tIc#_+C){hgb+d7{*V5
z@Ik2)WMM-Fgr5qPKad6C=Yje43=9If5P@>A00TpTKZFm8D@F(t97a<az#Nbt!2DHc
z{QYSB%V_*pX#9U5em&Ssuz3vP5ZmAum_hjPF+*E4z6%=P3ymLu#t%p1$1^f8NHQV~
zN~Rzyn0&U`oD&plpus*lG@fmX8k0f%<meWEF3Wg^_=1egyps6I6B^W^3H_3jp`o4u
z1H-R-JO>yUm_Vfg%v&PA?(r}&Ffg%z5Cao<95aQ1fgu7C2n@Uo496T4m>75_3$m%1
z>3Bee8Nqg>6lIpB7BR$DyFQ+v>K*s^NB_<4q+Q>$R^Pgls<nxEIo}`YtL@8{PTtX$
zP@ge#k4No-kgl_(0W4>D?0%*Q_WO!QrN-Ub=Uox>%C<ZFYC~6E1B2y0{b^_4KVXcu
zEXg|HJ)0@unpfDZ7fO6vR;td8;mI;HSu$t-+I88VHdMdjyY+o*-Z9R=%wz?7mQ^Bl
zM-%x3|7%Q);=XJl{M)OnChT$D!s4f#ucPx_T`CnV*La*M7w<l#TytmIbZy`4-^KFV
zCx)Bp2>AW|Q@tU{#s=g{g}*mdJ8w3Bzid7C*vVZ}?6Xo*Vv1Go<$fyOB~l(KGN~N^
Dw6~tz

delta 3684
zcmcaHfq6$C(}Yez2L@E&FmaZdPy#;#10w?iLxw6uGJ${MBRioBCJ<2_h#<oSlgUPm
zc0zxk!dsxie@rGfG1>`jFomcVfC_IgoqUPWPV@j&n1cx-4&pOU7G#s2?7*ZTbiy2>
z$$$eQd%}Ej5tBK{ks#+Vz#Qu^c^8wtiUUM#1w<(WNR}DKW?*1oV_;z5pDfF)BG~{@
z2vPu+Zvc}F7VMLwnB{#M7#K1b85oQi7#cus&1hs`=wM)A0J%?sfq{XeiGe|efq}uH
zg@GZ0k%57mfq@|aVgOhTNF7Lw7i^s%n>6c*W(EeiiGjkCe=utZEVYn#sMuz7ZI}C;
zdNuobe<vHTSTJ2M*__2v#;Eq8k%0kZhb9991ITWWT_B@nq5K7{3=A_EHw&^waZYkm
zl$+!Hmw|!d!wYq`T+Rg$_UB;6m7f<gP2^bOJlWn&m&u27a-N$W(^<yJ^W0?XFFatD
zS!K@2FeRC}A&8NI;esNA_(Ud_hD!_#3>O|>X4uKhz;NLqx8O>U9K#Cdzb=QxCTcEl
z{`-@ed1BTa=f4k_C09L2X88MiGSkE-=?p(57#f1~IT@yyKWYhj&>pM{GUp+)^r|L?
zh9Hh55cLceigFK@Gwjvo6kwPFGC+-!Vafw{ho7RHlfB%H7+ojVx=YtHaxzSDIxIF3
zWV2#}_(Wx>s~)8@{A2>TiXY*sBP<O;UQpM7TnBOm$W<U$Jh{xU^A`uh6hf}~4pPG`
zx#|;$c6VU@`H%zbnxA0ToaI1rje-PT*Ql`Jcg-=74FtmGFh~uv<f?rj8swU_VAn8B
zWu9E;E>%Ah**ORBJEw&mzjG#mtRUc=evlew$yHq-8swZxh;zO|oRf;|oCHa{q4R_T
zzjNY1RuFJbG)N7z<f<?b4RVev$K+o=#*8fBz%fL24TBV3*SK)uca0Xv1_G{82dQC}
zT%`n}L9P*m1kOiRuyg*hBPE>$_?@$a8^3eDu@g?npFwJvC0D%z(IDsCWd|i>riE-^
z=bS`#j)OGb&|%@l@0>#*D+q+nevlew$yK{RG{`wC*(bm9FlGe1#%wA(!xT^%0nX(Q
zq$eixH3Wgm86}p6OQ515>af_vU+D}#LAe=J9)ZfDDeRNw{H3@igA8PrT-C=u*~veO
zshNH9JbyQ)Z1%~|{N<#Q*%_uVJ0!?4JYmpJ=E##<sF^JH*qq^SGW%q?0Bx}i&VLOQ
z66AdN7(#mZ7(zfU;Mn3kIX=LW$%TFLtN=AVH+F_8ybcL+%mO)b4UP$N56vC^uH{UY
z(}Nn!=$HT&&u~bP(`29gE?_a!8t2KA10OP-Wt&_Uq{MWYZSt%jXRbrC6F+(S%P-pH
zJo$N0EYlpe$!@{3<n!4Wrbsv>$bG%+oonWhAh+O<>_i2}1UW~?1UZn>|G}b68Eliq
zLh6{D*fviKDPt6S&dM-_59A|(961Kg1i5$Z?71IVCo6>==VFqZ_{p6!H|db`<j-NM
zOv_m(ONA$yb+9r_5q3zBGj&Li`xfu;m%%YXE`dpIV$vb!zZqOva#|c&ax0X=<g-?l
z$uqM!{C()|@VAk5@}_VvrVQ52|HAzk#l%?|rf@^;d}z(^SAjD@PLXwTQsheI4U7(Q
z6Q9mz&t*E{{Fg;4Q|_f^ft)gn!{3kY41b@PGyHwTGFdOmkLd)<=DH|7My5WN$?Kxk
zw5nMcrZ71q$Z-nf$uTf6Y|wW|kn7aUlw-Ky{P%!nf}9^HI31JZs#qrfi&kSwXPK-P
zlg8x2GPy5Cj>(2)^17HZIXM=FDasBBa!<4!{yG|Z%5yp-$T>{RmCtfWkP~N_>=gS+
z`z|xXlujmxzbpbda@kA{pjvT54`+(pgLbZ54aX!oZO&{tRn1H}<}1#V<KtABZZJ=7
zi>qSV$vl}UzK^MgdGe%qFQx+K$@k)$m?D@bCndx&nKN%bl+ee>beU;#P*NAu2ByvT
zlKhyM<e4VxrD`!LGfj?4-Nf{kaWY?;7t=+?$x&&2OxqYIze>|%TFAIrDt!?nlQrYy
zeHmV2dW;NHcwm8Q&+zvEXM&s&<7Bx^C#GKvlan&rnBFpMewFFO$h3lCvR<|nQv<_f
zuWVDE2@E&DRr&<RJI<54vt^|k7$9NEAdn-c<&YrPz?mTD0}4-u$%nF&nBo{FtK}@^
znD9c~F_+<i^W@7p_A>id{$*yE$j`tqg^___f?|Qp#0LxwL3|7hQ|5k~oL6*cvqNzX
z<K#ajYbHM^Rp6?){vQm^uutYHi)2b+pPX8z%#`qFa#xufQ^mc>m&)ETc|4eWsQem}
zz~jlYDvmLAJe{0Wd5)9gIauHSp2?2&qLclqW-&c@G5J;1IVOSklXq1gW4iHia#PJQ
zrhw0ry=sq54CkB7)=<F2(6CvMt;e2M0os0=0BJ8BFz|!5>KGu+G)5){21Y#w2F3^m
z21ZcRZ5;yx<242bMh->>W;aF#<`6~(<|IZ2<~&9Q<|;-8<_<;%<{69(%!|Nn#L4{!
zG$x-puql}_o*@O)m;yD>m>Di)GB7iIFlJz8m|)Go%y7Yxftg`}8v`@L25%4%$-vBT
zAc29Ip&=E-)nj1bG>m5uGKyyqH;!kJF^OkTG>vCaGmB@?Hjig8u!v_cjc2f$%yvl0
z#Gc`S+y(vvtP9u=uw7u8z}mn!f&Bo>1r`PdM@9w)K1K!xSaaGF#Ajt-kTIW}bx5ur
z)}98nF8LW57(Pfs0sz$L2Jt7zLil>1kOGyU3J|_Il&_!+;e)~tq<(`cgzpKJ-=G8G
z2SND`1`vKUgkR6lU<45W1rx{s191pH2NDzD-~@$W6B>Uy8h<q!{}3Ag8XEsSm|qX>
zVt~vC2NEL#1HAPv$_RD{oG*jMS3=`!pz-z5_~vMQ8xUUr$wE-$oo}-35pz~JYx033
z%9Gz7ac7L5JmG-)<TD4>Pu4sp=w+#A0CqBj_;rux00RT#46uhFLL$HJ@qkir3j-4i
zs22fBexLx3fJifNGBC6*P+(%<nJmbrR<D`?6=p2R$jmEADatHMEn;xLFPanVuyMPE
zLSxG|zZ<(mwwGn^HvRSK@RsM2i>-n}oI2NQf9E;*DezR!*3<RellI)muy+nA@VFJg
z)|b6|_O4fJ6DKb$PhWX;b-==!6`5-HI-|n(mai&_UL5|q{RO}E#``D2>{u9QhnY*r
z{8zGJ&=lX0^0TzFW4`hpVF%;nnt6r~pJu-@=f1EgXlw1_+4<G?XLw#a@#i#y&b;g1
vChXc;lW!bLP~Z$%rfc;2vPQwgr;-Ns_j5znwyrw(CHCa@=2D)jg-%`o?tCxC

diff --git a/tests/phint optimization/phintdriver.f90 b/tests/phint optimization/phintdriver.f90
index e5caf4d..686307d 100644
--- a/tests/phint optimization/phintdriver.f90	
+++ b/tests/phint optimization/phintdriver.f90	
@@ -3,7 +3,8 @@ program phintdriver
   implicit none
 
   double precision :: phi, a, u, result
-  double precision :: phint, phint1, phint2, phint3, phint4, phint5, phint6, phint7, phint8
+  double precision :: phint, phint1, phint2, phint3, phint4, phint5
+  double precision :: phint6, phint7, phint8, phint9, phint10
   double precision :: step_phi, step_u
   integer :: i, j, number
   
@@ -114,4 +115,26 @@ program phintdriver
   call cpu_time(finish)
   print '("Time for phint8: ", f6.3, " s")',finish-start
   
+  call cpu_time(start)
+  do i = 0, number
+    phi = i * step_phi
+    do j = 0, number
+      u = 0.2 + j * step_u
+      result = phint9(phi, a, u)
+    enddo
+  enddo
+  call cpu_time(finish)
+  print '("Time for phint9: ", f6.3, " s")',finish-start
+  
+  call cpu_time(start)
+  do i = 0, number
+    phi = i * step_phi
+    do j = 0, number
+      u = 0.2 + j * step_u
+      result = phint10(phi, a, u)
+    enddo
+  enddo
+  call cpu_time(finish)
+  print '("Time for phint10: ", f6.3, " s")',finish-start
+  
   end program phintdriver
-- 
GitLab