summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--fourier.cpp15
-rw-r--r--fourier_integrator.cpp81
-rw-r--r--log-fourier.cpp14
-rw-r--r--log-fourier_integrator.cpp183
4 files changed, 189 insertions, 104 deletions
diff --git a/fourier.cpp b/fourier.cpp
index 95d0025..3821623 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -84,11 +84,20 @@ std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Rea
Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ) {
Real e = 0;
- for (unsigned i = 0; i < C.size() / 2; i++) {
- e += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
+ for (unsigned n = 0; n < C.size() / 4 -1; n++) {
+ Real h₂ₙ = M_PI * Δτ;
+ Real h₂ₙ₊₁ = M_PI * Δτ;
+ Real f₂ₙ = R[2*n] * df(λ, p, s, C[2*n]);
+ Real f₂ₙ₊₁ = R[2*n+1] * df(λ, p, s, C[2*n+1]);
+ Real f₂ₙ₊₂ = R[2*n+2] * df(λ, p, s, C[2*n+2]);
+ e += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
}
- return e;
+ return y * e;
}
std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) {
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 4f4eff5..96fbffd 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -8,9 +8,9 @@ int main(int argc, char* argv[]) {
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
- Real y₀ = 0;
- Real yₘₐₓ = 0.5;
- Real Δy = 0.05;
+ Real β₀ = 0;
+ Real βₘₐₓ = 0.5;
+ Real Δβ = 0.05;
unsigned log2n = 8;
Real τₘₐₓ = 20;
@@ -41,13 +41,13 @@ int main(int argc, char* argv[]) {
τ₀ = atof(optarg);
break;
case '0':
- y₀ = atof(optarg);
+ β₀ = atof(optarg);
break;
case 'y':
- yₘₐₓ = atof(optarg);
+ βₘₐₓ = atof(optarg);
break;
case 'd':
- Δy = atof(optarg);
+ Δβ = atof(optarg);
break;
case 'I':
maxIterations = (unsigned)atof(optarg);
@@ -68,8 +68,11 @@ int main(int argc, char* argv[]) {
Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
- Real z = 0.5;
- Real Γ₀ = 1 + τ₀ / 2;
+ Real Γ₀ = 1.0;
+ Real μ = Γ₀;
+ if (τ₀ > 0) {
+ μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+ }
std::vector<Real> C(2 * n);
std::vector<Real> R(2 * n);
@@ -78,42 +81,49 @@ int main(int argc, char* argv[]) {
std::vector<Complex> Ct;
std::vector<Complex> Rt;
- Real y = y₀;
+ Real β = β₀;
if (!loadData) {
// start from the exact solution for τ₀ = 0
for (unsigned i = 0; i < n; i++) {
Real τ = i * Δτ * M_PI;
if (τ₀ > 0) {
- C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+ C[i] = Γ₀ * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
} else {
- C[i] = Γ₀ / 2 * exp(-z * τ) / z;
+ C[i] = Γ₀ * exp(-μ * τ) / μ;
}
if (i > 0) {
C[2 * n - i] = C[i];
}
- R[i] = exp(-z * τ);
+ R[i] = exp(-μ * τ);
}
Ct = fft.fourier(C);
Rt = fft.fourier(R);
} else {
- std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+ std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
cfile.close();
for (unsigned i = 1; i < n; i++) {
C[2 * n - i] = C[i];
}
- std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+ std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
rfile.close();
Ct = fft.fourier(C);
Rt = fft.fourier(R);
- z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
+ μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);
}
- while (y += Δy, y <= yₘₐₓ) {
+ std::vector<Real> Cb = C;
+ std::vector<Real> Rb = R;
+ std::vector<Complex> Ctb = Ct;
+ std::vector<Complex> Rtb = Rt;
+ Real μb = μ;
+
+ Real fac = 1;
+ while (β <= βₘₐₓ) {
Real ΔC = 1;
Real ΔCprev = 1000;
unsigned it = 0;
@@ -123,12 +133,12 @@ int main(int argc, char* argv[]) {
for (unsigned i = 0; i < Rt.size(); i++) {
Real ω = i * Δω;
- Rt[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
+ Rt[i] = (1.0 + pow(β, 2) * RddfCt[i] * Rt[i]) / (μ + 1i * ω);
}
for (unsigned i = 0; i < Ct.size(); i++) {
Real ω = i * Δω;
- Ct[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
+ Ct[i] = (2 * Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(β, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (μ + 1i * ω);
}
std::vector<Real> Cnew = fft.inverse(Ct);
@@ -137,6 +147,8 @@ int main(int argc, char* argv[]) {
Rnew[i] = 0;
}
+ μ *= pow(tanh(C[0]-1)+1, 0.05);
+
ΔC = 0;
for (unsigned i = 0; i < Cnew.size() / 2; i++) {
ΔC += pow(Cnew[i] - C[i], 2);
@@ -151,8 +163,6 @@ int main(int argc, char* argv[]) {
R[i] += γ * (Rnew[i] - R[i]);
}
- z *= Cnew[0];
-
if (it % maxIterations == 0) {
if (ΔC < ΔCprev) {
γ = std::min(1.1 * γ, 1.0);
@@ -163,21 +173,40 @@ int main(int argc, char* argv[]) {
ΔCprev = ΔC;
}
- std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
-
+ std::cerr << "\x1b[2K" << "\r";
+ std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];
}
- Real e = energy(C, R, p, s, λ, y, Δτ);
+ if (std::isnan(C[0])) {
+ Δβ /= 2;
+ β -= Δβ;
+ C = Cb;
+ R = Rb;
+ Ct = Ctb;
+ Rt = Rtb;
+ μ = μb;
+ } else {
+
+ std::cerr << "\x1b[2K" << "\r";
- std::cerr << "y " << y << " " << e << " " << z << std::endl;
+ Real e = energy(C, R, p, s, λ, β, Δτ);
- std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+ std::cerr << "y " << β << " " << e << " " << μ << std::endl;
+
+ std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real));
outfile.close();
- std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+ std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
outfileR.close();
+ β += Δβ;
+ Cb = C;
+ Rb = R;
+ Rtb = Rt;
+ Ctb = Ct;
+ μb = μ;
+ }
}
return 0;
diff --git a/log-fourier.cpp b/log-fourier.cpp
index 1e2e7d9..7461a70 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -1,4 +1,5 @@
#include "log-fourier.hpp"
+#include <complex>
LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) {
τₛ = -0.5 * N;
@@ -40,7 +41,7 @@ Real LogarithmicFourierTransform::ν(unsigned n) const {
return exp(ω(n));
}
-Complex gamma(Complex z) {
+Complex Γ(Complex z) {
gsl_sf_result logΓ;
gsl_sf_result argΓ;
@@ -52,6 +53,7 @@ Complex gamma(Complex z) {
std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) {
std::vector<Complex> ĉ(N);
std::vector<Real> σs = {1};
+ /* c is either even or zero for negative arguments */
if (symmetric){
σs.push_back(-1);
}
@@ -65,7 +67,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
}
fftw_execute(a_to_â);
for (unsigned n = 0; n < pad*N; n++) {
- â[(pad*N / 2 + n) % (pad*N)] *= pow(1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n));
+ â[(pad*N / 2 + n) % (pad*N)] *= pow(1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n));
}
fftw_execute(â_to_a);
for (unsigned n = 0; n < N; n++) {
@@ -82,14 +84,18 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
for (Real σ : σs) {
for (unsigned n = 0; n < pad * N; n++) {
if (n < N) {
- a[n] = ĉ[n] * exp((1 - k) * ω(n));
+ if (σ < 0) {
+ a[n] = std::conj(ĉ[n]) * exp((1 - k) * ω(n));
+ } else {
+ a[n] = ĉ[n] * exp((1 - k) * ω(n));
+ }
} else {
a[n] = 0;
}
}
fftw_execute(a_to_â);
for (unsigned n = 0; n < pad*N; n++) {
- â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n));
+ â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n));
}
fftw_execute(â_to_a);
for (unsigned n = 0; n < N; n++) {
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 061cb08..b21f275 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -4,15 +4,18 @@
#include <iostream>
int main(int argc, char* argv[]) {
+ /* Model parameters */
unsigned p = 2;
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
+ /* Log-Fourier parameters */
unsigned log2n = 8;
Real Δτ = 0.1;
Real k = 0.1;
+ /* Iteration parameters */
Real ε = 1e-13;
Real γ = 1;
Real βₘₐₓ = 0.7;
@@ -58,96 +61,134 @@ int main(int argc, char* argv[]) {
LogarithmicFourierTransform fft(N, k, Δτ, 4);
- Real Γ₀ = 1.0;
- Real μ = Γ₀;
+ Real Γ₀ = 1.0 + τ₀;
+ Real μ = 1.0;
+/* Real μ = Γ₀;
if (τ₀ > 0) {
μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
}
+*/
- std::vector<Real> C(N);
- std::vector<Real> R(N);
- std::vector<Complex> Ct(N);
- std::vector<Complex> Rt(N);
+ std::vector<Real> Cₜ₋₁(N);
+ std::vector<Real> Rₜ₋₁(N);
+ std::vector<Complex> Ĉₜ₋₁(N);
+ std::vector<Complex> Ȓₜ₋₁(N);
- // start from the exact solution for β = 0
+ /* Start from the exact solution for β = 0 */
for (unsigned n = 0; n < N; n++) {
- if (τ₀ > 0) {
- C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+ if (τ₀ != 1) {
+ Cₜ₋₁[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
} else {
- C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+ Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n));
}
- R[n] = exp(-μ * fft.t(n));
+ Rₜ₋₁[n] = exp(-μ * fft.t(n));
- Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
- Rt[n] = 1.0 / (μ + 1i * fft.ν(n));
+ Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+ Ȓₜ₋₁[n] = 1.0 / (μ + 1i * fft.ν(n));
}
+ std::vector<Real> Cₜ = Cₜ₋₁;
+ std::vector<Real> Rₜ = Rₜ₋₁;
+ std::vector<Complex> Ĉₜ = Ĉₜ₋₁;
+ std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
+
+ Real fac = 1;
Real β = 0;
while (β < βₘₐₓ) {
- Real ΔC = 100;
- while (ΔC > ε) {
- std::vector<Real> RddfC(N);
- std::vector<Real> dfC(N);
- for (unsigned n = 0; n < N; n++) {
- RddfC[n] = R[n] * ddf(λ, p, s, C[n]);
- dfC[n] = df(λ, p, s, C[n]);
- }
- std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
- std::vector<Complex> dfCt = fft.fourier(dfC, true);
-
- std::vector<Complex> Rtnew(N);
- std::vector<Complex> Ctnew(N);
-
- for (unsigned n = 0; n < N; n++) {
- Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n));
- Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n));
-// Ctnew[n] = - 2 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
- }
-
- std::vector<Real> Cnew = fft.inverse(Ctnew);
- std::vector<Real> Rnew = fft.inverse(Rtnew);
-
- ΔC = 0;
- for (unsigned i = 0; i < N; i++) {
- ΔC += std::norm(Ct[i] - Ctnew[i]);
- ΔC += std::norm(Rt[i] - Rtnew[i]);
+ Real ΔC = 100;
+ Real ΔC₀ = 100;
+ unsigned it = 0;
+ while (ΔC > ε) {
+ std::vector<Real> dfC(N);
+ std::vector<Real> RddfC(N);
+ for (unsigned n = 0; n < N; n++) {
+ RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+ dfC[n] = df(λ, p, s, Cₜ[n]);
+ }
+ std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+ std::vector<Complex> dfCt = fft.fourier(dfC, true);
+
+ std::vector<Complex> Ȓₜ₊₁(N);
+ std::vector<Complex> Ĉₜ₊₁(N);
+
+ for (unsigned n = 0; n < N; n++) {
+ Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+ }
+
+ std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+
+ for (unsigned n = 0; n < N; n++) {
+ RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
+ }
+ RddfCt = fft.fourier(RddfC, false);
+ for (unsigned n = 0; n < N; n++) {
+ Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
+ }
+ std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+
+ μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
+
+ ΔC = 0;
+ for (unsigned i = 0; i < N; i++) {
+ ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
+ ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
+ }
+ ΔC = sqrt(ΔC) / (2*N);
+
+ for (unsigned i = 0; i < N; i++) {
+ Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+ Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+ Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+ Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
+ }
+
+ /*
+ if (ΔC < ΔC₀) {
+ ΔC₀ = ΔC;
+ it = 0;
+ γ = std::min(1.001 * γ, 1.0);
+ } else {
+ it++;
+ }
+
+ if (it > 100) {
+ γ = std::max(0.5 * γ, 1e-3);
+ it = 0;
+ ΔC₀ = ΔC;
+ }
+ */
+
+ std::cerr << "\x1b[2K" << "\r";
+ std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
}
- ΔC = sqrt(ΔC) / (2*N);
- for (unsigned i = 0; i < N; i++) {
- C[i] += γ * (Cnew[i] - C[i]);
- R[i] += γ * (Rnew[i] - R[i]);
- Ct[i] += γ * (Ctnew[i] - Ct[i]);
- Rt[i] += γ * (Rtnew[i] - Rt[i]);
+ /* Integrate the energy using Simpson's rule */
+ Real E = 0;
+ for (unsigned n = 0; n < N/2-1; n++) {
+ Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n);
+ Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+ Real f₂ₙ = Rₜ[2*n] * df(λ, p, s, Cₜ[2*n]);
+ Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]);
+ Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]);
+ E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+ (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+ + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+ );
}
-
- μ *= C[0];
-
-// std::cerr << ΔC << std::endl;
- }
-
- /* Integrate the energy using Simpson's rule */
- Real E = 0;
- for (unsigned i = 0; i < N/2-1; i++) {
- Real h₂ᵢ = fft.t(2*i+1) - fft.t(2*i);
- Real h₂ᵢ₊₁ = fft.t(2*i+2) - fft.t(2*i+1);
- Real f₂ᵢ = R[2*i] * df(λ, p, s, C[2*i]);
- Real f₂ᵢ₊₁ = R[2*i+1] * df(λ, p, s, C[2*i+1]);
- Real f₂ᵢ₊₂ = R[2*i+2] * df(λ, p, s, C[2*i+2]);
- E += (h₂ᵢ + h₂ᵢ₊₁) / 6 * (
- (2 - h₂ᵢ₊₁ / h₂ᵢ) * f₂ᵢ
- + pow(h₂ᵢ + h₂ᵢ₊₁, 2) / (h₂ᵢ * h₂ᵢ₊₁) * f₂ᵢ₊₁
- + (2 - h₂ᵢ / h₂ᵢ₊₁) * f₂ᵢ₊₂
- );
- }
- E *= β;
-
- std::cerr << β << " " << μ << " " << Ct[0].real() << " " << E << std::endl;
- β += Δβ;
+ E *= β;
+
+ std::cerr << "\x1b[2K" << "\r";
+ std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
+ β += Δβ;
+ Cₜ₋₁ = Cₜ;
+ Rₜ₋₁ = Rₜ;
+ Ĉₜ₋₁ = Ĉₜ;
+ Ȓₜ₋₁ = Ȓₜ;
}
for (unsigned i = 0; i < N; i++) {
- std::cout << fft.t(i) << " " << C[i] << std::endl;
+ std::cout << fft.t(i) << " " << Cₜ[i] << std::endl;
}
return 0;