Mercurial > op > rk_sakamoti

changeset 5:0ca8347e4c47

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add rk_Lorenz.R and rksub.f90
author "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
date Tue, 23 Jul 2024 03:31:21 +0900
parents 8026ceb4de7a
children d90925ba642e
files Makefile rk_Lorenz.R rkmain.f90 rksub.f90
diffstat 4 files changed, 45 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile	Tue Jul 23 01:48:31 2024 +0900
+++ b/Makefile	Tue Jul 23 03:31:21 2024 +0900
@@ -13,7 +13,7 @@
 	gfortran -c rkmain.f90
 
 rkmod.o: rkmod.f90 rungekutta.o
-	gfortran -c rkmod.f90
+	gfortran -c rkmod.f90 -fPIC
 
 rungekutta.o: rungekutta.f90
 	gfortran -c rungekutta.f90
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rk_Lorenz.R	Tue Jul 23 03:31:21 2024 +0900
@@ -0,0 +1,29 @@
+dll <- paste("rksub", .Platform$dynlib.ext, sep="")
+dyn.load(dll)
+
+max_i <- as.integer(5000)
+
+r <- .Fortran("rk_sakamoti",
+	as.double(c(0.5, 0.5, 0.5)), 
+	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", 
+    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", 
+    xlab = colnames(m)[2],
+    ylab = colnames(m)[3],
+    zlab = colnames(m)[4],
+    cex.main = 3,
+    cex.lab = 2)
+dev.off()
+
--- a/rkmain.f90	Tue Jul 23 01:48:31 2024 +0900
+++ b/rkmain.f90	Tue Jul 23 03:31:21 2024 +0900
@@ -3,10 +3,13 @@
   implicit none
   real(kind=dpkind),dimension(3) :: x = [0.5, 0.5, 0.5]
   integer :: max_i = 5000
-  real(kind=dpkind), allocatable, dimension(:, :) :: r
+!  real(kind=dpkind), allocatable, dimension(:, :) :: r
+  real(kind=dpkind), dimension(5000, 4) :: r
 
-  allocate(r(max_i, 4))
+!  allocate(r(max_i, 4))
   call calc(x, max_i, r)
-  print *, r
+  open(17,file='rkmain.dat', form='unformatted', access = "stream", status='replace')
+  write(17) r
+  close(17)
 end program
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rksub.f90	Tue Jul 23 03:31:21 2024 +0900
@@ -0,0 +1,9 @@
+subroutine rk_sakamoti(x, max_i, r)
+  use rkmod
+  implicit none
+  real(kind=dpkind),dimension(3),intent(in) :: x
+  integer, intent(in) :: max_i
+  real(kind=dpkind), dimension(max_i, 4),intent(out) :: r
+  call calc(x, max_i, r)
+end subroutine 
+