Update fun.f90

suggestion for optimization

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