! Author: Andreas S. Krogen
! Computes the partition function z, internal energy U, heat capacity Ch and
! magnetisation per spin m as functions of temperature T in the 1D Ising model.
! The boundary conditions are periodic.
!
! Written as part of a numerical exercise in the course TFY4230 at NTNU (2016).
! Compile with: gfortran -O3 -o ising1d ising_1d_ring.f95

program ising_1d_ring
    implicit none

    ! ================ Declaration of variables ================
    ! Constants
    real, parameter :: kb = 1.38064852E-23 ! Boltzmann constant

    ! Model parameters
    real :: h, J_spin
    integer :: N
    integer :: Ntotal ! Ntotal = N*N = total number of lattice points

    ! Other variables
    integer, dimension(:), allocatable :: L ! Array of spin configurations
    integer :: alloc_stat
    ! npts: number of points at which z, U, etc. are evaluated.
    integer, parameter :: npts = 100
    real, parameter :: Tmin = 0.01, Tmax = 8.0
    real(kind=8), dimension(1:npts) :: z_arr, U_arr, H_squared_arr, Ch_arr
    real(kind=8), dimension(1:npts) :: ln_z_arr, m_arr, T_arr
    real(kind=8) :: T, Tstep, z, beta, E, bfactor
    integer :: i, j, ti
    integer(kind=8) :: bitfield, bits

    ! Filenames
    character(len=32) :: ln_z_fname, U_fname, m_fname, Ch_fname

    ! ================ Setting model parameters and stuff ================
    h = 0.0
    ! Iterate over N
    do N = 2, 5
    Ntotal = N*N
    J_spin = 1.0

    ! Allocate memory for lattice
    allocate (L(0:Ntotal-1), stat=alloc_stat)
    if (alloc_stat /= 0) then
        print *, 'Error: memory allocation.'
        stop
    endif
    ! Generate evenly spaced array of T values
    ! and set the elements of the other arrays to zero
    Tstep = (Tmax - Tmin) / (npts-1)
    T = Tmin
    do ti=1, npts
       T_arr(ti) = T
       T = T + Tstep
       z_arr(ti) = 0
       ln_z_arr(ti) = 0
       U_arr(ti) = 0
       H_squared_arr(ti) = 0
       Ch_arr(ti) = 0
       m_arr(ti) = 0
    enddo 
   
    ! ================ And now, let's crunch some numbers... ================
    !
    ! Iterate over each microstate
    ! Generate all 2**Ntotal permutations of spin, bit hacking style
    do bitfield=0, 2**Ntotal - 1
        bits = bitfield
        do i=0, Ntotal - 1
            ! map lsb: 0 --> -1, 1 --> 1, then do a right shift
            L(i) = merge(-1, 1, iand(bits, 1) == 0)
            bits = ishft(bits, -1) ! logical shift right
        enddo 

        ! Iterate over temperatures T
        do ti=1, npts
            T = T_arr(ti)
            beta = 1.0 / T ! Omit k to avoid enormously huge numbers
            E = calc_H()
            bfactor = exp(-beta*E)
            z_arr(ti) = z_arr(ti) + bfactor
            U_arr(ti) = U_arr(ti) + E * bfactor
            H_squared_arr(ti) = H_squared_arr(ti) + E**2 * bfactor 
            m_arr(ti) = m_arr(ti) + abs(sum(L)) * bfactor
        enddo
    enddo 

    ! Normalisation
    ln_z_arr = log(z_arr)
    U_arr = U_arr / z_arr
    H_squared_arr = H_squared_arr / z_arr
    m_arr = m_arr / (Ntotal*z_arr)
    Ch_arr = (H_squared_arr - U_arr**2) / T_arr**2 / Ntotal
    U_arr = U_arr / Ntotal

    ! ================ Save results ================
    ! Filenames
    write (ln_z_fname, "(A5,I0.2,A1,F3.1,A7)") "lnz_N", Ntotal, "h", h, "_1d.dat"
    write (U_fname, "(A3,I0.2,A1,F3.1,A7)") "U_N", Ntotal, "h", h, "_1d.dat"
    write (Ch_fname, "(A4,I0.2,A1,F3.1,A7)") "Ch_N", Ntotal, "h", h, "_1d.dat"
    write (m_fname, "(A3,I0.2,A1,F3.1,A7)") "m_N", Ntotal, "h", h, "_1d.dat"

    ! Save data to files
    call save_data(ln_z_fname, T_arr, ln_z_arr)
    call save_data(U_fname, T_arr, U_arr)
    call save_data(Ch_fname, T_arr, Ch_arr)
    call save_data(m_fname, T_arr, m_arr)

    deallocate(L)

    enddo ! for each N

! ================ Subroutines and functions ================
contains
    ! Returns the energy of the current microstate
    function calc_H() result(E)
        implicit none
        real(8) :: E
        integer :: i, j, Sij

        ! Contribution from neighbour spins
        Sij = 0
        ! Iterating over each spin
        do i=0, Ntotal-1
            ! Periodic BCs
            Sij = Sij + L(i) * L(mod(i+1, Ntotal)) 
        enddo

        E = -1 *(J_spin * Sij + sum(L) * h)
    end function calc_H

    ! Saves points (x, y) to file.
    subroutine save_data(filename, xs, ys)
        implicit none
        character(len=32), intent(in) :: filename
        real(8), dimension(1:), intent(in) :: xs, ys
        integer :: size_x, size_y, i
        
        size_x = size(xs)
        size_y = size(ys)
        if (size_x /= size_y) then
            print *, "Error: the size of xs does not match the size of ys."
        else
            open(unit=1, file=filename, action='write')
            do i=1, size_x
                write(1, *) xs(i), ' ', ys(i)
            enddo
            close(unit=1)
        endif
    end subroutine save_data

end program ising_1d_ring
