From 562f67dabddabd60baaa0787702ae88548acfd49 Mon Sep 17 00:00:00 2001
From: "Martin D. Weinberg" <weinberg@astro.umass.edu>
Date: Wed, 3 Jan 2024 14:21:20 -0500
Subject: [PATCH] Compute velocity second moments in VelocityBasis by default
 [no ci]

---
 coefs/FieldBasis.cc | 19 ++++++++++++++++---
 1 file changed, 16 insertions(+), 3 deletions(-)

diff --git a/coefs/FieldBasis.cc b/coefs/FieldBasis.cc
index 55173ec60..5b4ee1436 100644
--- a/coefs/FieldBasis.cc
+++ b/coefs/FieldBasis.cc
@@ -640,7 +640,7 @@ namespace BasisClasses
     double vr  = (u*x + v*y)/R;
     double vp  = (u*y - v*x)/R;
     
-    return {vr, w, vp};
+    return {vr, w, vp, vr*vr, w*w, vp*vp};
   }
 
   std::vector<double> sphVel(double mass,
@@ -656,13 +656,14 @@ namespace BasisClasses
     double vt  = (u*z*x + v*z*y - w*R)/R/r;
     double vp  = (u*y - v*x)/R;
 
-    return {vr, vt, vp};
+    return {vr, vt, vp, vr*vr, vt*vt, vp*vp};
   }
 
   std::vector<double> crtVel(double mass,
 			     FieldBasis::PS3& pos, FieldBasis::PS3& vel)
   {
-    return {vel[0], vel[1], vel[2]};
+    auto [u, v, w] = vel;
+    return {u, v, w, u*u, v*v, w*w};
   }
 
 
@@ -679,21 +680,33 @@ namespace BasisClasses
       fieldLabels.push_back("v_R");
       fieldLabels.push_back("v_z");
       fieldLabels.push_back("v_p");
+      fieldLabels.push_back("v_R^2");
+      fieldLabels.push_back("v_z^2");
+      fieldLabels.push_back("v_p^2");
       fieldFunc = cylVel;
     } else if (coordinates == Coord::Cartesian) {
       fieldLabels.push_back("v_x");
       fieldLabels.push_back("v_y");
       fieldLabels.push_back("v_z");
+      fieldLabels.push_back("v_x^2");
+      fieldLabels.push_back("v_y^2");
+      fieldLabels.push_back("v_z^2");
       fieldFunc = crtVel;
     } else if (coordinates == Coord::None) {
       fieldLabels.push_back("v_x");
       fieldLabels.push_back("v_y");
       fieldLabels.push_back("v_z");
+      fieldLabels.push_back("v_x^2");
+      fieldLabels.push_back("v_y^2");
+      fieldLabels.push_back("v_z^2");
       fieldFunc = crtVel;
     } else {
       fieldLabels.push_back("v_r");
       fieldLabels.push_back("v_t");
       fieldLabels.push_back("v_p");
+      fieldLabels.push_back("v_r^2");
+      fieldLabels.push_back("v_t^2");
+      fieldLabels.push_back("v_p^2");
       fieldFunc = sphVel;
     }