Afegir sessió 3 del lab de fenòmens
Change-Id: Ic21d679decc844e3c2899d6da4cca48301dbae0c
diff --git a/quad10/fenomens/lab/p3/MC-2.f90 b/quad10/fenomens/lab/p3/MC-2.f90
new file mode 100644
index 0000000..20efdca
--- /dev/null
+++ b/quad10/fenomens/lab/p3/MC-2.f90
@@ -0,0 +1,117 @@
+program main
+ use, intrinsic :: iso_c_binding, only: sp=>c_float, dp=>c_double
+ implicit none
+ integer*4, PARAMETER :: L = 32
+ integer*2 :: S(1:L, 1:L)
+ real*8 :: MAGNE, ENERG, TEMP, E, DIFE, DELTA, M, SUMA, W(-8:8), genrand_real2
+ integer*4 :: PBC(0:L+1), I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0
+ character :: TEMPSTRING*100, str*20
+ character(:), allocatable :: NOM
+
+ integer :: SUMI
+ real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM
+
+ if (command_argument_count() < 1) then
+ print *, "Error: please supply a command-line argument with the"
+ print *, "temperature, via './mc1 {TEMPERATURE}', where {TEMPERATURE}"
+ print *, "is the temperature."
+ stop
+ endif
+
+ ! Inicialitzem algunes variables sobre el problema
+ call get_command_argument(1, TEMPSTRING) ! Ex: 2.0d0
+ TEMPSTRING = trim(adjustl(TEMPSTRING))
+ read (TEMPSTRING, *) TEMP
+ NSEED = 10
+ SEED0 = 117654
+ MCTOT = 10000
+ MCINI = 1000
+ MCD = 10
+ N = L*L
+ NOM = "out_data/" // "SIM-L-" // trim(str(L)) // "-TEMP-" // &
+ trim(str(int(1000*TEMP))) // "-MCTOT-" // trim(str(MCTOT))
+
+ ! Obrim el fitxer on escriurem els resultats a cada iteració
+ open(unit = 12, file = NOM // ".out")
+
+ ! Inicialitzem variable que ens fa latent la periodicitat
+ PBC(0) = L
+ PBC(L + 1) = 1
+ do I = 1, L
+ PBC(I) = I
+ enddo
+
+ ! Cache dels valors de l'exponencial
+ do I = -8, 8
+ W(I) = exp(-float(I)/TEMP)
+ enddo
+
+ 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, PBC)
+ M = MAGNE(S, L)
+ print *, "Energia inicial:", E, "Magnet. inicial:", M
+ 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 = INT(genrand_real2()*L) + 1
+ J = INT(genrand_real2()*L) + 1
+ SUMA = S(PBC(I - 1), J) + S(PBC(I + 1), J) + S(I, PBC(J - 1)) + S(I, PBC(J + 1))
+ DIFE = 2*SUMA*S(I, J)
+
+ if (DIFE > 0) then
+ DELTA = genrand_real2()
+ 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, 1000) == 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*SUME
+ VARM = SUMM2 - SUMM*SUMM
+
+ ! Guardem a un fiter els promitjos
+ open(unit = 13, file = NOM // ".res")
+ write(13, *) L, TEMP, SUMI, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM
+
+ ! Tanquem els fitxers de sortida
+ close(12)
+ close(13)
+endprogram main