Add extraQuanities to FieldlinesState from cdf source

modules/fieldlinessequence/util/fieldlinesstate.cpp

@@ -122,6 +122,155 @@ bool FieldlinesState::addLinesFromKameleon(ccmc::Kameleon* kameleon,
 }
 #endif // OPENSPACE_MODULE_KAMELEON_ENABLED
 
+#ifdef OPENSPACE_MODULE_KAMELEON_ENABLED
+void FieldlinesState::loadExtrasIntoKameleon(ccmc::Kameleon* kameleon,
+                                             std::vector<std::string>& xtraScalarVars,
+                                             std::vector<std::string>& xtraMagVars) {
+    // Load the existing SCALAR variables into kameleon.
+    // Remove non-existing variables from vector
+    for (int i = 0; i < xtraScalarVars.size(); i++) {
+        std::string& str = xtraScalarVars[i];
+        bool isSuccesful = kameleon->doesVariableExist(str) && kameleon->loadVariable(str);
+        if (!isSuccesful &&
+            (_model == fls::Model::BATSRUS && (str == T_AS_P_OVER_RHO || str == "T" ))) {
+            LDEBUG("BATSRUS doesn't contain variable T for temperature. Trying to "
+                   << "calculate it using the ideal gas law: T = pressure/density");
+            const std::string P = "p", R = "rho";
+            isSuccesful = kameleon->doesVariableExist(P) && kameleon->loadVariable(P)
+                       && kameleon->doesVariableExist(R) && kameleon->loadVariable(R);
+            str = T_AS_P_OVER_RHO;
+        }
+        if (!isSuccesful) {
+            LWARNING("FAILED TO LOAD EXTRA VARIABLE: '" << str << "'. Ignoring it!");
+            xtraScalarVars.erase(xtraScalarVars.begin() + i);
+            --i;
+        } else {
+            _extraQuantityNames.push_back(str);
+        }
+    }
+
+    // Load the existing magnitude variables (should be provided in multiple of 3)
+    // into kameleon. Remove non-existing variables from vector
+    if (xtraMagVars.size() % 3 == 0) {
+        for (int i = 0; i < static_cast<int>(xtraMagVars.size()); i += 3) {
+            std::string s1 = xtraMagVars[i];
+            std::string s2 = xtraMagVars[i+1];
+            std::string s3 = xtraMagVars[i+2];
+            bool isSuccesful = kameleon->doesVariableExist(s1) &&
+                               kameleon->doesVariableExist(s2) &&
+                               kameleon->doesVariableExist(s3) &&
+                               kameleon->loadVariable(s1) &&
+                               kameleon->loadVariable(s2) &&
+                               kameleon->loadVariable(s3);
+            std::string name = "Magnitude of (" + s1 + ", "+ s2 + ", "+ s3 + ")";
+            if (isSuccesful && _model == fls::Model::BATSRUS && s1 == "jx" && s2 == "jy"
+                    && s3 == "jz") {
+                // CCMC isn't really interested in the magnitude of current, but by the
+                // magnitude of the part of the current's vector that is parallel to the
+                // magnetic field => ensure that the magnetic variables are loaded
+                isSuccesful =  kameleon->doesVariableExist("bx") &&
+                               kameleon->doesVariableExist("by") &&
+                               kameleon->doesVariableExist("bz") &&
+                               kameleon->loadVariable("bx") &&
+                               kameleon->loadVariable("by") &&
+                               kameleon->loadVariable("bz");
+                name = J_PARALLEL_B;
+            }
+            if (!isSuccesful) {
+                LWARNING("FAILED TO LOAD AT LEAST ONE OF THE MAGNITUDE VARIABLES: "
+                        << s1 << ", " << s2 <<  " & " << s3
+                        << ". Removing ability to store corresponding magnitude!");
+                xtraMagVars.erase(xtraMagVars.begin() + i, xtraMagVars.begin() + i + 3);
+                i -= 3;
+            } else {
+                _extraQuantityNames.push_back(name);
+            }
+        }
+    } else {
+        // WRONG NUMBER OF MAGNITUDE VARIABLES.. REMOVE ALL!
+        xtraMagVars.clear();
+        LWARNING("Wrong number of variables provided for storing magnitudes. "
+                << "Expects multiple of 3 but " << xtraMagVars.size()
+                << " are provided");
+    }
+}
+#endif // OPENSPACE_MODULE_KAMELEON_ENABLED
+
+#ifdef OPENSPACE_MODULE_KAMELEON_ENABLED
+/**
+ * Loops through _vertexPositions and extracts corresponding 'extraQuantities' att each
+ * position from the kameleon object using a ccmc::interpolator.
+ * Note that the positions MUST be unaltered (NOT scaled NOR converted to different
+ * coordinate system)!
+ *
+ * @param kameleon raw pointer to an already opened Kameleon object
+ * @param xtraScalarVars vector of strings. Strings should be names of a scalar quantities
+ * to load into _extraQuantites; such as: "T" for temperature or "rho" for density.
+ * @param xtraMagVars vector of strings. Size must be multiple of 3. Strings should be
+ * names of the components needed to calculate magnitude. E.g. {"ux", "uy", "uz"} will
+ * calculate: sqrt(ux*ux + uy*uy + uz*uz). Magnitude will be stored in _extraQuantities
+ */
+void FieldlinesState::addExtraQuantities(ccmc::Kameleon* kameleon,
+                                         std::vector<std::string>& xtraScalarVars,
+                                         std::vector<std::string>& xtraMagVars) {
+
+    loadExtrasIntoKameleon(kameleon, xtraScalarVars, xtraMagVars);
+
+    const size_t N_XTRA_SCALARS = xtraScalarVars.size();
+    const size_t N_XTRA_MAGNITUDES = xtraMagVars.size() / 3;
+
+    _extraQuantities.resize(N_XTRA_SCALARS + N_XTRA_MAGNITUDES);
+
+    std::unique_ptr<ccmc::Interpolator> interpolator =
+        std::make_unique<ccmc::KameleonInterpolator>(kameleon->model);
+
+    // ------ Extract all the extraQuantities from kameleon and store in state! ------ //
+    for (const glm::vec3& P : _vertexPositions) {
+        // Load the scalars!
+        for (size_t i = 0; i < N_XTRA_SCALARS; i++) {
+            float val;
+            if (xtraScalarVars[i] == T_AS_P_OVER_RHO) {
+                val = interpolator->interpolate("p", P.x, P.y, P.z);
+                val *= TO_KELVIN;
+                val /= interpolator->interpolate("rho", P.x, P.y, P.z);
+            } else {
+                val = interpolator->interpolate(xtraScalarVars[i], P.x, P.y, P.z);
+
+                // When measuring density in ENLIL CCMC multiply by the radius^2
+                if (xtraScalarVars[i] == "rho" && _model == fls::Model::ENLIL) {
+                    val *= std::pow(P.x * fls::A_U_TO_METER, 2.0f);
+                }
+            }
+            _extraQuantities[i].push_back(val);
+        }
+        // Calculate and store the magnitudes!
+        for (size_t i = 0; i < N_XTRA_MAGNITUDES; ++i) {
+            const size_t IDX = i*3;
+
+            const float X = interpolator->interpolate(xtraMagVars[IDX]  , P.x, P.y, P.z);
+            const float Y = interpolator->interpolate(xtraMagVars[IDX+1], P.x, P.y, P.z);
+            const float Z = interpolator->interpolate(xtraMagVars[IDX+2], P.x, P.y, P.z);
+            float val;
+            // When looking at the current's magnitude in Batsrus, CCMC staff are
+            // only interested in the magnitude parallel to the magnetic field
+            if (_extraQuantityNames[N_XTRA_SCALARS + i] == J_PARALLEL_B) {
+                const glm::vec3 NORM_MAGNETIC =  glm::normalize(glm::vec3(
+                        interpolator->interpolate("bx", P.x, P.y, P.z),
+                        interpolator->interpolate("by", P.x, P.y, P.z),
+                        interpolator->interpolate("bz", P.x, P.y, P.z)));
+                // Magnitude of the part of the current vector that's parallel to
+                // the magnetic field vector!
+                val = glm::dot(glm::vec3(X,Y,Z), NORM_MAGNETIC);
+
+            } else {
+                val = std::sqrt(X*X + Y*Y + Z*Z);
+            }
+            _extraQuantities[i + N_XTRA_SCALARS].push_back(val);
+        }
+    }
+}
+#endif // OPENSPACE_MODULE_KAMELEON_ENABLED
+
 #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED
 /**
 * Converts all glm::vec3 in _vertexPositions from spherical (radius, latitude, longitude)