100 Days of Code Log R1D11

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

program logisticmap
    use, intrinsic :: iso_fortran_env, only : dp=>real64, error_unit
    implicit integer(i-k)
    integer, parameter :: l=2001, b = 2001
    real(dp) :: r(l), x0(b), x(l, b)
    integer :: fileno, istat, n
    character(len=1024) :: msg
    print *, "No of iterations: "
    read *, n
    x0 = linspace(0._dp, 1._dp, b)
    r = linspace(3.5_dp, 4._dp, l)

    forall (i=1:l)
        x(i,:) = lm(x0, r(i), n)
    end forall
    print *, "Calculations done"

    open(newunit=fileno, file='logisticmap.dat', iostat=istat, iomsg=msg)
    if (istat /= 0) then
        write(error_unit, '(*(A))') trim(msg)
        stop 1
    end if
    do i=1,l
        write(fileno, *) r(i), x(i, :)
    end do
    close(fileno)
    call execute_command_line('gnuplot -p logisticmap.plt')

    contains
    elemental function lm(x0, r, n) result(xn)
        integer, intent(in) :: n
        real(dp), intent(in) :: x0, r
        real(dp) :: xn
        integer :: i

        xn = x0
        do i = 1, n
            xn = r*xn*(1-xn)
        end do
    end function lm

    function linspace(x0, xn, n) result(x)
        integer, intent(in) :: n
        real(dp), intent(in) :: x0, xn
        real(dp) :: x(n), step
        step = (xn-x0)/(n-1)
        do concurrent(i=1:n)
            x(i) = x0 + step*(i-1)
        end do
    end function linspace
end program logisticmap