 #endif
         // inversely rotate momentum vector.
 #if NDIM == 3
-        pjsoln[1] = matinv[0][0]*mom[0] + matinv[0][1]*mom[1]
-                  + matinv[0][2]*mom[2];
-        pjsoln[2] = matinv[1][0]*mom[0] + matinv[1][1]*mom[1]
-                  + matinv[1][2]*mom[2];
-        pjsoln[3] = matinv[2][0]*mom[0] + matinv[2][1]*mom[1]
-                  + matinv[2][2]*mom[2];
-#else
-        pjsoln[1] = matinv[0][0]*mom[0] + matinv[0][1]*mom[1];
-        pjsoln[2] = matinv[1][0]*mom[0] + matinv[1][1]*mom[1];
+        pjsoln[1] = mvt[0][0]*mom[0] + mvt[0][1]*mom[1]
+                  + mvt[0][2]*mom[2];
+        pjsoln[2] = mvt[1][0]*mom[0] + mvt[1][1]*mom[1]
+                  + mvt[1][2]*mom[2];
+        pjsoln[3] = mvt[2][0]*mom[0] + mvt[2][1]*mom[1]
+                  + mvt[2][2]*mom[2];
+#else
+        pjsoln[1] = mvt[0][0]*mom[0] + mvt[0][1]*mom[1];
+        pjsoln[2] = mvt[1][0]*mom[0] + mvt[1][1]*mom[1];
 #endif
         // set solutions.
         pjsoln[0] = pisoln[0];

        mat[1][2] = pndcrd[2] - pfccnd[2];
        len = sqrt(mat[1][0]*mat[1][0] + mat[1][1]*mat[1][1]
                 + mat[1][2]*mat[1][2]);
        mat[1][0] = matinv[0][1] = mat[1][0]/len;
        mat[1][1] = matinv[1][1] = mat[1][1]/len;
        mat[1][2] = matinv[2][1] = mat[1][2]/len;
        mat[2][0] = matinv[0][2] = mat[0][1]*mat[1][2] - mat[0][2]*mat[1][1];
        mat[2][1] = matinv[1][2] = mat[0][2]*mat[1][0] - mat[0][0]*mat[1][2];
        mat[2][2] = matinv[2][2] = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
#else
        mat[1][0] = matinv[0][1] =  pfcnml[1];
        mat[1][1] = matinv[1][1] = -pfcnml[0];
#endif
        // rotate momentum vector.
#if NDIM == 3

        pjdsoln = exd->dsoln + jcl*NEQ*NDIM;
        // coordinate transformation and set transformed vectors.
        pfcnml = exd->fcnml + ifc*NDIM;
        mat[0][0] = mvt[0][0] = pfcnml[0];
        mat[0][1] = mvt[1][0] = pfcnml[1];
#if NDIM == 3
        mat[0][2] = mvt[2][0] = pfcnml[2];
        pfcnds = exd->fcnds + ifc*(FCMND+1);
        pndcrd = exd->ndcrd + pfcnds[1]*NDIM;
        pfccnd = exd->fccnd + ifc*NDIM;
