diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 7103f14..aff88c9 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -955,11 +955,9 @@ contains
     real(wp_) :: gr2
     real(wp_), dimension(3) :: dgr2
     
-!    if(igrad.eq.1) then
-      gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
-      dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
-!    end if
-    fk1 = yp
+    gr2  = dot_product(dgr, dgr)
+    dgr2 = 2 * matmul(ddgr, dgr)
+    fk1  = yp
 
     yy = y + fk1*hh
     call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2,igrad)
@@ -1025,8 +1023,8 @@ contains
     real(wp_), dimension(3,3) :: derbv
     real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_
 
-    gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
-    dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
+    gr2  = dot_product(dgr, dgr)
+    dgr2 = 2 * matmul(ddgr, dgr)
     call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi)
     call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
        dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
@@ -1259,7 +1257,7 @@ contains
 ! local variables
     integer :: jv
     real(wp_) ::  xx,yy,zz
-    real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm
+    real(wp_) :: csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm
     real(wp_), dimension(3) :: dbtot,bvc
     real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
     real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi
@@ -1367,10 +1365,9 @@ contains
     dbv(:,3) = dbvdc(:,3)
 
 !     B magnitude and derivatives 
-    b2tot = bv(1)**2 + bv(2)**2 + bv(3)**2
-    btot  = sqrt(b2tot)
+    btot  = norm2(bv)
 
-    dbtot = (bv(1)*dbv(1,:) + bv(2)*dbv(2,:) + bv(3)*dbv(3,:))/btot
+    dbtot = matmul(bv, dbv)/btot
 
     yg = btot/Bres
 
@@ -1417,8 +1414,8 @@ contains
     real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm
     real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv
 
-    an2  = anv(1)*anv(1) + anv(2)*anv(2) + anv(3)*anv(3)
-    anpl = anv(1)*bv(1)  + anv(2)*bv(2)  + anv(3)*bv(3)
+    an2  = dot_product(anv, anv)
+    anpl = dot_product(anv, bv)
 
     anpl2 = anpl**2
     dnl   = one - anpl2
@@ -1467,7 +1464,7 @@ contains
       end if
     end if
 
-    bdotgr = bv(1)*dgr(1) + bv(2)*dgr(2) + bv(3)*dgr(3)
+    bdotgr = dot_product(bv, dgr)
     do iv=1,3
       dbgr(iv)     = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv)                 &
                    + dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv)                 &
@@ -1481,22 +1478,23 @@ contains
                *(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl)
     derdnv   = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv
 
-    derdnm   = sqrt(derdnv(1)**2 + derdnv(2)**2 + derdnv(3)**2)
+    derdnm   = norm2(derdnv)
 
     derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl     &
              + two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg               &
                                           + half*yg*dfdiadyg               &
                                           + half*anpl*dfdiadnpl)
 
+    ! integration variable
     if (idst == 0) then
-!   integration variable: s
+      ! optical path: s
       denom = derdnm
     else if (idst == 1) then
-!   integration variable: c*t
+      ! "time": c*t
       denom = -derdom
     else
-!   integration variable: Sr
-      denom = anv(1)*derdnv(1) + anv(2)*derdnv(2) + anv(3)*derdnv(3)
+      ! real eikonal: S_R
+      denom = dot_product(anv, derdnv)
     end if
 
 !   coefficient for integration in s