diff --git a/source/f90/fun.f90 b/source/f90/fun.f90
index d01e398da5d731b7f22e6997860d71a47f4c104e..73df08d284045a0f43005a0b328567993735ff24 100644
--- a/source/f90/fun.f90
+++ b/source/f90/fun.f90
@@ -22,5 +22,8 @@ double precision function fun(phi)
   fun = dsqrt((1.0d0 - elleps + elleps * sinphi**2) *   &
               (1.0d0 - sinpsi * sinphi) *               &
               (1.0d0 + sinpsi * sinphi))
+!  optimized but untested version
+!  sinphi2 = dsin(phi)**2
+!  fun = dsqrt((1.0d0 + elleps*(sinphi2 - 1.0d0)) * (1.0d0 - sinpsi**2 * sinphi2))
   return
 end function fun