diff -r 925250ea21a4 rk_Lorenz.R
--- a/rk_Lorenz.R	Sun Jul 21 06:20:46 2024 +0900
+++ b/rk_Lorenz.R	Tue Jul 23 02:35:22 2024 +0900
@@ -1,25 +1,25 @@
 dll <- paste("rksub", .Platform$dynlib.ext, sep="")
 dyn.load(dll)
 
-max_i <- as.integer(5000)
+max_i <- as.integer(5000*10)
 
 r <- .Fortran("rk_sakamoti",
-	as.double(c(0.5, 0.5, 0.5)), 
-	max_i,
+	as.double(c(1.0, 1.0, 1.0)), 
+		max_i,
         matrix(0, max_i, 4))
 dyn.unload(dll)
 
 m <- r[[3]]
 colnames(m) <- c("t", "x", "y", "z")
 
-png(filename = "Lorenz system.png", 
+png(filename = "Langford.png", 
     width = 800, height = 600, bg="white", type="cairo")
 library(scatterplot3d)
 scatterplot3d(m[,2], m[,3], m[,4], 
     highlight.3d = FALSE, 
     col.axis = "blue", col.grid = "lightblue", 
     type = "l", color = "purple", lwd = 2,
-    main = "Lorenz system", 
+    main = "Langford Attractor", 
     xlab = colnames(m)[2],
     ylab = colnames(m)[3],
     zlab = colnames(m)[4],
@@ -27,3 +27,16 @@
     cex.lab = 2)
 dev.off()
 
+subplot <- function(m, n){
+	m2 <- m[n, ]
+	scatterplot3d(m2[,2], m2[,3], m2[,4], 
+	    highlight.3d = FALSE, 
+	    col.axis = "blue", col.grid = "lightblue", 
+	    type = "l", color = "purple", lwd = 2,
+	    main = "Langford Attractor", 
+	    xlab = colnames(m)[2],
+	    ylab = colnames(m)[3],
+	    zlab = colnames(m)[4],
+	    cex.main = 3,
+	    cex.lab = 2)
+}
diff -r 925250ea21a4 rkmod.f90
--- a/rkmod.f90	Sun Jul 21 06:20:46 2024 +0900
+++ b/rkmod.f90	Tue Jul 23 02:35:22 2024 +0900
@@ -10,10 +10,15 @@
       real(kind=dpkind) :: time
       real(kind=dpkind),dimension(:) :: wk
       real(kind=dpkind),dimension(size(wk)) :: rhside
-      real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
-      rhside(1)=a*(wk(2)-wk(1))
-      rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
-      rhside(3)=wk(1)*wk(2)-c*wk(3)
+!      real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
+!      rhside(1)=a*(wk(2)-wk(1))
+!      rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
+!      rhside(3)=wk(1)*wk(2)-c*wk(3)
+      real(kind=dpkind),parameter :: a = 0.6
+      rhside(1) = (wk(3) - 0.7)*wk(1) - 3.5*wk(2)
+      rhside(2) = 3.5*wk(1) + (wk(3) - 0.7)*wk(2)
+      rhside(3) = a + wk(3) - wk(3)**3 / 3 - (wk(1)**2 + wk(2)**2)*(1 + 0.25*wk(3))
+
     end function
     !Rから呼び出すサブルーチン
     subroutine calc(x, max_i, r)