--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Lorenz2Langford.diff	Tue Jul 23 02:36:26 2024 +0900
@@ -0,0 +1,74 @@
+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)