diff --git a/.gitignore b/.gitignore
index 831406f..de33913 100644
--- a/.gitignore
+++ b/.gitignore
@@ -26,3 +26,8 @@ perf.*
 tmp/
 
 .vscode/
+
+*.s
+*.ii
+*.yaml
+*.i
diff --git a/Few_body_using_ERK/ERK.hpp b/Few_body_using_ERK/ERK.hpp
index 0845132..06dfae9 100644
--- a/Few_body_using_ERK/ERK.hpp
+++ b/Few_body_using_ERK/ERK.hpp
@@ -12,7 +12,7 @@
 
 namespace mp = boost::multiprecision; // for sqrt
 
-using float_mp = mp::number<mp::cpp_dec_float<50>>;
+using float_mp = double; // mp::number<mp::cpp_dec_float<50>>;
 
 namespace mino2357{
 
diff --git a/Few_body_using_ERK/main.cpp b/Few_body_using_ERK/main.cpp
index f8e4de9..e9683b3 100644
--- a/Few_body_using_ERK/main.cpp
+++ b/Few_body_using_ERK/main.cpp
@@ -44,22 +44,25 @@ int main() {
 	auto intv = float_mp("1.0e-5");
 	auto write_time = float_mp("0.1") * mino2357::min_dt<float_mp>();
 	auto time = float_mp("0.0");
-	auto end_time = float_mp("1.0e2");
+	auto end_time = float_mp("1.0e-3");
 	mino2357::RKF78<float_mp> rkf78(tol, tol);
 
 	for(;;){
 		rkf78.Integrate(time, dt, x);
-		if(time > write_time){
+		//if(time > write_time){
+			/*
 			std::cerr << time << " " << angular_momentum(x) << std::endl;
 		
 			for(size_t i=0; i<num/2; i++){
 				std::cout << x.vec[i] << " ";
 			}
 			std::cout << std::endl;
-			write_time = time + intv;
-		}
+			*/
+			//write_time = time + intv;
+		//}
 		if(time > end_time) {
 			break;
 		}
 	}
+	std::cout << "finished" << std::endl;
 }
diff --git a/Few_body_using_ERK/makefile b/Few_body_using_ERK/makefile
index 9f34a18..f3cf195 100644
--- a/Few_body_using_ERK/makefile
+++ b/Few_body_using_ERK/makefile
@@ -1,6 +1,6 @@
 CXX = /usr/bin/g++
 TARGET = main
-CXXFLAGS = -Wall -Wextra -O2 -march=native -mtune=native -std=c++2b 
+CXXFLAGS = -S -save-temps -Wall -Wextra -O2 -march=native -mtune=native -std=c++2b 
 LDLFLAGS = -lstdc++
 SRCS = main.cpp
 OBJS = $(SRCS:.cpp=.o)
@@ -11,4 +11,4 @@ run: all
 	./$(TARGET)
 
 clean:
-	rm $(TARGET)
+	rm $(TARGET)
\ No newline at end of file
diff --git a/Simple_Solar_System_double_precision/main2 b/Simple_Solar_System_double_precision/main2
deleted file mode 100755
index 36459ea..0000000
Binary files a/Simple_Solar_System_double_precision/main2 and /dev/null differ
diff --git a/Simple_Solar_System_double_precision/main2.cpp b/Simple_Solar_System_double_precision/main2.cpp
index 0db12c8..0b6b870 100644
--- a/Simple_Solar_System_double_precision/main2.cpp
+++ b/Simple_Solar_System_double_precision/main2.cpp
@@ -10,7 +10,7 @@
 //  Sun     |  1.0       |    0.0
 //  Mercury |  1.660e-7  |    0.387
 //  Venus   |  2.447e-6  |    0.723
-//  Eeath   |  3.003e-6  |    1.000
+//  Earth   |  3.003e-6  |    1.000
 //  Mars    |  3.226e-7  |    1.524
 //  Jupiter |  9.543e-4  |    5.204
 //  Saturn  |  2.857e-4  |    9.582
@@ -46,7 +46,7 @@ namespace mino2357 {
             ret.vec[i+1] = float_mp(0.0);
             for(size_t j=num/2; j<num; j+=2){
                 if(j!=i){
-		    auto m_j = m[(j-num/2)/2];
+				    auto m_j = m[(j-num/2)/2];
                     auto x = u.vec[i-num/2] - u.vec[j-num/2];
                     auto y = u.vec[i+1-num/2] - u.vec[j+1-num/2];;
                     auto r = sqrt(x*x + y*y);
@@ -62,7 +62,7 @@ namespace mino2357 {
 
 template <typename T>
 inline std::tuple<T, T> eccentricity_vector(const T& pos_x, const T& pos_y, const T& vel_x, const T& vel_y) noexcept {
-        T r = sqrt(pos_x * pos_x + pos_y * pos_y);
+	T r = sqrt(pos_x * pos_x + pos_y * pos_y);
 	T v = sqrt(vel_x * vel_x + vel_y * vel_y);
 	T mu = T(1.0);
 	T rv = pos_x * vel_x + pos_y * vel_y;
@@ -140,6 +140,7 @@ int main() {
 
 	for(size_t i=0;;i++){
 		rkf78.Integrate(mino2357::time, mino2357::dt, x);
+		/*
 		if(i%mino2357::intv == 0){
 			// position
 			std::cout << mino2357::time / (float_mp(2.0) * float_mp(3.141592653));
@@ -148,13 +149,13 @@ int main() {
 			}
 			std::cout << std::endl;
 			// velocity
-			/*
-			std::cout << mino2357::time / (float_mp(2.0) * float_mp(3.141592653));
-			for(size_t i=mino2357::num/2; i<mino2357::num; i+=2){
-				std::cout << " " << x.vec[i] - x.vec[0] << " " << x.vec[i+1] - x.vec[1];
-			}
-			std::cout << std::endl;
-			*/
+			
+			//std::cout << mino2357::time / (float_mp(2.0) * float_mp(3.141592653));
+			//for(size_t i=mino2357::num/2; i<mino2357::num; i+=2){
+			//	std::cout << " " << x.vec[i] - x.vec[0] << " " << x.vec[i+1] - x.vec[1];
+			//}
+			//std::cout << std::endl;
+			
 			// eccentricity_vector
 			std::cerr << mino2357::time / (float_mp(2.0) * float_mp(3.141592653));
 			for(size_t i=2; i<mino2357::num/2; i+=2){
@@ -169,5 +170,9 @@ int main() {
 
 			mino2357::write_time = mino2357::time + mino2357::intv;
 		}
+		*/
+		if(mino2357::time > 1.0){
+			break;
+		}
 	}
 }
diff --git a/Simple_Solar_System_double_precision_MERCURY/main2 b/Simple_Solar_System_double_precision_MERCURY/main2
deleted file mode 100755
index 8850315..0000000
Binary files a/Simple_Solar_System_double_precision_MERCURY/main2 and /dev/null differ
diff --git a/Tsuchinshan-ATLAS-2024-10-20/.gitignore b/Tsuchinshan-ATLAS-2024-10-20/.gitignore
new file mode 100644
index 0000000..fa52ec2
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/.gitignore
@@ -0,0 +1,5 @@
+*.log
+*.out
+main2
+log_data/
+exp/
diff --git a/Tsuchinshan-ATLAS-2024-10-20/ERK.hpp b/Tsuchinshan-ATLAS-2024-10-20/ERK.hpp
new file mode 100644
index 0000000..42f005e
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/ERK.hpp
@@ -0,0 +1,291 @@
+#pragma once
+
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <iomanip>
+#include <string>
+#include <random>
+#include "vector.hpp"
+#include "parameter.hpp"
+
+namespace mino2357{
+    template <typename T>
+    constexpr T ratio(int a, int b){
+        return static_cast<T>(a) / static_cast<T>(b);
+    }
+
+    template <typename T>
+    constexpr T alpha(){
+        return static_cast<T>(0.8);
+    }
+    
+    template <typename T>
+    class ButcherRKF78{
+    public:
+        T table[13][13];
+        T order8[13];
+        T order9[13];
+
+        constexpr ButcherRKF78();
+        
+        constexpr T operator()(int i, int j){
+            return table[i][j];
+        };
+
+        constexpr T o8(int i){
+            return order8[i];
+        };
+        
+        constexpr T o9(int i){
+            return order9[i];
+        };
+    };
+
+    template <typename T>
+    constexpr ButcherRKF78<T>::ButcherRKF78(){
+        T zero = 0;
+        
+        table[0][0] = zero;
+
+        table[1][0] = ratio<T>(2, 27);
+        table[0][1] = zero;
+        
+        table[2][0] = ratio<T>(1, 36);
+        table[2][1] = ratio<T>(1, 12);
+        table[2][2] = zero;
+        
+        table[3][0] = ratio<T>(1, 24);
+        table[3][1] = zero;
+        table[3][2] = ratio<T>(1, 8);
+        table[3][3] = zero;
+        
+        table[4][0] = ratio<T>(5, 12);
+        table[4][1] = zero;
+        table[4][2] = ratio<T>(-25, 16);
+        table[4][3] = ratio<T>(25, 16);
+        table[4][4] = zero;
+        
+        table[5][0] = ratio<T>(1, 20);
+        table[5][1] = zero;
+        table[5][2] = zero;
+        table[5][3] = ratio<T>(1, 4);
+        table[5][4] = ratio<T>(1, 5);
+        table[5][5] = zero;
+        
+        table[6][0] = ratio<T>(-25, 108);
+        table[6][1] = zero;
+        table[6][2] = zero;
+        table[6][3] = ratio<T>(125, 108);
+        table[6][4] = ratio<T>(-65, 27);
+        table[6][5] = ratio<T>(125, 54);
+        table[6][6] = zero;
+        
+        table[7][0] = ratio<T>(31, 300);
+        table[7][1] = zero;
+        table[7][2] = zero;
+        table[7][3] = zero;
+        table[7][4] = ratio<T>(61, 225);
+        table[7][5] = ratio<T>(-2, 9);
+        table[7][6] = ratio<T>(13, 900);
+        table[7][7] = zero;
+        
+        table[8][0] = ratio<T>(2, 1);
+        table[8][1] = zero;
+        table[8][2] = zero;
+        table[8][3] = ratio<T>(-53, 6);
+        table[8][4] = ratio<T>(704, 45);
+        table[8][5] = ratio<T>(-107, 9);
+        table[8][6] = ratio<T>(67, 90);
+        table[8][7] = ratio<T>(3, 1);
+        table[8][8] = zero;
+        
+        table[9][0] = ratio<T>(-91, 108);
+        table[9][1] = zero;
+        table[9][2] = zero;
+        table[9][3] = ratio<T>(23, 108);
+        table[9][4] = ratio<T>(-976, 135);
+        table[9][5] = ratio<T>(311, 54);
+        table[9][6] = ratio<T>(-19, 60);
+        table[9][7] = ratio<T>(17, 6);
+        table[9][8] = ratio<T>(-1, 12);
+        table[9][9] = zero;
+        
+        table[10][0]  = ratio<T>(2383, 4100);
+        table[10][1]  = zero;
+        table[10][2]  = zero;
+        table[10][3]  = ratio<T>(-341, 164);
+        table[10][4]  = ratio<T>(4496, 1025);
+        table[10][5]  = ratio<T>(-301, 82);
+        table[10][6]  = ratio<T>(2133, 4100);
+        table[10][7]  = ratio<T>(45, 82);
+        table[10][8]  = ratio<T>(45, 164);
+        table[10][9]  = ratio<T>(18, 41);
+        table[10][10] = zero;
+        
+        table[11][0]  = ratio<T>(3, 205);
+        table[11][1]  = zero;
+        table[11][2]  = zero;
+        table[11][3]  = zero;
+        table[11][4]  = zero;
+        table[11][5]  = ratio<T>(-6, 41);
+        table[11][6]  = ratio<T>(-3, 205);
+        table[11][7]  = ratio<T>(-3, 41);
+        table[11][8]  = ratio<T>(3, 41);
+        table[11][9]  = ratio<T>(6, 41);
+        table[11][10] = zero;
+        table[11][11] = zero;
+        
+        table[12][0]  = ratio<T>(-1777, 4100);
+        table[12][1]  = zero;
+        table[12][2]  = zero;
+        table[12][3]  = ratio<T>(-341, 164);
+        table[12][4]  = ratio<T>(4496, 1025);
+        table[12][5]  = ratio<T>(-289, 82);
+        table[12][6]  = ratio<T>(2193, 4100);
+        table[12][7]  = ratio<T>(51, 82);
+        table[12][8]  = ratio<T>(33, 164);
+        table[12][9]  = ratio<T>(12, 41);
+        table[12][10] = zero;
+        table[12][11] = ratio<T>(1, 1);
+        table[12][12] = zero;
+
+        order8[0]  = ratio<T>(41, 840);
+        order8[1]  = zero;
+        order8[2]  = zero;
+        order8[3]  = zero;
+        order8[4]  = zero;
+        order8[5]  = ratio<T>(34, 105);
+        order8[6]  = ratio<T>(9, 35);
+        order8[7]  = ratio<T>(9, 35);
+        order8[8]  = ratio<T>(9, 280);
+        order8[9]  = ratio<T>(9, 280);
+        order8[10] = ratio<T>(41, 840);
+        order8[11] = zero;
+        order8[12] = zero;
+        
+        order9[0]  = zero;
+        order9[1]  = zero;
+        order9[2]  = zero;
+        order9[3]  = zero;
+        order9[4]  = zero;
+        order9[5]  = ratio<T>(34, 105);
+        order9[6]  = ratio<T>(9, 35);
+        order9[7]  = ratio<T>(9, 35);
+        order9[8]  = ratio<T>(9, 280);
+        order9[9]  = ratio<T>(9, 280);
+        order9[10] = zero;
+        order9[11] = ratio<T>(41, 840);
+        order9[12] = ratio<T>(41, 840);
+    }
+
+    template <typename T>
+    class RKF78{
+    public:
+        T crt_h;
+        T next_h;
+        T A_Tol;
+        T R_Tol;//今は使わない
+
+        constexpr RKF78(T, T);
+        
+        inline constexpr void Integrate(T& t, T& dt, mino2357::vector<T>& x) noexcept;
+    };
+
+    template <typename T>
+    constexpr RKF78<T>::RKF78(T at,T rt){
+        A_Tol = at;
+        R_Tol = rt;
+    }
+
+    template <typename T>
+    inline constexpr void RKF78<T>::Integrate(T& t, T& dt, mino2357::vector<T>& x) noexcept{
+        crt_h = dt;
+
+        mino2357::vector<T> x9(x.vec.size()), x8(x.vec.size()), temp(x.vec.size());
+        T R, delta;
+        R = 0;
+
+        mino2357::vector<T> k0(x.vec.size()), k1(x.vec.size()), k2(x.vec.size()), k3(x.vec.size()), k4(x.vec.size()), k5(x.vec.size()), k6(x.vec.size()), k7(x.vec.size()), k8(x.vec.size()), k9(x.vec.size()), k10(x.vec.size()), k11(x.vec.size()), k12(x.vec.size());
+
+        ButcherRKF78<T> f;
+
+        k0 = func<T>(x);
+        k1 = func<T>(x + crt_h * f(1, 0) * k0);
+        k2 = func<T>(x + crt_h * f(2, 0) * k0 + crt_h * f(2, 1) * k1);
+        k3 = func<T>(x + crt_h * f(3, 0) * k0 + crt_h * f(3, 1) * k1 + crt_h * f(3, 2) * k2);
+        k4 = func<T>(x + crt_h * f(4, 0) * k0 + crt_h * f(4, 1) * k1 + crt_h * f(4, 2) * k2 + crt_h * f(4, 3) * k3);
+        k5 = func<T>(x + crt_h * f(5, 0) * k0 + crt_h * f(5, 1) * k1 + crt_h * f(5, 2) * k2 + crt_h * f(5, 3) * k3 + crt_h * f(5, 4) * k4);
+        k6 = func<T>(x + crt_h * f(6, 0) * k0 + crt_h * f(6, 1) * k1 + crt_h * f(6, 2) * k2 + crt_h * f(6, 3) * k3 + crt_h * f(6, 4) * k4
+                       + crt_h * f(6, 5) * k5);
+        k7 = func<T>(x + crt_h * f(7, 0) * k0 + crt_h * f(7, 1) * k1 + crt_h * f(7, 2) * k2 + crt_h * f(7, 3) * k3 + crt_h * f(7, 4) * k4
+                       + crt_h * f(7, 5) * k5 + crt_h * f(7, 6) * k6);
+        k8 = func<T>(x + crt_h * f(8, 0) * k0 + crt_h * f(8, 1) * k1 + crt_h * f(8, 2) * k2 + crt_h * f(8, 3) * k3 + crt_h * f(8, 4) * k4
+                       + crt_h * f(8, 5) * k5 + crt_h * f(8, 6) * k6 + crt_h * f(8, 7) * k7);
+        k9 = func<T>(x + crt_h * f(9, 0) * k0 + crt_h * f(9, 1) * k1 + crt_h * f(9, 2) * k2 + crt_h * f(9, 3) * k3 + crt_h * f(9, 4) * k4
+                       + crt_h * f(9, 5) * k5 + crt_h * f(9, 6) * k6 + crt_h * f(9, 7) * k7 + crt_h * f(9, 8) * k8);
+        k10= func<T>(x + crt_h * f(10, 0) * k0 + crt_h * f(10, 1) * k1 + crt_h * f(10, 2) * k2 + crt_h * f(10, 3) * k3 + crt_h * f(10, 4) * k4
+                       + crt_h * f(10, 5) * k5 + crt_h * f(10, 6) * k6 + crt_h * f(10, 7) * k7 + crt_h * f(10, 8) * k8 + crt_h * f(10, 9) * k9);
+        k11= func<T>(x + crt_h * f(11, 0) * k0 + crt_h * f(11, 1) * k1 + crt_h * f(11, 2) * k2 + crt_h * f(11, 3) * k3 + crt_h * f(11, 4) * k4
+                       + crt_h * f(11, 5) * k5 + crt_h * f(11, 6) * k6 + crt_h * f(11, 7) * k7 + crt_h * f(11, 8) * k8 + crt_h * f(11, 9) * k9
+                       + crt_h * f(11, 10) * k10);
+        k12= func<T>(x + crt_h * f(12, 0) * k0 + crt_h * f(12, 1) * k1 + crt_h * f(12, 2) * k2 + crt_h * f(12, 3) * k3 + crt_h * f(12, 4) * k4
+                       + crt_h * f(12, 5) * k5 + crt_h * f(12, 6) * k6 + crt_h * f(12, 7) * k7 + crt_h * f(12, 8) * k8 + crt_h * f(12, 9) * k9
+                       + crt_h * f(12, 10) * k10 + crt_h * f(12, 11) * k11);
+        
+        x8 = x + crt_h * (f.o8(0) * k0
+                        + f.o8(1) * k1
+                        + f.o8(2) * k2
+                        + f.o8(3) * k3
+                        + f.o8(4) * k4
+                        + f.o8(5) * k5
+                        + f.o8(6) * k6
+                        + f.o8(7) * k7
+                        + f.o8(8) * k8
+                        + f.o8(9) * k9
+                        + f.o8(10) * k10
+                        + f.o8(11) * k11
+                        + f.o8(12) * k12);
+        
+        x9 = x + crt_h * (f.o9(0) * k0
+                        + f.o9(1) * k1
+                        + f.o9(2) * k2
+                        + f.o9(3) * k3
+                        + f.o9(4) * k4
+                        + f.o9(5) * k5
+                        + f.o9(6) * k6
+                        + f.o9(7) * k7
+                        + f.o9(8) * k8
+                        + f.o9(9) * k9
+                        + f.o9(10) * k10
+                        + f.o9(11) * k11
+                        + f.o9(12) * k12);
+
+        temp = x8 - x9; // todo!
+
+        for(size_t i=0; i<temp.vec.size(); ++i){
+            R += temp.vec[i] * temp.vec[i];
+        }
+        R = R / crt_h;
+
+        delta = sqrt(R);
+       
+        if(delta > A_Tol){
+            //std::cerr << Retry  << t <<   << dt << std::endl;
+            dt = crt_h * pow(alpha<T>() * A_Tol / delta, ratio<T>(1, 5));
+            return;
+        }
+
+        x = x9;
+
+        t += crt_h;
+        
+        next_h = crt_h * pow(alpha<T>() * A_Tol / delta, ratio<T>(1, 8));
+
+        dt = next_h;
+
+	if(dt > max_dt){
+		dt = max_dt;
+	}
+    }
+}
diff --git a/Tsuchinshan-ATLAS-2024-10-20/angular_momentum.gp b/Tsuchinshan-ATLAS-2024-10-20/angular_momentum.gp
new file mode 100644
index 0000000..e2c91bc
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/angular_momentum.gp
@@ -0,0 +1,6 @@
+set xlabel "yr"
+set ylabel "angular momentum"
+plot "log.log" u 1:102 w l title "Lx"
+replot "log.log" u 1:103 w l title "Ly"
+replot "log.log" u 1:104 w l title "Lz"
+replot "log.log" u 1:105 w l title "||L||"
diff --git a/Tsuchinshan-ATLAS-2024-10-20/eccentricity.gp b/Tsuchinshan-ATLAS-2024-10-20/eccentricity.gp
new file mode 100644
index 0000000..4311ab4
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/eccentricity.gp
@@ -0,0 +1,15 @@
+# set terminal dumb # size 150, 50
+intv=100000
+start=1
+#set xr [0.9999:1.0001]
+set xlabel "time [Myr]"
+set ylabel "eccentricity"
+plot   "log.log" every intv::start u ($1/1000000):21 w l title "Mercury"
+replot "log.log" every intv::start u ($1/1000000):31 w l title "Venus"
+replot "log.log" every intv::start u ($1/1000000):41 w l title "Earth"
+replot "log.log" every intv::start u ($1/1000000):51 w l title "Mars"
+replot "log.log" every intv::start u ($1/1000000):61 w l title "Jupiter"
+replot "log.log" every intv::start u ($1/1000000):71 w l title "Saturn"
+replot "log.log" every intv::start u ($1/1000000):81 w l title "Uranus"
+replot "log.log" every intv::start u ($1/1000000):91 w l title "Neptune"
+replot "log.log" every intv::start u ($1/1000000):101 w l title "*Pluto"
diff --git a/Tsuchinshan-ATLAS-2024-10-20/eccentricity2.gp b/Tsuchinshan-ATLAS-2024-10-20/eccentricity2.gp
new file mode 100644
index 0000000..5fbd831
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/eccentricity2.gp
@@ -0,0 +1,16 @@
+# set terminal dumb size 100, 50
+intv=100
+start=3
+end=1000000000
+set xr [0:0.2]
+set xlabel "Myr"
+set ylabel "eccentricity"
+#plot   "log.log" every intv::start u 1:21 w l title "Mercury"
+#replot "log.log" every intv::start u 1:31 w l title "Venus"
+#replot "log.log" every intv::start u 1:41 w l title "Earth"
+#replot "log.log" every intv::start u 1:51 w l title "Mars"
+plot "log.log" every intv::start::end u ($1/1000000):61 w p ps 0.01 title "Jupiter"
+replot "log.log" every intv::start::end u ($1/1000000):71 w p ps 0.01 title "Saturn"
+#replot "log.log" every intv::start u 1:81 w l title "Uranus"
+#replot "log.log" every intv::start u 1:91 w l title "Neptune"
+#replot "log.log" every intv::start u 1:101 w l title "*Pluto"
diff --git a/Tsuchinshan-ATLAS-2024-10-20/jpl_2023-09-01.data b/Tsuchinshan-ATLAS-2024-10-20/jpl_2023-09-01.data
new file mode 100644
index 0000000..6a0838f
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/jpl_2023-09-01.data
@@ -0,0 +1,139 @@
+
+x.vec[0] = 0.0 / au;
+x.vec[1] = 0.0 / au;
+x.vec[2] = 0.0 / au;
+x.vec[3] = -1.404774254690487E+07 / au;
+x.vec[4] = -6.825208888976701E+07 / au;
+x.vec[5] = -4.289153234359771E+06 / au;
+x.vec[6] = 4.608792459711022E+07 / au;
+x.vec[7] = -9.859187757812914E+07 / au;
+x.vec[8] = -4.013349395839728E+06 / au;
+x.vec[9] = 1.330220204950462E+08 / au;
+x.vec[10] = 6.704428213911822E+07 / au;
+x.vec[11] = -4.869615619417280E+03 / au;
+x.vec[12] = 6.649302608051282E+07 / au;
+x.vec[13] = 2.186169120955794E+08 / au;
+x.vec[14] = 2.950343934283048E+06 / au;
+x.vec[15] = 2.385041315415040E+08 / au;
+x.vec[16] = 7.182563535822021E+08 / au;
+x.vec[17] = -8.319703965834230E+06 / au;
+x.vec[18] = 1.406365269217089E+09 / au;
+x.vec[19] = -3.235355848531194E+08 / au;
+x.vec[20] = -5.034784972917840E+07 / au;
+x.vec[21] = 1.696652737323582E+09 / au;
+x.vec[22] = 2.385017003636673E+09 / au;
+x.vec[23] = -1.314104590030003E+07 / au;
+x.vec[24] = 4.469410600175806E+09 / au;
+x.vec[25] = -1.293917449213956E+08 / au;
+x.vec[26] = -1.003300431929548E+08 / au;
+x.vec[27] = 9.755484379394056E+07 / au;
+x.vec[28] = -5.078485333345367E+06 / au;
+x.vec[29] = 3.513104848711042E+07 / au;
+
+x.vec[30] = 0.0 / au;
+x.vec[31] = 0.0 / au;
+x.vec[32] = 0.0 / au;
+x.vec[33] = 3.277813191168440E+06 / au / (2.0 * PI) * year;
+x.vec[34] = -6.360278146631564E+05 / au / (2.0 * PI) * year;
+x.vec[35] = -3.526256010095223E+05 / au / (2.0 * PI) * year;
+x.vec[36] = 2.720842478216296E+06 / au / (2.0 * PI) * year;
+x.vec[37] = 1.271078843814805E+06 / au / (2.0 * PI) * year;
+x.vec[38] = -1.395383786068834E+05 / au / (2.0 * PI) * year;
+x.vec[39] = -1.199359489952496E+06 / au / (2.0 * PI) * year;
+x.vec[40] = 2.288045813113915E+06 / au / (2.0 * PI) * year;
+x.vec[41] = -1.848630523101974E+02 / au / (2.0 * PI) * year;
+x.vec[42] = -1.923719964482550E+06 / au / (2.0 * PI) * year;
+x.vec[43] = 7.871516467196010E+05 / au / (2.0 * PI) * year;
+x.vec[44] = 6.367679199693510E+04 / au / (2.0 * PI) * year;
+x.vec[45] = -1.085443281912170E+06 / au / (2.0 * PI) * year;
+x.vec[46] = 4.092938588278859E+05 / au / (2.0 * PI) * year;
+x.vec[47] = 2.258591609536928E+04 / au / (2.0 * PI) * year;
+x.vec[48] = 1.401829990042232E+05 / au / (2.0 * PI) * year;
+x.vec[49] = 8.122858624024559E+05 / au / (2.0 * PI) * year;
+x.vec[50] = -1.976963001796068E+04 / au / (2.0 * PI) * year;
+x.vec[51] = -4.848243855055465E+05 / au / (2.0 * PI) * year;
+x.vec[52] = 3.142563737252478E+05 / au / (2.0 * PI) * year;
+x.vec[53] = 7.456670936695803E+03 / au / (2.0 * PI) * year;
+x.vec[54] = 9.552707840392186E+03 / au / (2.0 * PI) * year;
+x.vec[55] = 4.728243111876977E+05 / au / (2.0 * PI) * year;
+x.vec[56] = -9.988672364090331E+03 / au / (2.0 * PI) * year;
+x.vec[57] = 1.913451980110258E+06 / au / (2.0 * PI) * year;
+x.vec[58] = -2.739124864901621E+06 / au / (2.0 * PI) * year;
+x.vec[59] = 2.814734584202529E+06 / au / (2.0 * PI) * year;
+
+
+	// position
+	x.vec[0]  = 0.0 / au; // sol
+	x.vec[1]  = 0.0 / au;
+	x.vec[2]  = 0.0 / au;
+	x.vec[3]  = 4.749853696679540E+07 / au; // sui X = 4.749853696679540E+07 Y =-3.737724061730301E+07 Z =-7.411266476235436E+06
+	x.vec[4]  = -3.737724061730301E+07 / au;
+	x.vec[5]  = -7.411266476235436E+06 / au;
+	x.vec[6]  = 1.067832390650615E+08 / au; // kinn X = 1.067832390650615E+08 Y =-1.974951870317106E+07 Z =-6.432625455056737E+06
+	x.vec[7]  = -1.974951870317106E+07 / au;
+	x.vec[8]  = -6.432625455056737E+06 / au;
+	x.vec[9]  = 1.400066731258702E+08 / au; // ti X = 1.400066731258702E+08 Y =-5.656864340673850E+07 Z = 2.358117319162935E+03
+	x.vec[10]  = -5.656864340673850E+07 / au;
+	x.vec[11]  = 2.358117319162935E+03 / au;
+	x.vec[12]  = -2.317149312291399E+08 / au; // ka X =-2.317149312291399E+08 Y =-7.301765439010516E+07 Z = 4.153589884319007E+06
+	x.vec[13]  = -7.301765439010516E+07 / au;
+	x.vec[14]  = 4.153589884319007E+06 / au;
+	x.vec[15]  = 6.125858099714541E+08 / au; // moku X = 6.125858099714541E+08 Y = 4.198544407080007E+08 Z =-1.544953710371283E+07
+	x.vec[16]  = 4.198544407080007E+08 / au;
+	x.vec[17]  = -1.544953710371283E+07 / au;
+	x.vec[18]  = 1.309454381823805E+09 / au; // do X = 1.309454381823805E+09 Y =-6.485533471689560E+08 Z =-4.083518400133377E+07
+	x.vec[19]  = -6.485533471689560E+08 / au;
+	x.vec[20]  = -4.083518400133377E+07 / au;
+	x.vec[21]  = 1.891786630777107E+09 / au; // ten X = 1.891786630777107E+09 Y = 2.246534492281295E+09 Z =-1.618145404476953E+07
+	x.vec[22]  = 2.246534492281295E+09 / au;
+	x.vec[23]  = -1.618145404476953E+07 / au;
+	x.vec[24]  = 4.461097261292484E+09 / au; // kai X = 4.461097261292484E+09 Y =-3.255589171394691E+08 Z =-9.610267007640028E+07
+	x.vec[25]  = -3.255589171394691E+08 / au;
+	x.vec[26]  = -9.610267007640028E+07 / au;
+	x.vec[27]  = 2.523199306620363E+09 / au; // mei X = 2.523199306620363E+09 Y =-4.554189299791788E+09 Z =-2.421678207220302E+08
+	x.vec[28]  = -4.554189299791788E+09 / au;
+	x.vec[29]  = -2.421678207220302E+08 / au;
+ 
+	// velocity [km/day]
+	x.vec[30] = 0.0 / au; // sol
+	x.vec[31] = 0.0 / au;
+	x.vec[32] = 0.0 / au;
+	x.vec[33] = 1.772352879533610E+06 / au; // sui VX= 1.772352879533610E+06 VY= 3.503420726740979E+06 VZ= 1.237342409601002E+05
+	x.vec[34] = 3.503420726740979E+06 / au;
+	x.vec[35] = 1.237342409601002E+05 / au;
+	x.vec[36] = 5.355688919820316E+05 / au; // kinn VX= 5.355688919820316E+05 VY= 2.961698410495598E+06 VZ= 9.767588002410719E+03
+	x.vec[37] = 2.961698410495598E+06 / au;
+	x.vec[38] = 9.767588002410719E+03 / au;
+	x.vec[39] = 9.220286795888582E+05 / au; // ti VX= 9.220286795888582E+05 VY= 2.375474361150964E+06 VZ=-2.180159483080331E+02
+	x.vec[40] = 2.375474361150964E+06 / au;
+	x.vec[41] = -2.180159483080331E+02 / au;
+	x.vec[42] = 7.072649115996205E+05 / au; // ka VX= 7.072649115996205E+05 VY=-1.817932890099840E+06 VZ=-5.544867514259787E+04
+	x.vec[43] = -1.817932890099840E+06 / au;
+	x.vec[44] = -5.544867514259787E+04 / au;
+	x.vec[45] = -6.520016889714854E+05 / au; // moku VX=-6.520016889714854E+05 VY= 9.852604412872446E+05 VZ= 1.049049729723986E+04
+	x.vec[46] = 9.852604412872446E+05 / au;
+	x.vec[47] = 1.049049729723986E+04 / au;
+	x.vec[48] = 3.236801726147626E+05 / au; // do VX= 3.236801726147626E+05 VY= 7.473526635600305E+05 VZ=-2.592763479803463E+04
+	x.vec[49] = 7.473526635600305E+05 / au;
+	x.vec[50] = -2.592763479803463E+04 / au;
+	x.vec[51] = -4.550220547965139E+05 / au; // ten VX=-4.550220547965139E+05 VY= 3.526668388808242E+05 VZ= 7.193044943207628E+03
+	x.vec[52] = 3.526668388808242E+05 / au;
+	x.vec[53] = 7.193044943207628E+03 / au;
+	x.vec[54] = 3.057160253033973E+04 / au; // kai VX= 3.057160253033973E+04 VY= 4.723624786308346E+05 VZ=-1.040664801597096E+04
+	x.vec[55] = 4.723624786308346E+05 / au;
+	x.vec[56] = -1.040664801597096E+04 / au;
+	x.vec[57] = 4.215862004260509E+05 / au; // mei VX= 4.215862004260509E+05 VY= 1.243989357680785E+05 VZ=-1.368765412288774E+05
+	x.vec[58] = 1.243989357680785E+05 / au;
+	x.vec[59] = -1.368765412288774E+05 / au;
+
+# Mass
+Sun Mass, 10^23 kg    = 19885000.0
+Mass x10^23 (kg)      =     3.302 sui
+Mass x10^23 (kg)      =    48.685 kinn
+Mass x10^23 (kg)      =     59.7219 ti
+Mass x10^23 (kg)      =     6.4171 ka
+Mass x10^23 (kg)      =    18981.8722 moku
+Mass x10^23 (kg)      =     5683.4 moku
+Mass x10^23 (kg)      =    868.13 do
+Mass x10^23 (kg)      =   1024.09 tenn
+Mass x10^23 (kg)      =     0.1307 kai
\ No newline at end of file
diff --git a/Tsuchinshan-ATLAS-2024-10-20/main2.cpp b/Tsuchinshan-ATLAS-2024-10-20/main2.cpp
new file mode 100644
index 0000000..4b92a56
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/main2.cpp
@@ -0,0 +1,184 @@
+#include <iostream>
+#include <iomanip>
+#include <tuple>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include "vector.hpp"
+#include "parameter.hpp"
+#include "ERK.hpp"
+
+constexpr double PI=3.1415926535897932384626433832795028841971;
+
+const std::vector<float_mp> m = {
+	float_mp(1.0),
+	float_mp(3.302 / 19885000.0),
+	float_mp(48.685 / 19885000.0),
+	float_mp(59.7219 / 19885000.0),
+	float_mp(6.4171 / 19885000.0),
+	float_mp(18981.8722 / 19885000.0),
+	float_mp(5683.4 / 19885000.0),
+	float_mp(868.13 / 19885000.0),
+	float_mp(1024.09 / 19885000.0),
+	float_mp(5.0289163e-18 / 19885000.0), // 10^13
+};
+
+namespace mino2357 {
+    // f: R^N -> R^N
+    template <typename T>
+    inline mino2357::vector<T> func(const mino2357::vector<T>& u) noexcept {
+        auto ret = u;
+        auto num = u.vec.size();
+	for(size_t i=0; i<num/2; ++i){
+            ret.vec[i] = u.vec[i+num/2];
+        }
+        for(size_t i=num/2; i<num; i+=3){
+            ret.vec[i] = float_mp(0.0);
+            ret.vec[i+1] = float_mp(0.0);
+			ret.vec[i+2] = float_mp(0.0);
+            for(size_t j=num/2; j<num; j+=3){
+                if(j!=i){
+				    auto m_j = m[(j-num/2)/3];
+                    auto x = u.vec[i-num/2] - u.vec[j-num/2];
+                    auto y = u.vec[i+1-num/2] - u.vec[j+1-num/2];
+					auto z = u.vec[i+2-num/2] - u.vec[j+2-num/2];
+                    auto r = sqrt(x*x + y*y + z*z);
+                    ret.vec[i]   += - m_j * x / (r * r * r);
+                    ret.vec[i+1] += - m_j * y / (r * r * r);
+		    		ret.vec[i+2] += - m_j * z / (r * r * r);
+                }
+            }
+        }
+        return ret;
+    }
+}
+
+
+template <typename T>
+inline std::tuple<T, T, T> eccentricity_vector(const T& pos_x, const T& pos_y, const T& pos_z, const T& vel_x, const T& vel_y, const T& vel_z) noexcept {
+	T r = sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
+	T v = sqrt(vel_x * vel_x + vel_y * vel_y + vel_z * vel_z);
+	T mu = T(1.0);
+	T rv = pos_x * vel_x + pos_y * vel_y + pos_z * vel_z;
+	T ret_x = (v * v / mu - T(1.0) / r) * pos_x - (rv / mu) * vel_x;
+	T ret_y = (v * v / mu - T(1.0) / r) * pos_y - (rv / mu) * vel_y;
+	T ret_z = (v * v / mu - T(1.0) / r) * pos_z - (rv / mu) * vel_z;
+	return {ret_x, ret_y, ret_z};
+}
+
+template <typename T>
+inline std::tuple<T, T, T> angular_momentum(const mino2357::vector<T>& x) noexcept {
+	auto ret_x = T(0.0);
+	auto ret_y = T(0.0);
+	auto ret_z = T(0.0);
+	for(size_t i=0; i<mino2357::num/2; i+=3){
+		auto r_x = x.vec[i];
+		auto r_y = x.vec[i+1];
+		auto r_z = x.vec[i+2];
+		auto v_x = x.vec[i+mino2357::num/2];
+		auto v_y = x.vec[i+mino2357::num/2+1];
+		auto v_z = x.vec[i+mino2357::num/2+2];
+		auto m_m   = m[i/3];
+		ret_x += r_y * m_m * v_z - r_z * m_m * v_y;
+		ret_y += r_z * m_m * v_x - r_x * m_m * v_z;
+		ret_z += r_x * m_m * v_y - r_y * m_m * v_x;
+	}
+	return {ret_x, ret_y, ret_z};
+}
+
+int main() {
+	std::cout << std::fixed << std::setprecision(std::numeric_limits<float_mp>::digits10 + 1);
+	std::cerr << std::fixed << std::setprecision(std::numeric_limits<float_mp>::digits10 + 1);
+
+	mino2357::vector<float_mp> x(mino2357::num);
+	
+	auto au = 149597870.700; // [km]
+	auto year = 365.25636; //365.24218944;
+
+	// position
+	x.vec[0] = 0.0 / au;
+	x.vec[1] = 0.0 / au;
+	x.vec[2] = 0.0 / au;
+	x.vec[3] = -1.404774254690487E+07 / au;
+	x.vec[4] = -6.825208888976701E+07 / au;
+	x.vec[5] = -4.289153234359771E+06 / au;
+	x.vec[6] = 4.608792459711022E+07 / au;
+	x.vec[7] = -9.859187757812914E+07 / au;
+	x.vec[8] = -4.013349395839728E+06 / au;
+	x.vec[9] = 1.330220204950462E+08 / au;
+	x.vec[10] = 6.704428213911822E+07 / au;
+	x.vec[11] = -4.869615619417280E+03 / au;
+	x.vec[12] = 6.649302608051282E+07 / au;
+	x.vec[13] = 2.186169120955794E+08 / au;
+	x.vec[14] = 2.950343934283048E+06 / au;
+	x.vec[15] = 2.385041315415040E+08 / au;
+	x.vec[16] = 7.182563535822021E+08 / au;
+	x.vec[17] = -8.319703965834230E+06 / au;
+	x.vec[18] = 1.406365269217089E+09 / au;
+	x.vec[19] = -3.235355848531194E+08 / au;
+	x.vec[20] = -5.034784972917840E+07 / au;
+	x.vec[21] = 1.696652737323582E+09 / au;
+	x.vec[22] = 2.385017003636673E+09 / au;
+	x.vec[23] = -1.314104590030003E+07 / au;
+	x.vec[24] = 4.469410600175806E+09 / au;
+	x.vec[25] = -1.293917449213956E+08 / au;
+	x.vec[26] = -1.003300431929548E+08 / au;
+	x.vec[27] = 9.755484379394056E+07 / au;
+	x.vec[28] = -5.078485333345367E+06 / au;
+	x.vec[29] = 3.513104848711042E+07 / au;
+
+	// velocity
+	x.vec[30] = 0.0 / au;
+	x.vec[31] = 0.0 / au;
+	x.vec[32] = 0.0 / au;
+	x.vec[33] = 3.277813191168440E+06 / au / (2.0 * PI) * year;
+	x.vec[34] = -6.360278146631564E+05 / au / (2.0 * PI) * year;
+	x.vec[35] = -3.526256010095223E+05 / au / (2.0 * PI) * year;
+	x.vec[36] = 2.720842478216296E+06 / au / (2.0 * PI) * year;
+	x.vec[37] = 1.271078843814805E+06 / au / (2.0 * PI) * year;
+	x.vec[38] = -1.395383786068834E+05 / au / (2.0 * PI) * year;
+	x.vec[39] = -1.199359489952496E+06 / au / (2.0 * PI) * year;
+	x.vec[40] = 2.288045813113915E+06 / au / (2.0 * PI) * year;
+	x.vec[41] = -1.848630523101974E+02 / au / (2.0 * PI) * year;
+	x.vec[42] = -1.923719964482550E+06 / au / (2.0 * PI) * year;
+	x.vec[43] = 7.871516467196010E+05 / au / (2.0 * PI) * year;
+	x.vec[44] = 6.367679199693510E+04 / au / (2.0 * PI) * year;
+	x.vec[45] = -1.085443281912170E+06 / au / (2.0 * PI) * year;
+	x.vec[46] = 4.092938588278859E+05 / au / (2.0 * PI) * year;
+	x.vec[47] = 2.258591609536928E+04 / au / (2.0 * PI) * year;
+	x.vec[48] = 1.401829990042232E+05 / au / (2.0 * PI) * year;
+	x.vec[49] = 8.122858624024559E+05 / au / (2.0 * PI) * year;
+	x.vec[50] = -1.976963001796068E+04 / au / (2.0 * PI) * year;
+	x.vec[51] = -4.848243855055465E+05 / au / (2.0 * PI) * year;
+	x.vec[52] = 3.142563737252478E+05 / au / (2.0 * PI) * year;
+	x.vec[53] = 7.456670936695803E+03 / au / (2.0 * PI) * year;
+	x.vec[54] = 9.552707840392186E+03 / au / (2.0 * PI) * year;
+	x.vec[55] = 4.728243111876977E+05 / au / (2.0 * PI) * year;
+	x.vec[56] = -9.988672364090331E+03 / au / (2.0 * PI) * year;
+	x.vec[57] = 1.913451980110258E+06 / au / (2.0 * PI) * year;
+	x.vec[58] = -2.739124864901621E+06 / au / (2.0 * PI) * year;
+	x.vec[59] = 2.814734584202529E+06 / au / (2.0 * PI) * year;
+
+	mino2357::RKF78<float_mp> rkf78(mino2357::tol, mino2357::tol);
+
+	for(size_t i=0;;i++){
+		rkf78.Integrate(mino2357::time, mino2357::dt, x);
+		if(i%mino2357::intv == 0){
+			// time
+			std::cerr << mino2357::time / (float_mp(2.0) * float_mp(PI)) << std::endl; // (float_mp(2.0) * float_mp(PI));
+			std::cout << mino2357::time / (float_mp(2.0) * float_mp(PI)) << " "; // (float_mp(2.0) * float_mp(PI));
+			// eccentricity_vector
+			for(size_t i=0; i<mino2357::num/2; i+=3){
+				auto r_x = x.vec[i] - x.vec[0];
+				auto r_y = x.vec[i+1] - x.vec[1];
+				auto r_z = x.vec[i+2] - x.vec[2];
+				auto v_x = x.vec[i+mino2357::num/2];
+				auto v_y = x.vec[i+mino2357::num/2+1];
+				auto v_z = x.vec[i+mino2357::num/2+2];
+				std::cout << " " << r_x << " " << r_y << " " << r_z << " " << v_x << " " << v_y << " " << v_z;
+				auto [e_x, e_y, e_z] = eccentricity_vector<float_mp>(r_x, r_y, r_z, v_x, v_y, v_z);
+				std::cout << " " << e_x << " " << e_y << " " << e_z << " " << sqrt(e_x * e_x + e_y * e_y + e_z * e_z);
+			}
+			auto [l_x, l_y, l_z] = angular_momentum<float_mp>(x);
+			std::cout << " " << l_x << " " << l_y << " " << l_z << " " << sqrt(l_x*l_x+l_y*l_y+l_z*l_z) << std::endl;
+		}
+	}
+}
diff --git a/Tsuchinshan-ATLAS-2024-10-20/makefile b/Tsuchinshan-ATLAS-2024-10-20/makefile
new file mode 100644
index 0000000..f768c02
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/makefile
@@ -0,0 +1,14 @@
+CXX = /usr/bin/g++
+TARGET = main2
+CXXFLAGS = -Wall -Wextra -O2 -march=native -mtune=native -std=c++2b 
+LDLFLAGS = -lstdc++
+SRCS = main2.cpp
+OBJS = $(SRCS:.cpp=.o)
+
+all: $(TARGET)
+
+run: all
+	./$(TARGET)
+
+clean:
+	rm $(TARGET)
diff --git a/Tsuchinshan-ATLAS-2024-10-20/orbit_line.gp b/Tsuchinshan-ATLAS-2024-10-20/orbit_line.gp
new file mode 100644
index 0000000..0a80351
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/orbit_line.gp
@@ -0,0 +1,22 @@
+# set size square
+intv=10
+start=1
+end=1000
+set title "Orbits of C/2023 A3 (Tsuchinshan-ATLAS) from 2024-10-20, mass = 1.0e13"
+set xlabel "x [au]"
+set ylabel "y [au]"
+set zlabel "z [au]"
+set xr [-10:10]
+set yr [-10:10]
+set zr [-10:10]
+splot "log.log"  every intv::start::end u 12:13:14 w l title "Mercury"
+replot "log.log" every intv::start::end u 22:23:24 w l title "Venus"
+replot "log.log" every intv::start::end u 32:33:34 w l title "Earth"
+replot "log.log" every intv::start::end u 42:43:44 w l title "Mars"
+replot "log.log" every intv::start::end u 52:53:54 w l title "Jupiter"
+replot "log.log" every intv::start::end u 62:63:64 w l title "Saturn"
+replot "log.log" every intv::start::end u 72:73:74 w l title "Uranus"
+replot "log.log" every intv::start::end u 82:83:84 w l title "Neptune"
+replot "log.log" every intv::start::end u 92:93:94 w l title "C/2023 A3 (Tsuchinshan-ATLAS)"
+pause 60
+reread
\ No newline at end of file
diff --git a/Tsuchinshan-ATLAS-2024-10-20/orbit_line_2d.gp b/Tsuchinshan-ATLAS-2024-10-20/orbit_line_2d.gp
new file mode 100644
index 0000000..268220c
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/orbit_line_2d.gp
@@ -0,0 +1,23 @@
+# set size square
+intv=10000
+start=1
+end=10009000
+set title "Orbits of C/2023 A3 (Tsuchinshan-ATLAS) from 2024-10-20, mass = 1.0e13"
+set xlabel "x [au]"
+set ylabel "y [au]"
+set size square
+set grid
+set xr [-200:200]
+set yr [-200:200]
+# plot   "log.log" every intv::start::end u  2:3  w l title "Sun"
+plot   "log.log" every intv::start::end u 12:13 w l title "Mercury"
+replot "log.log" every intv::start::end u 22:23 w l title "Venus"
+replot "log.log" every intv::start::end u 32:33 w l title "Earth"
+replot "log.log" every intv::start::end u 42:43 w l title "Mars"
+replot "log.log" every intv::start::end u 52:53 w l title "Jupiter"
+replot "log.log" every intv::start::end u 62:63 w l title "Saturn"
+replot "log.log" every intv::start::end u 72:73 w l title "Uranus"
+replot "log.log" every intv::start::end u 82:83 w l title "Neptune"
+replot "log.log" every intv::start::end u 92:93 w l title "C/2023 A3 (Tsuchinshan-ATLAS)"
+pause 360
+reread
\ No newline at end of file
diff --git a/Tsuchinshan-ATLAS-2024-10-20/parameter.hpp b/Tsuchinshan-ATLAS-2024-10-20/parameter.hpp
new file mode 100644
index 0000000..2a89033
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/parameter.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+using float_mp = double;
+
+namespace mino2357 {
+	constexpr size_t body = 10;
+	constexpr size_t num = 6*body; // body [n]*3[dim]*2[pos, vel]
+	constexpr size_t intv = 100;
+	auto max_dt = float_mp(1.0e-4);
+	auto tol = float_mp(1.0e-10);
+	auto dt = float_mp(1.0e-6);
+	auto time = float_mp(0.0); // zero
+}
+
diff --git a/Tsuchinshan-ATLAS-2024-10-20/test.gp b/Tsuchinshan-ATLAS-2024-10-20/test.gp
new file mode 100644
index 0000000..17a55ab
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/test.gp
@@ -0,0 +1 @@
+plot "log.log" every 1::1 u 1:(sqrt(($12-$22)**2+($13-$23)**2+($14-$24)**2)+sqrt(($22-$32)**2+($23-$33)**2+($24-$34)**2)+sqrt(($32-$42)**2+($33-$43)**2+($34-$44)**2)+sqrt(($42-$52)**2+($43-$53)**2+($44-$54)**2)+sqrt(($52-$62)**2+($53-$63)**2+($54-$64)**2)+sqrt(($62-$72)**2+($63-$73)**2+($64-$74)**2)+sqrt(($72-$82)**2+($73-$83)**2+($74-$84)**2)+sqrt(($82-$92)**2+($83-$93)**2+($84-$94)**2)) w l
diff --git a/Tsuchinshan-ATLAS-2024-10-20/vector.hpp b/Tsuchinshan-ATLAS-2024-10-20/vector.hpp
new file mode 100644
index 0000000..2aa9ae5
--- /dev/null
+++ b/Tsuchinshan-ATLAS-2024-10-20/vector.hpp
@@ -0,0 +1,67 @@
+#pragma once
+
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <iomanip>
+#include <vector>
+
+namespace mino2357{
+    template <typename T=double>
+    class vector{
+    public:
+        std::vector<T> vec;
+        vector(size_t size) noexcept;
+        inline T norm() noexcept;
+    };
+
+    template <typename T>
+    vector<T>::vector(size_t size) noexcept{
+        vec.resize(size);
+    }
+
+    template <typename T>
+    inline T vector<T>::norm() noexcept{
+        T ret = 0.0;
+        for(size_t i=0; i<vec.size(); ++i) {
+            ret += this->vec[i] * this->vec[i];
+        }
+        return sqrt(ret);
+    }
+
+    template <typename T>
+    inline vector<T> operator+(const vector<T>& x, const vector<T>& y) noexcept{
+        vector<T> ret(x.vec.size());
+        for(size_t i=0; i<ret.vec.size(); ++i) {
+            ret.vec[i] = x.vec[i] + y.vec[i];
+        }
+        return ret;
+    }
+
+    template <typename T>
+    inline vector<T> operator-(const vector<T>& x, const vector<T>& y) noexcept{
+        vector<T> ret(x.vec.size());
+        for(size_t i=0; i<ret.vec.size(); ++i) {
+            ret.vec[i] = x.vec[i] - y.vec[i];
+        }
+        return ret;
+    }
+
+    template <typename T>
+    inline T operator*(const vector<T>& x, const vector<T>& y) noexcept{
+        T ret = 0.0;
+        for(size_t i=0; i<x.vec.size(); ++i) {
+            ret += x.vec[i] * y.vec[i];
+        }
+        return ret;
+    }
+
+    template <typename T>
+    inline vector<T> operator*(const T a, const vector<T>& y) noexcept{
+        vector<T> ret(y.vec.size());
+        for(size_t i=0; i<ret.vec.size(); ++i) {
+            ret.vec[i] = a * y.vec[i];
+        }
+        return ret;
+    }
+}