Afegir sessió 4 del lab de fenomens
Change-Id: I68a9c725acd6cde3fd43413fe673e5ec7794f7e6
diff --git a/quad10/fenomens/lab/p4/.gitignore b/quad10/fenomens/lab/p4/.gitignore
new file mode 100644
index 0000000..a1a0b85
--- /dev/null
+++ b/quad10/fenomens/lab/p4/.gitignore
@@ -0,0 +1,2 @@
+mc2.dat
+data_out
diff --git a/quad10/fenomens/lab/p4/MC-2.f90 b/quad10/fenomens/lab/p4/MC-2.f90
new file mode 100644
index 0000000..3a2260d
--- /dev/null
+++ b/quad10/fenomens/lab/p4/MC-2.f90
@@ -0,0 +1,116 @@
+program main
+ use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double
+ use iso_fortran_env
+ implicit none
+ integer*4 :: L
+ real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, V1, V2, V3, V4, W(-8:8)
+ integer*4 :: I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0, PBC
+ integer :: genrand_int31
+ character(100) :: NOM
+
+ integer :: SUMI
+ real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM, CV, CHI
+
+ ! This should be replaced by L but it doesn't work dynamically, so we're
+ ! setting 64 as the upper limit which should be higher than L.
+ integer*2 :: S(1:64, 1:64)
+
+ ! Inicialitzem algunes variables sobre el problema
+ namelist /DADES/ L, NOM, TEMP, NSEED, SEED0, MCTOT, MCINI, MCD
+ open(10, file="mc2.dat")
+ read(10, DADES)
+ close(10)
+
+ N = L*L
+
+ ! Obrim el fitxer on escriurem els resultats a cada iteració
+ open(unit = 12, file = "data_out/" // trim(NOM) // ".out")
+
+ ! Cache dels valors de l'exponencial
+ do I = 0, 8
+ W(I) = (2**30 - 1 + 2**30)*exp(-float(I)/TEMP)
+ enddo
+
+ SUMI = 0
+ SUME = 0
+ SUME2 = 0
+
+ SUMM = 0
+ SUMAM = 0
+ SUMM2 = 0
+
+ do SEED = SEED0, SEED0 + NSEED - 1
+ print *, "Seed: ", SEED
+
+ ! Inicialitzem la matriu d'spins aleatòriament
+ call WRITECONFIGSEED(S, L, SEED)
+
+ E = ENERG(S, L)
+ M = MAGNE(S, L)
+ !write (12, *) "0", E, M
+
+ ! Iterem amb el mètode de Montecarlo
+ do IMC = 1, MCTOT
+ do IPAS = 1, N
+ ! LOS LOOPS HACEN COSAS
+ I = MOD(genrand_int31(), L) + 1
+ J = MOD(genrand_int31(), L) + 1
+ V1 = S(PBC(I - 1, L), J)
+ V2 = S(PBC(I + 1, L), J)
+ V3 = S(I, PBC(J - 1, L))
+ V4 = S(I, PBC(J + 1, L))
+ SUMA = V1 + V2 + V3 + V4
+ DIFE = 2*SUMA*S(I, J)
+
+ if (DIFE > 0) then
+ DELTA = genrand_int31()
+ if (DELTA >= W(int(DIFE))) then
+ ! NO ho acceptem
+ cycle
+ endif
+ endif
+
+ ! Sí ho acceptem:
+ S(I, J) = -S(I, J)
+ E = E + DIFE
+ M = M + 2*S(I, J)
+ enddo
+
+ if (mod(IMC, 10000) == 0) then
+ print *, "Iter:", IMC, "Energia:", E, "Magn:", INT(M)
+ endif
+ !write (12, *) IMC, E, M
+
+ if (IMC > MCINI .and. mod(IMC, MCD) == 0) then
+ SUMI = SUMI + 1
+
+ SUME = SUME + E
+ SUME2 = SUME2 + E*E
+
+ SUMM = SUMM + M
+ SUMAM = SUMAM + abs(M)
+ SUMM2 = SUMM2 + M*M
+ endif
+ enddo
+ enddo
+
+ ! Normalitzem els promitjos
+ SUME = SUME/SUMI
+ SUME2 = SUME2/SUMI
+ SUMM = SUMM/SUMI
+ SUMAM = SUMAM/SUMI
+ SUMM2 = SUMM2/SUMI
+ VARE = SUME2 - SUME**2
+ VARM = SUMM2 - SUMM**2
+ CV = (SUME2 - SUME**2)/real(N*TEMP**2)
+ CHI = (SUMM2 - SUMAM**2)/real(N*TEMP)
+
+ ! Guardem a un fiter els promitjos
+ open(unit = 13, file = "data_out/" // trim(NOM) // ".res")
+ !write(13, *) "L, T, <E>, <E**2>, Var(E), <M>, <M**2>, Var(M), C_V, CHI"
+ write(13, *) L, TEMP, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM, CV, CHI
+
+ ! Tanquem els fitxers de sortida
+ close(12)
+ close(13)
+endprogram main
diff --git a/quad10/fenomens/lab/p4/Makefile b/quad10/fenomens/lab/p4/Makefile
new file mode 100644
index 0000000..761b3f7
--- /dev/null
+++ b/quad10/fenomens/lab/p4/Makefile
@@ -0,0 +1,27 @@
+.PHONY: all make_out_folder
+FC=gfortran
+FC_FLAGS=-O3 -Wall #-g -fbacktrace -fcheck=all
+
+all: MC-2.o
+make_out_folder:
+ mkdir -p out data_out
+MC-2.o: make_out_folder mt19937ar.o writeconfigseed.o magne.o energ.o str.o pbc.o MC-2.f90
+ $(FC) $(FC_FLAGS) -c MC-2.f90
+ $(FC) $(FC_FLAGS) MC-2.o writeconfigseed.o mt19937ar.o magne.o energ.o str.o pbc.o -o out/mc2
+mt19937ar.o: mt19937ar.f
+ $(FC) $(FC_FLAGS) -c mt19937ar.f
+writeconfig.o: writeconfig.f90
+ $(FC) $(FC_FLAGS) -c writeconfig.f90
+writeconfigseed.o: writeconfigseed.f90
+ $(FC) $(FC_FLAGS) -c writeconfigseed.f90
+magne.o: magne.f90
+ $(FC) $(FC_FLAGS) -c magne.f90
+energ.o: energ.f90
+ $(FC) $(FC_FLAGS) -c energ.f90
+str.o: str.f90
+ $(FC) $(FC_FLAGS) -c str.f90
+pbc.o: pbc.f90
+ $(FC) $(FC_FLAGS) -c pbc.f90
+
+clean:
+ rm -rf out data_out *.o
diff --git a/quad10/fenomens/lab/p4/README.md b/quad10/fenomens/lab/p4/README.md
new file mode 100644
index 0000000..9406d07
--- /dev/null
+++ b/quad10/fenomens/lab/p4/README.md
@@ -0,0 +1,9 @@
+# Pràctica 4: Promitjos en funció de T. Dependència amb L. Extrapolació a L tendint a inf.
+**Data: 12 de maig de 2022**
+
+## Instruccions per compilar
+Per compilar, executeu `make all`. Els executables es trobaran a la carpeta `out`.
+
+## Programes
+- `dependenciaEnT.bash`: genera gràfiques de la capacitat calorífica i la susceptibilitat magnètica en funció de la temperatura.
+- `dependenciaEnL.bash`: el mateix que el programa anterior, però a les gràfiques apareixen les corbes per Ls diferents.
diff --git a/quad10/fenomens/lab/p4/dependenciaEnL.bash b/quad10/fenomens/lab/p4/dependenciaEnL.bash
new file mode 100755
index 0000000..05badba
--- /dev/null
+++ b/quad10/fenomens/lab/p4/dependenciaEnL.bash
@@ -0,0 +1,17 @@
+#!/bin/bash
+l_list="8 16 32 64"
+out_folder="data_out/depEnL"
+
+rm -rf $out_folder
+mkdir -p $out_folder
+
+for l in $l_list; do
+ echo "============================="
+ echo "Executant per L = $l"
+ echo "============================="
+
+ ./dependenciaEnT.bash $l
+ cp data_out/dep_en_T.dat "$out_folder/dep_en_T_L$l.dat"
+done
+
+./graphs/dependenciaEnL.gnu
diff --git a/quad10/fenomens/lab/p4/dependenciaEnLPost.sh b/quad10/fenomens/lab/p4/dependenciaEnLPost.sh
new file mode 100755
index 0000000..d9230fd
--- /dev/null
+++ b/quad10/fenomens/lab/p4/dependenciaEnLPost.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+./graphs/dependenciaEnL.gnu
diff --git a/quad10/fenomens/lab/p4/dependenciaEnT.bash b/quad10/fenomens/lab/p4/dependenciaEnT.bash
new file mode 100755
index 0000000..e9f5a6a
--- /dev/null
+++ b/quad10/fenomens/lab/p4/dependenciaEnT.bash
@@ -0,0 +1,34 @@
+#!/bin/bash
+temps="1.8 2.0 2.15 2.25 2.3 2.35 2.45 2.6 3.0"
+out="data_out/dep_en_T.dat"
+out_partial="partial_dep_en_T"
+l=${1:-32}
+
+rm -f $out
+touch $out
+
+for temp in $temps; do
+ echo "-----------------------------"
+ echo "Executant per T = $temp"
+ echo "-----------------------------"
+
+ cat <<EOF > mc2.dat
+&DADES
+L=$l
+NOM="$out_partial",
+TEMP=$temp,
+NSEED=10,
+SEED0=117654,
+MCTOT=20000,
+MCINI=2000,
+MCD=20
+&END
+EOF
+ ./out/mc2 $temp
+ cat data_out/$out_partial.res >> $out
+done
+
+rm -f data_out/$out_partial.res
+rm -f data_out/$out_partial.out
+
+./graphs/dependenciaEnT.gnu
diff --git a/quad10/fenomens/lab/p4/dependenciaEnTPost.sh b/quad10/fenomens/lab/p4/dependenciaEnTPost.sh
new file mode 100755
index 0000000..067eecd
--- /dev/null
+++ b/quad10/fenomens/lab/p4/dependenciaEnTPost.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+./graphs/dependenciaEnT.gnu
diff --git a/quad10/fenomens/lab/p4/energ.f90 b/quad10/fenomens/lab/p4/energ.f90
new file mode 100644
index 0000000..ba7841c
--- /dev/null
+++ b/quad10/fenomens/lab/p4/energ.f90
@@ -0,0 +1,13 @@
+real*8 function ENERG(S, L)
+ integer*2 :: S(1:64, 1:64)
+ integer*4 :: I, J, L, PBC
+ real*8 :: ENE
+ ENE = 0.0d0
+ do I = 1, L
+ do J = 1, L
+ ENE = ENE - S(I, J)*S(PBC(I + 1, L), J) - S(I, J)*S(I, PBC(J + 1, L))
+ enddo
+ enddo
+ ENERG = ENE
+ return
+endfunction
diff --git a/quad10/fenomens/lab/p4/generateTemps.js b/quad10/fenomens/lab/p4/generateTemps.js
new file mode 100644
index 0000000..6352a62
--- /dev/null
+++ b/quad10/fenomens/lab/p4/generateTemps.js
@@ -0,0 +1,8 @@
+let b = 1.4
+let s = ""
+while (b <= 3.4) {
+ s += b.toPrecision(2) + " ";
+ b += 0.1
+}
+s += b.toPrecision(2) + " ";
+console.log(s)
diff --git a/quad10/fenomens/lab/p4/generateTemps2.js b/quad10/fenomens/lab/p4/generateTemps2.js
new file mode 100644
index 0000000..035b98c
--- /dev/null
+++ b/quad10/fenomens/lab/p4/generateTemps2.js
@@ -0,0 +1,8 @@
+let b = 2.1
+let s = ""
+while (b <= 2.4) {
+ s += b.toPrecision(3) + " ";
+ b += 0.02
+}
+s += b.toPrecision(3) + " ";
+console.log(s)
diff --git a/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu b/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu
new file mode 100755
index 0000000..b7d239b
--- /dev/null
+++ b/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu
@@ -0,0 +1,15 @@
+#!/usr/bin/env -S gnuplot -c
+outputfile = 'data_out/dep_en_L' # Nom de la imatge resultant (sense extensió)
+datafilepre = 'data_out/depEnL/dep_en_T_L'
+datafilepost = '.dat'
+LS="8 16 32 64"
+
+set terminal svg dashed size 600, 1200 font "Computer Modern,Tinos,Helvetica,15"
+set output outputfile.'.svg'
+
+set multiplot layout 2,1
+
+set title "Capacitat calorífica"
+plot for [L in LS] datafilepre . L . datafilepost using 2:10 with points title "L = ".L
+set title "Susceptibilitat magnètica"
+plot for [L in LS] datafilepre . L . datafilepost using 2:11 with points title "L = ".L
diff --git a/quad10/fenomens/lab/p4/graphs/dependenciaEnT.gnu b/quad10/fenomens/lab/p4/graphs/dependenciaEnT.gnu
new file mode 100755
index 0000000..18a435a
--- /dev/null
+++ b/quad10/fenomens/lab/p4/graphs/dependenciaEnT.gnu
@@ -0,0 +1,15 @@
+#!/usr/bin/env -S gnuplot -c
+outputfile = 'data_out/dep_en_T' # Nom de la imatge resultant (sense extensió)
+datafile = 'data_out/dep_en_T.dat'
+L=32
+TEMPS="1500 1800 2500 3500 4500"
+
+set terminal svg dashed size 600, 1200 font "Computer Modern,Tinos,Helvetica,15"
+set output outputfile.'.svg'
+
+set multiplot layout 2,1
+
+set title "Capacitat calorífica"
+plot datafile using 2:10 with points
+set title "Susceptibilitat magnètica"
+plot datafile using 2:11 with points
diff --git a/quad10/fenomens/lab/p4/magne.f90 b/quad10/fenomens/lab/p4/magne.f90
new file mode 100644
index 0000000..be56bc9
--- /dev/null
+++ b/quad10/fenomens/lab/p4/magne.f90
@@ -0,0 +1,17 @@
+real*8 function MAGNE(S, L)
+ implicit none
+ integer*4, intent(in) :: L
+ integer*2, intent(out) :: S(1:64,1:64)
+ integer*4 :: i, j
+ real*8 :: mag
+ mag = 0d0
+
+ do i = 1, L
+ do j = 1, L
+ mag = mag + S(i, j)
+ enddo
+ enddo
+
+ magne = mag
+ return
+end function MAGNE
diff --git a/quad10/fenomens/lab/p4/mt19937ar.f b/quad10/fenomens/lab/p4/mt19937ar.f
new file mode 100644
index 0000000..be81d16
--- /dev/null
+++ b/quad10/fenomens/lab/p4/mt19937ar.f
@@ -0,0 +1,282 @@
+c
+c A C-program for MT19937, with initialization improved 2002/1/26.
+c Coded by Takuji Nishimura and Makoto Matsumoto.
+c
+c Before using, initialize the state by using init_genrand(seed)
+c or init_by_array(init_key, key_length).
+c
+c Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+c All rights reserved.
+c Copyright (C) 2005, Mutsuo Saito,
+c All rights reserved.
+c
+c Redistribution and use in source and binary forms, with or without
+c modification, are permitted provided that the following conditions
+c are met:
+c
+c 1. Redistributions of source code must retain the above copyright
+c notice, this list of conditions and the following disclaimer.
+c
+c 2. Redistributions in binary form must reproduce the above copyright
+c notice, this list of conditions and the following disclaimer in the
+c documentation and/or other materials provided with the distribution.
+c
+c 3. The names of its contributors may not be used to endorse or promote
+c products derived from this software without specific prior written
+c permission.
+c
+c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+c "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+c LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+c A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+c CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+c EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+c PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+c PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+c LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+c NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+c SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+c
+c
+c Any feedback is very welcome.
+c http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+c email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+c
+c-----------------------------------------------------------------------
+c FORTRAN77 translation by Tsuyoshi TADA. (2005/12/19)
+c
+c ---------- initialize routines ----------
+c subroutine init_genrand(seed): initialize with a seed
+c subroutine init_by_array(init_key,key_length): initialize by an array
+c
+c ---------- generate functions ----------
+c integer function genrand_int32(): signed 32-bit integer
+c integer function genrand_int31(): unsigned 31-bit integer
+c double precision function genrand_real1(): [0,1] with 32-bit resolution
+c double precision function genrand_real2(): [0,1) with 32-bit resolution
+c double precision function genrand_real3(): (0,1) with 32-bit resolution
+c double precision function genrand_res53(): (0,1) with 53-bit resolution
+c
+c This program uses the following non-standard intrinsics.
+c ishft(i,n): If n>0, shifts bits in i by n positions to left.
+c If n<0, shifts bits in i by n positions to right.
+c iand (i,j): Performs logical AND on corresponding bits of i and j.
+c ior (i,j): Performs inclusive OR on corresponding bits of i and j.
+c ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
+c
+c-----------------------------------------------------------------------
+c initialize mt(0:N-1) with a seed
+c-----------------------------------------------------------------------
+ subroutine init_genrand(s)
+ integer s
+ integer N
+ integer DONE
+ integer ALLBIT_MASK
+ parameter (N=624)
+ parameter (DONE=123456789)
+ integer mti,initialized
+ integer mt(0:N-1)
+ common /mt_state1/ mti,initialized
+ common /mt_state2/ mt
+ common /mt_mask1/ ALLBIT_MASK
+c
+ call mt_initln
+ mt(0)=iand(s,ALLBIT_MASK)
+ do 100 mti=1,N-1
+ mt(mti)=1812433253*
+ & ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti
+ mt(mti)=iand(mt(mti),ALLBIT_MASK)
+ 100 continue
+ initialized=DONE
+c
+ return
+ end
+c-----------------------------------------------------------------------
+c initialize by an array with array-length
+c init_key is the array for initializing keys
+c key_length is its length
+c-----------------------------------------------------------------------
+ subroutine init_by_array(init_key,key_length)
+ integer init_key(0:*)
+ integer key_length
+ integer N
+ integer ALLBIT_MASK
+ integer TOPBIT_MASK
+ parameter (N=624)
+ integer i,j,k
+ integer mt(0:N-1)
+ common /mt_state2/ mt
+ common /mt_mask1/ ALLBIT_MASK
+ common /mt_mask2/ TOPBIT_MASK
+c
+ call init_genrand(19650218)
+ i=1
+ j=0
+ do 100 k=max(N,key_length),1,-1
+ mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1664525)
+ & +init_key(j)+j
+ mt(i)=iand(mt(i),ALLBIT_MASK)
+ i=i+1
+ j=j+1
+ if(i.ge.N)then
+ mt(0)=mt(N-1)
+ i=1
+ endif
+ if(j.ge.key_length)then
+ j=0
+ endif
+ 100 continue
+ do 200 k=N-1,1,-1
+ mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1566083941)-i
+ mt(i)=iand(mt(i),ALLBIT_MASK)
+ i=i+1
+ if(i.ge.N)then
+ mt(0)=mt(N-1)
+ i=1
+ endif
+ 200 continue
+ mt(0)=TOPBIT_MASK
+c
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on [0,0xffffffff]-interval
+c-----------------------------------------------------------------------
+ function genrand_int32()
+ integer genrand_int32
+ integer N,M
+ integer DONE
+ integer UPPER_MASK,LOWER_MASK,MATRIX_A
+ integer T1_MASK,T2_MASK
+ parameter (N=624)
+ parameter (M=397)
+ parameter (DONE=123456789)
+ integer mti,initialized
+ integer mt(0:N-1)
+ integer y,kk
+ integer mag01(0:1)
+ common /mt_state1/ mti,initialized
+ common /mt_state2/ mt
+ common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
+ common /mt_mag01/ mag01
+c
+ if(initialized.ne.DONE)then
+ call init_genrand(21641)
+ endif
+c
+ if(mti.ge.N)then
+ do 100 kk=0,N-M-1
+ y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK))
+ mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
+ 100 continue
+ do 200 kk=N-M,N-1-1
+ y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK))
+ mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
+ 200 continue
+ y=ior(iand(mt(N-1),UPPER_MASK),iand(mt(0),LOWER_MASK))
+ mt(kk)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
+ mti=0
+ endif
+c
+ y=mt(mti)
+ mti=mti+1
+c
+ y=ieor(y,ishft(y,-11))
+ y=ieor(y,iand(ishft(y,7),T1_MASK))
+ y=ieor(y,iand(ishft(y,15),T2_MASK))
+ y=ieor(y,ishft(y,-18))
+c
+ genrand_int32=y
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on [0,0x7fffffff]-interval
+c-----------------------------------------------------------------------
+ function genrand_int31()
+ integer genrand_int31
+ integer genrand_int32
+ genrand_int31=int(ishft(genrand_int32(),-1))
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on [0,1]-real-interval
+c-----------------------------------------------------------------------
+ function genrand_real1()
+ double precision genrand_real1,r
+ integer genrand_int32
+ r=dble(genrand_int32())
+ if(r.lt.0.d0)r=r+2.d0**32
+ genrand_real1=r/4294967295.d0
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on [0,1)-real-interval
+c-----------------------------------------------------------------------
+ function genrand_real2()
+ double precision genrand_real2,r
+ integer genrand_int32
+ r=dble(genrand_int32())
+ if(r.lt.0.d0)r=r+2.d0**32
+ genrand_real2=r/4294967296.d0
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on (0,1)-real-interval
+c-----------------------------------------------------------------------
+ function genrand_real3()
+ double precision genrand_real3,r
+ integer genrand_int32
+ r=dble(genrand_int32())
+ if(r.lt.0.d0)r=r+2.d0**32
+ genrand_real3=(r+0.5d0)/4294967296.d0
+ return
+ end
+c-----------------------------------------------------------------------
+c generates a random number on [0,1) with 53-bit resolution
+c-----------------------------------------------------------------------
+ function genrand_res53()
+ double precision genrand_res53
+ integer genrand_int32
+ double precision a,b
+ a=dble(ishft(genrand_int32(),-5))
+ b=dble(ishft(genrand_int32(),-6))
+ if(a.lt.0.d0)a=a+2.d0**32
+ if(b.lt.0.d0)b=b+2.d0**32
+ genrand_res53=(a*67108864.d0+b)/9007199254740992.d0
+ return
+ end
+c-----------------------------------------------------------------------
+c initialize large number (over 32-bit constant number)
+c-----------------------------------------------------------------------
+ subroutine mt_initln
+ integer ALLBIT_MASK
+ integer TOPBIT_MASK
+ integer UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
+ integer mag01(0:1)
+ common /mt_mask1/ ALLBIT_MASK
+ common /mt_mask2/ TOPBIT_MASK
+ common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
+ common /mt_mag01/ mag01
+CC TOPBIT_MASK = Z'80000000'
+CC ALLBIT_MASK = Z'ffffffff'
+CC UPPER_MASK = Z'80000000'
+CC LOWER_MASK = Z'7fffffff'
+CC MATRIX_A = Z'9908b0df'
+CC T1_MASK = Z'9d2c5680'
+CC T2_MASK = Z'efc60000'
+ TOPBIT_MASK=1073741824
+ TOPBIT_MASK=ishft(TOPBIT_MASK,1)
+ ALLBIT_MASK=2147483647
+ ALLBIT_MASK=ior(ALLBIT_MASK,TOPBIT_MASK)
+ UPPER_MASK=TOPBIT_MASK
+ LOWER_MASK=2147483647
+ MATRIX_A=419999967
+ MATRIX_A=ior(MATRIX_A,TOPBIT_MASK)
+ T1_MASK=489444992
+ T1_MASK=ior(T1_MASK,TOPBIT_MASK)
+ T2_MASK=1875247104
+ T2_MASK=ior(T2_MASK,TOPBIT_MASK)
+ mag01(0)=0
+ mag01(1)=MATRIX_A
+ return
+ end
diff --git a/quad10/fenomens/lab/p4/pbc.f90 b/quad10/fenomens/lab/p4/pbc.f90
new file mode 100644
index 0000000..09cfff3
--- /dev/null
+++ b/quad10/fenomens/lab/p4/pbc.f90
@@ -0,0 +1,11 @@
+integer*4 function PBC(N, L)
+ integer*4, intent(in) :: N, L
+ if (N == 0) then
+ PBC = L
+ else if (N == L + 1) then
+ PBC = 1
+ else
+ PBC = N
+ endif
+ return
+endfunction
diff --git a/quad10/fenomens/lab/p4/str.f90 b/quad10/fenomens/lab/p4/str.f90
new file mode 100644
index 0000000..31cd529
--- /dev/null
+++ b/quad10/fenomens/lab/p4/str.f90
@@ -0,0 +1,6 @@
+character(len=20) function str(k)
+ ! "Convert an integer to string."
+ integer, intent(in) :: k
+ write (str, *) k
+ str = adjustl(str)
+end function str
diff --git a/quad10/fenomens/lab/p4/writeconfig.f90 b/quad10/fenomens/lab/p4/writeconfig.f90
new file mode 100644
index 0000000..b356d17
--- /dev/null
+++ b/quad10/fenomens/lab/p4/writeconfig.f90
@@ -0,0 +1,20 @@
+subroutine WRITECONFIG(S, L)
+ implicit none
+ integer*4, intent(in) :: L
+ integer*2, intent(out) :: S(1:64,1:64)
+ integer*4 :: SEED, i, j
+ real*8 :: genrand_real2
+
+ SEED = 23456
+ call init_genrand(SEED)
+
+ do i = 1, L
+ do j = 1, L
+ if (genrand_real2() < .5d0) then
+ S(i, j) = 1
+ else
+ S(i, j) = -1
+ endif
+ enddo
+ enddo
+end subroutine WRITECONFIG
diff --git a/quad10/fenomens/lab/p4/writeconfigseed.f90 b/quad10/fenomens/lab/p4/writeconfigseed.f90
new file mode 100644
index 0000000..3d4bd13
--- /dev/null
+++ b/quad10/fenomens/lab/p4/writeconfigseed.f90
@@ -0,0 +1,21 @@
+subroutine WRITECONFIGSEED(S, L, SEED)
+ use iso_fortran_env
+ implicit none
+ integer*4, intent(in) :: L
+ integer*4, intent(in) :: SEED
+ integer*2, intent(out) :: S(1:64,1:64)
+ integer*4 :: i, j
+ integer(int32) :: genrand_int31
+
+ call init_genrand(SEED)
+
+ do i = 1, L
+ do j = 1, L
+ if (genrand_int31() < 2**30) then
+ S(i, j) = 1
+ else
+ S(i, j) = -1
+ endif
+ enddo
+ enddo
+end subroutine WRITECONFIGSEED