Skip to content

File defectsobservables.h

File List > code_source > cosmolattice > include > CosmoInterface > definitions > defectsmodule > defectsobservables.h

Go to the documentation of this file

#ifndef COSMOINTERFACE_DEFINITIONS_WINDINGNUMBER_H
#define COSMOINTERFACE_DEFINITIONS_WINDINGNUMBER_H

/* This file is part of CosmoLattice, available at www.cosmolattice.net .
 *   Copyright Daniel G. Figueroa, Adrien Florio, Francisco Torrenti and Wessel Valkenburg.
 *   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Jorge Baeza-Ballesteros,  Year: 2026

#include "TempLat/util/constants.h"

#include "TempLat/util/rangeiteration/tagliteral.h"
#include "TempLat/lattice/algebra/operators/power.h"
#include "TempLat/lattice/algebra/operators/operators.h"
#include "TempLat/lattice/algebra/helpers/doeval.h"

#include "TempLat/util/rangeiteration/make_list_tag.h"
#include "TempLat/util/rangeiteration/sum_in_range.h"

namespace TempLat
{
  class DefectsObservables
  {
  public:
    template <typename Model> static auto computeWindingNumberLengthScalarSinglet(const Model &model)
    {
      // return 2. / 3. * average(WindingNumberLengthScalarSinglet(model.fldS(0_c), model.fldS(1_c)));
      auto f = Complexify(model.fldS(0_c), model.fldS(1_c));
      return 2. / 3. *
             average(abs(arg2(f * conj(shift<1_c>(f))) + arg2(shift<1_c>(f) * conj(shift<2_c>(shift<1_c>(f)))) +
                         arg2(shift<2_c>(shift<1_c>(f)) * conj(shift<2_c>(f))) + arg2(shift<2_c>(f) * conj(f))) +
                     abs(arg2(f * conj(shift<2_c>(f))) + arg2(shift<2_c>(f) * conj(shift<2_c>(shift<3_c>(f)))) +
                         arg2(shift<2_c>(shift<3_c>(f)) * conj(shift<3_c>(f))) + arg2(shift<3_c>(f) * conj(f))) +
                     abs(arg2(f * conj(shift<3_c>(f))) + arg2(shift<3_c>(f) * conj(shift<3_c>(shift<1_c>(f)))) +
                         arg2(shift<3_c>(shift<1_c>(f)) * conj(shift<1_c>(f))) + arg2(shift<1_c>(f) * conj(f)))) /
             static_cast<typename Model::FloatType>(2.) / Constants::pi<typename Model::FloatType> *
             pow<3>(GetNGrid::get(model)) * model.dx;
    }

    template <typename Model> static auto computeWindingNumberLengthComplexScalarU1(const Model &model)
    {
      const auto &f = model.fldCS(0_c);
      const auto &a1 = model.fldU1(0_c)(1_c) * model.dx * model.gQ_CsU1(0_c, 0_c);
      const auto &a2 = model.fldU1(0_c)(2_c) * model.dx * model.gQ_CsU1(0_c, 0_c);
      const auto &a3 = model.fldU1(0_c)(3_c) * model.dx * model.gQ_CsU1(0_c, 0_c);
      const auto &expa1 = complexPhase(model.fldU1(0_c)(1_c) * model.dx * model.gQ_CsU1(0_c, 0_c));
      const auto &expa2 = complexPhase(model.fldU1(0_c)(2_c) * model.dx * model.gQ_CsU1(0_c, 0_c));
      const auto &expa3 = complexPhase(model.fldU1(0_c)(3_c) * model.dx * model.gQ_CsU1(0_c, 0_c));

      return 2. / 3. *
             average(abs(arg2(f * expa1 * conj(shift<1_c>(f))) - a1 +
                         arg2(shift<1_c>(f * expa2) * conj(shift<2_c>(shift<1_c>(f)))) - shift<1_c>(a2) +
                         arg2(shift<2_c>(shift<1_c>(f)) * conj(shift<2_c>(f * expa1))) + shift<2_c>(a1) +
                         arg2(shift<2_c>(f) * conj(f * expa2)) + a2) +
                     abs(arg2(f * expa2 * conj(shift<2_c>(f))) - a2 +
                         arg2(shift<2_c>(f * expa3) * conj(shift<2_c>(shift<3_c>(f)))) - shift<2_c>(a3) +
                         arg2(shift<2_c>(shift<3_c>(f)) * conj(shift<3_c>(f * expa2))) + shift<3_c>(a2) +
                         arg2(shift<3_c>(f) * conj(f * expa3)) + a3) +
                     abs(arg2(f * expa3 * conj(shift<3_c>(f))) - a3 +
                         arg2(shift<3_c>(f * expa1) * conj(shift<3_c>(shift<1_c>(f)))) - shift<3_c>(a1) +
                         arg2(shift<3_c>(shift<1_c>(f)) * conj(shift<1_c>(f * expa3))) + shift<1_c>(a3) +
                         arg2(shift<1_c>(f) * conj(f * expa1)) + a1)) /
             static_cast<typename Model::FloatType>(2.) / Constants::pi<typename Model::FloatType> *
             pow<3>(GetNGrid::get(model)) * model.dx;
    }

    template <typename Model> static auto computeAreaParameter(const Model &model)
    {
      const auto &f = model.fldS(0_c);
      return model.dx * model.dx *
             average(Total(i, 1_c, 3_c, heaviside(-f * shift<i>(f))) * abs(Total(j, 1_c, 3_c, shift<j>(f) - f)) /
                     Total(j, 1_c, 3_c, abs(shift<j>(f) - f))) *
             pow<3>(GetNGrid::get(model));
    }
  };

} // namespace TempLat

#endif