!***********************************************
!***********************************************
! 易辛二维度模拟正负电荷的2D阵列进化
! 每个很会“翻转”与邻居一致
!***********************************************
! SINOMACH
!
! 易辛二维度模拟正负电荷的2D阵列进化
! 每个很会“翻转”与邻居一致
!
! 周勇，PotsyYZ，GoldbergJr.,JohnBktJr.
!
! Sept.2022
!***********************************************
Program main

!*********************************************************************72
!
!c MAIN is the main program for ISING_2D_SIMULATION.
!
! Usage:
!
! ising_2d_simulation m n iterations thresh seed
!
! * M, N, the number of rows and columns.
! * ITERATIONS, the number of iterations.
! * THRESH, the threshhold.
! * SEED, a seed for the random number generator.
!
Implicit None

Integer i
Integer iterations
Integer m
Integer n
Double Precision prob(5)
Integer seed
Integer step
Double Precision thresh

Save prob

Data prob/0.97D+00, 0.86D+00, 0.51D+00, 0.16D+00, 0.03D+00/

Call timestamp()
Write (, '(a)') ' '
Write (, '(a)') 'ISING_2D_SIMULATION'
Write (, '(a)') ' FORTRAN95,2013&18 version'
Write (, '(a)') ' Monte Carlo simulation of a 2D Ising model.'
!
! Get input from commandline or user.
!
Call get_input(m, n, iterations, thresh, seed)

Write (, '(a)') ' '
Write (, '(a,i8)') ' The number of rows is M = ', m
Write (, '(a,i8)') ' The number of columns is N = ', n
Write (, '(a,i8)') ' The number of iterations taken is ITERATIONS = ', iterations
Write (, '(a,g14.6)') ' The threshhold THRESH = ', thresh
Write (, '(a,i12)') ' The seed SEED = ', seed
Write (, '(a)') ' '
Write (, '(a)') ' The transition probability table, d on the number of'
Write (, '(a)') ' neighbors with the same spin.'
Write (, '(a)') ' '
Write (, '(a)') ' 1 2 3 4 5'
Write (, '(a)') ' '
Write (, '(7f10.4)')(prob(i), i=1, 5)
!
! Do the simulation.
! I have to do this in a subroutine, so that I can dimension C1 on the fly.
!
Call handle(m, n, iterations, prob, thresh, seed)
!
! Terminate.
!
Write (, '(a)') ' '
Write (, '(a)') 'ISING_2D_SIMULATION'
Write (, '(a)') ' Normal end of execution.'

Write (*, '(a)') ' '
Call timestamp()

Stop
End Program main
Subroutine ch_cap(ch)

!*********************************************************************72
!
!c CH_CAP capitalizes a single character.
!
! Parameters:
!
! Input/output, character CH, the character to capitalize.
!
Implicit None

Character ch
Integer itemp

itemp = ichar(ch)

If (97<=itemp .And. itemp<=122) Then
ch = char(itemp-32)
End If

Return
End Subroutine ch_cap
Function ch_eqi(c1, c2)

!*********************************************************************72
!
!c CH_EQI is a case insensitive comparison of two characters for equality.
!
! Example:
!
! CH_EQI ( 'A', 'a' ) is TRUE.
!
! Parameters:
!
! Input, character C1, C2, the characters to compare.
!
! Output, logical CH_EQI, the result of the comparison.
!
Implicit None

Character c1
Character c1_cap
Character c2
Character c2_cap
Logical ch_eqi

c1_cap = c1
c2_cap = c2

Call ch_cap(c1_cap)
Call ch_cap(c2_cap)

If (c1_cap==c2_cap) Then
ch_eqi = .True.
Else
ch_eqi = .False.
End If

Return
End Function ch_eqi
Subroutine ch_to_digit(c, digit)

!*********************************************************************72
!
!c CH_TO_DIGIT returns the integer value of a 10 digit.
!
! Example:
!
! C DIGIT
! --- -----
! '0' 0
! '1' 1
! ... ...
! '9' 9
! ' ' 0
! 'X' -1
!
! Parameters:
!
! Input, character C, the decimal digit, '0' through '9' or blank
! are legal.
!
! Output, integer DIGIT, the corresponding integer value. If C was
! 'illegal', then DIGIT is -1.
!
Implicit None

Character c
Integer digit

If (lge(c,'0') .And. lle(c,'9')) Then

digit = ichar(c) - 48

Else If (c==' ') Then

digit = 0

Else

digit = -1

End If

Return
End Subroutine ch_to_digit
Subroutine get_input(m, n, iterations, thresh, seed)

!*********************************************************************72
!
!c GET_INPUT gets input parameters.
!
! Parameters:
!
! Output, integer M, N, the number of rows and columns.
!
! Output, integer ITERATIONS, the number of iterations.
!
! Output, double precision THRESH, the threshhold.
!
! Output, integer SEED, a seed for the random
! number generator.
!
Implicit None

Integer arg_num
Integer iarg
Integer iargc
Integer ierror
Integer iterations
Integer last
Integer m
Integer n
Integer seed
Character *(80) string
Double Precision thresh

arg_num = iargc()

If (1<=arg_num) Then
iarg = 1
Call getarg(iarg, string)
Call s_to_i4(string, m, ierror, last)
Else If (.True.) Then
m = 15
Else
Write (, '(a)') ' '
Write (, '(a)') ' Enter M, the number of rows:'
Read (*, *) m
End If

If (2<=arg_num) Then
iarg = 2
Call getarg(iarg, string)
Call s_to_i4(string, n, ierror, last)
Else If (.True.) Then
n = 20
Else
Write (, '(a)') ' '
Write (, '(a)') ' Enter N, the number of columns:'
Read (*, *) n
End If

If (3<=arg_num) Then
iarg = 3
Call getarg(iarg, string)
Call s_to_i4(string, iterations, ierror, last)
Else If (.True.) Then
iterations = 35
Else
Write (, '(a)') ' '
Write (, '(a)') ' Enter the number of iterations to take.'
Read (*, *) iterations
End If

If (4<=arg_num) Then
iarg = 4
Call getarg(iarg, string)
Call s_to_r8(string, thresh, ierror, last)
Else If (.True.) Then
thresh = 0.50D+00
Else
Write (, '(a)') ' '
Write (, '(a)') ' Enter the threshhold.'
Read (*, *) thresh
End If

If (5<=arg_num) Then
iarg = 5
Call getarg(iarg, string)
Call s_to_i4(string, seed, ierror, last)
Else If (.True.) Then
seed = 123456789
Else
Write (, '(a)') ' '
Write (, '(a)') ' Enter the random number seed.'
Read (*, *) seed
End If

Return
End Subroutine get_input
Subroutine get_unit(iunit)

!*********************************************************************72
!
!c GET_UNIT returns a free FORTRAN unit number.
!
! Discussion:
!
! A "free" FORTRAN unit number is a value between 1 and 99 which
! is not currently associated with an I/O device. A free FORTRAN unit
! number is needed in order to open a file with the OPEN command.
!
! If IUNIT = 0, then no free FORTRAN unit could be found, although
! all 99 units were checked (except for units 5, 6 and 9, which
! are commonly reserved for console I/O).
!
! Otherwise, IUNIT is a value between 1 and 99, representing a
! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6
! are special, and will never return those values.
!
! Parameters:
!
! Output, integer IUNIT, the free unit number.
!
Implicit None

Integer i
Integer iunit
Logical value

iunit = 0

Do i = 1, 99

If (i/=5 .And. i/=6 .And. i/=9) Then

  Inquire (Unit=i, Opened=value, Err=10)

  If (.Not. value) Then
    iunit = i
    Return
  End If

End If

10 Continue

End Do

Return
End Subroutine get_unit
Subroutine handle(m, n, iterations, prob, thresh, seed)

!*********************************************************************72
!
!c HANDLE does initialization, computation and output.
!
! Discussion:
!
! The main program picks up the dimensions M and N, so we have to
! call this subroutine so we can essentially allocate C1 on the fly.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns.
!
! Input, integer ITERATIONS, the number of iterations.
!
! Input, double precision PROB(1:5). PROB(I) represents the probability
! that the spin of a given cell will be reversed, given that it has I
! immediate neighbors (including itself) with spin the same as its own.
!
! Input, double precision THRESH, the threshhold.
!
! Input/output, integer SEED, a seed for the random number
! generator.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer iterations
Double Precision prob(1:5)
Integer seed
Double Precision thresh
!
! Initialize the system.
!
Call ising_2d_initialize(m, n, thresh, seed, c1)
!
! Write the initial state to a gnuplot file.
!
Call plot_file(m, n, c1, 'Initial Configuration', 'ising_2d_initial.txt', &
'ising_2d_initial.png')
!
! Do the simulation.
!
Call transition(m, n, iterations, prob, seed, c1)
!
! Write the final state to a gnuplot file.
!
Call plot_file(m, n, c1, 'Final Configuration', 'ising_2d_final.txt',&
'ising_2d_final.png')

Return
End Subroutine handle
Function i4_modp(i, j)

!*********************************************************************72
!
!c I4_MODP returns the nonnegative remainder of integer division.
!
! Discussion:
!
! If
! NREM = I4_MODP ( I, J )
! NMULT = ( I - NREM ) / J
! then
! I = J * NMULT + NREM
! where NREM is always nonnegative.
!
! The MOD function computes a result with the same sign as the
! quantity being divided. Thus, suppose you had an angle A,
! and you wanted to ensure that it was between 0 and 360.
! Then mod(A,360) would do, if A was positive, but if A
! was negative, your result would be between -360 and 0.
!
! On the other hand, I4_MODP(A,360) is between 0 and 360, always.
!
! Example:
!
! I J MOD I4_MODP Factorization
!
! 107 50 7 7 107 = 2 * 50 + 7
! 107 -50 7 7 107 = -2 * -50 + 7
! -107 50 -7 43 -107 = -3 * 50 + 43
! -107 -50 -7 43 -107 = 3 * -50 + 43
!
! Parameters:
!
! Input, integer I, the number to be divided.
!
! Input, integer J, the number that divides I.
!
! Output, integer I4_MODP, the nonnegative remainder when I is
! divided by J.
!
Implicit None

Integer i
Integer i4_modp
Integer j
Integer value

If (j==0) Then
Write (, '(a)') ' '
Write (, '(a)') 'I4_MODP - Fatal error!'
Write (*, '(a,i8)') ' Illegal divisor J = ', j
Stop
End If

value = mod(i, j)

If (value<0) Then
value = value + abs(j)
End If

i4_modp = value

Return
End Function i4_modp
Function i4_wrap(ival, ilo, ihi)

!*********************************************************************72
!
!c I4_WRAP forces an I4 to lie between given limits by wrapping.
!
! Example:
!
! ILO = 4, IHI = 8
!
! I Value
!
! -2 8
! -1 4
! 0 5
! 1 6
! 2 7
! 3 8
! 4 4
! 5 5
! 6 6
! 7 7
! 8 8
! 9 4
! 10 5
! 11 6
! 12 7
! 13 8
! 14 4
!
! Parameters:
!
! Input, integer IVAL, an integer value.
!
! Input, integer ILO, IHI, the desired bounds for the integer value.
!
! Output, integer I4_WRAP, a "wrapped" version of IVAL.
!
Implicit None

Integer i4_modp
Integer i4_wrap
Integer ihi
Integer ilo
Integer ival
Integer jhi
Integer jlo
Integer value
Integer wide

jlo = min(ilo, ihi)
jhi = max(ilo, ihi)

wide = jhi - jlo + 1

If (wide==1) Then
value = jlo
Else
value = jlo + i4_modp(ival-jlo, wide)
End If

i4_wrap = value

Return
End Function i4_wrap
Subroutine ising_2d_agree(m, n, c1, c5)

!*********************************************************************72
!
!c ISING_2D_AGREE returns the number of neighbors agreeing with each cell.
!
! Discussion:
!
! The count includes the cell itself, so it is between 1 and 5.
!
! Parameters:
!
! Input, integer M, N, the number of cells in each
! spatial dimension.
!
! Input, integer C1(M,N), an array of 1's and -1's.
!
! Output, integer C5(M,N), the number of neighbors
! that agree. 1, 2, 3, 4, or 5.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer c5(m, n)
Integer i
Integer i4_wrap
Integer im
Integer ip
Integer j
Integer jm
Integer jp

Do j = 1, n
jp = i4_wrap(j+1, 1, n)
jm = i4_wrap(j-1, 1, n)
Do i = 1, m
ip = i4_wrap(i+1, 1, m)
im = i4_wrap(i-1, 1, m)
c5(i, j) = c1(i, j) + c1(ip, j) + c1(im, j) + c1(i, jm) + c1(i, jp)
End Do
End Do

Do j = 1, n
Do i = 1, m
If (0<c1(i,j)) Then
c5(i, j) = (5+c5(i,j))/2
Else
c5(i, j) = (5-c5(i,j))/2
End If
End Do
End Do

Return
End Subroutine ising_2d_agree
Subroutine ising_2d_initialize(m, n, thresh, seed, c1)

!*********************************************************************72
!
!c ISING_2D_INITIALIZE initializes the Ising array.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns.
!
! Input, double precision THRESH, the threshhold.
!
! Input/output, integer SEED, a seed for the random
! number generator.
!
! Output, integer C1(M,N), the initial Ising array.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer i
Integer j
Double Precision r(m, n)
Double Precision thresh
Integer seed

Call r8mat_uniform_01(m, n, seed, r)

Do j = 1, n
Do i = 1, m
If (r(i,j)<=thresh) Then
c1(i, j) = -1
Else
c1(i, j) = +1
End If
End Do
End Do

Return
End Subroutine ising_2d_initialize
Subroutine ising_2d_stats(step, m, n, c1)

!*********************************************************************72
!
!c ISING_2D_STATS prints information about the current step.
!
! Parameters:
!
! Input, integer STEP, the step number.
!
! Input, integer M, N, the number of rows and columns.
!
! Input, integer C1(M,N), the current state of the system.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer i
Integer j
Integer pos_count
Double Precision pos_percent
Integer step
Integer neg_count
Double Precision neg_percent

If (step==0) Then
Write (, '(a)') ' '
Write (, '(a)') ' Step Positives Negatives'
Write (, '(a)') ' # % # %'
Write (, '(a)') ' '
End If

pos_count = 0
Do j = 1, n
Do i = 1, m
If (0<c1(i,j)) Then
pos_count = pos_count + 1
End If
End Do
End Do

neg_count = mn - pos_count
pos_percent = dble(100pos_count)/dble(mn)
neg_percent = dble(100neg_count)/dble(m*n)

Write (*, '(2x,i4,2x,i6,2x,f7.3,2x,i6,2x,f7.3)') step, pos_count, &
pos_percent, neg_count, neg_percent

Return
End Subroutine ising_2d_stats
Subroutine neighbor_2d_stats(step, m, n, c1, c5)

!*********************************************************************72
!
!c NEIGHBOR_2D_STATS prints neighbor statistics about the current step.
!
! Parameters:
!
! Input, integer STEP, the step number.
!
! Input, integer M, N, the number of rows and columns.
!
! Input, integer C1(M,N), the current state of the system.
!
! Input, integer C5(M,N), the number of agreeable neighbors.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer c5(m, n)
Integer i
Integer j
Integer k
Integer stats(-5:5)
Integer step

If (step==0) Then
Write (, '(a)') ' '
Write (, '(a)') ' Step Neighborhood Charge:'
Write (, '(a)') ' -5 -4 -3 -2 -1' //&
' +1 +2 +3 +4 +5'
Write (, '(a)') ' '
End If

Do i = -5, 5
stats(i) = 0
End Do

Do j = 1, n
Do i = 1, m
stats(c5(i,j)) = stats(c5(i,j)) + 1
End Do
End Do
Write (*, '(2x,i4,10(2x,i4))') step, (stats(i), i=-5, -1), (stats(i), i=1, 5)

Return
End Subroutine neighbor_2d_stats
Subroutine plot_file(m, n, c1,title, plot_filename, png_filename)

!*********************************************************************72
!
!c PLOT_FILE writes the current configuration to a GNUPLOT plot file.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns.
!
! Input, integer C1(M,N), the current state of the system.
!
! Input, character * ( * ) , a for the plot.
!
! Input, character * ( * ) PLOT_FILENAME, a name for the GNUPLOT
! command file to be created.
!
! Input, character * ( * ) PNG_FILENAME, the name of the PNG graphics
! file to be created.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer i
Integer j
Character () plot_filename
Integer plot_unit
Character () png_filename
Double Precision ratio
Character () title
Integer x1
Integer x2
Integer y1
Integer y2

Call get_unit(plot_unit)

Open (Unit=plot_unit, File=plot_filename, Status='replace')

ratio = dble(n)/dble(m)

Write (plot_unit, '(a)') 'set term png'
Write (plot_unit, '(a)') 'set output "' // trim(png_filename) // '"'
Write (plot_unit, '(a,i4,a)') 'set xrange [ 0 : ', m, ' ]'
Write (plot_unit, '(a,i4,a)') 'set yrange [ 0 : ', n, ' ]'
Write (plot_unit, '(a)') 'set nokey'
Write (plot_unit, '(a)') 'set title"' // trim(title ) // '"'
Write (plot_unit, '(a)') 'unset tics'

Write (plot_unit, '(a,g15.7)') 'set size ratio ', ratio
Do j = 1, n
y1 = j - 1
y2 = j
Do i = 1, m
x1 = m - i
x2 = m + 1 - i
If (c1(i,j)<0) Then
Write (plot_unit, '(a,i3,a,i3,a,i3,a,i3,a)') &
'set rectangle from ', x1, ',', y1, ' to ', x2, ',', y2, &
' fc rgb "blue"'
Else
Write (plot_unit, '(a,i3,a,i3,a,i3,a,i3,a)')&
'set rectangle from ', x1, ',', y1, ' to ', x2, ',', y2,&
' fc rgb "red"'
End If
End Do
End Do

Write (plot_unit, '(a)') 'plot 1'
Write (plot_unit, '(a)') 'quit'

Close (Unit=plot_unit)

Write (, '(a)') ' '
Write (, '(a)') ' Created the gnuplot graphics file "' // &
trim(plot_filename) // '".'

Return
End Subroutine plot_file
Subroutine r8mat_uniform_01(m, n, seed, r)

!*********************************************************************72
!
!c R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT.
!
! Discussion:
!
! An R8MAT is an array of R8's.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns in the array.
!
! Input/output, integer SEED, the "seed" value, which should NOT be 0.
! On output, SEED has been updated.
!
! Output, double precision R(M,N), the array of pseudorandom values.
!
Implicit None

Integer m
Integer n

Integer i
Integer j
Integer k
Integer seed
Double Precision r(m, n)

If (seed==0) Then
Write (, '(a)') ' '
Write (, '(a)') 'R8MAT_UNIFORM_01 - Fatal error!'
Write (*, '(a)') ' Input value of SEED = 0.'
Stop
End If

Do j = 1, n

Do i = 1, m

  k = seed/127773

  seed = 16807*(seed-k*127773) - k*2836

  If (seed<0) Then
    seed = seed + 2147483647
  End If

  r(i, j) = dble(seed)*4.656612875D-10

End Do

End Do

Return
End Subroutine r8mat_uniform_01
Function s_len_trim(s)

!*******************************************************************72
!
!c S_LEN_TRIM returns the length of a string to the last nonblank.
!
! Parameters:
!
! Input, character() S, a string.
!
! Output, integer S_LEN_TRIM, the length of the string to the last nonblank.
!
Implicit None

Integer i
Character () s
Integer s_len_trim

Do i = len(s), 1, -1

If (s(i:i)/=' ') Then
  s_len_trim = i
  Return
End If

End Do

s_len_trim = 0

Return
End Function s_len_trim
Subroutine s_to_i4(s, ival, ierror, length)

!*********************************************************************72
!
!c S_TO_I4 reads an I4 from a string.
!
! Parameters:
!
! Input, character * ( * ) S, a string to be examined.
!
! Output, integer IVAL, the integer value read from the string.
! If the string is blank, then IVAL will be returned 0.
!
! Output, integer IERROR, an error flag.
! 0, no error.
! 1, an error occurred.
!
! Output, integer LENGTH, the number of characters of S
! used to make IVAL.
!
Implicit None

Character c
Integer i
Integer ierror
Integer isgn
Integer istate
Integer ival
Integer length
Character () s
Integer s_len_trim

ierror = 0
istate = 0
isgn = 1
ival = 0

Do i = 1, s_len_trim(s)

c = s(i:i)

!
! Haven't read anything.
!
If (istate==0) Then

If (c==' ') Then

  Else If (c=='-') Then
    istate = 1
    isgn = -1
  Else If (c=='+') Then
    istate = 1
    isgn = +1
  Else If (lle('0',c) .And. lle(c,'9')) Then
    istate = 2
    ival = ichar(c) - ichar('0')
  Else
    ierror = 1
    Return
  End If

!
! Have read the sign, expecting digits.
!
Else If (istate==1) Then

If (c==' ') Then

  Else If (lle('0',c) .And. lle(c,'9')) Then
    istate = 2
    ival = ichar(c) - ichar('0')
  Else
    ierror = 1
    Return
  End If

!
! Have read at least one digit, expecting more.
!
Else If (istate==2) Then

If (lle('0',c) .And. lle(c,'9')) Then
    ival = 10*ival + ichar(c) - ichar('0')
  Else
    ival = isgn*ival
    length = i - 1
    Return
  End If

End If

End Do
!
! If we read all the characters in the string, see if we're OK.
!
If (istate==2) Then
ival = isgn*ival
length = s_len_trim(s)
Else
ierror = 1
length = 0
End If

Return
End Subroutine s_to_i4
Subroutine s_to_r8(s, dval, ierror, length)

!*********************************************************************72
!
!c S_TO_R8 reads an R8 from a string.
!
! Discussion:
!
! The routine will read as many characters as possible until it reaches
! the end of the string, or encounters a character which cannot be
! part of the number.
!
! Legal input is:
!
! 1 blanks,
! 2 '+' or '-' sign,
! 2.5 blanks
! 3 integer part,
! 4 decimal point,
! 5 fraction part,
! 6 'E' or 'e' or 'D' or 'd', exponent marker,
! 7 exponent sign,
! 8 exponent integer part,
! 9 exponent decimal point,
! 10 exponent fraction part,
! 11 blanks,
! 12 final comma or semicolon,
!
! with most quantities optional.
!
! Example:
!
! S DVAL
!
! '1' 1.0
! ' 1 ' 1.0
! '1A' 1.0
! '12,34,56' 12.0
! ' 34 7' 34.0
! '-1E2ABCD' -100.0
! '-1X2ABCD' -1.0
! ' 2E-1' 0.2
! '23.45' 23.45
! '-4.2E+2' -420.0
! '17d2' 1700.0
! '-14e-2' -0.14
! 'e2' 100.0
! '-12.73e-9.23' -12.73 * 10.0^(-9.23)
!
! Parameters:
!
! Input, character * ( * ) S, the string containing the
! data to be read. Reading will begin at position 1 and
! terminate at the end of the string, or when no more
! characters can be read to form a legal real. Blanks,
! commas, or other nonnumeric data will, in particular,
! cause the conversion to halt.
!
! Output, double precision DVAL, the value read from the string.
!
! Output, integer IERROR, error flag.
! 0, no errors occurred.
! 1, 2, 6 or 7, the input number was garbled. The
! value of IERROR is the last type of input successfully
! read. For instance, 1 means initial blanks, 2 means
! a plus or minus sign, and so on.
!
! Output, integer LENGTH, the number of characters read
! to form the number, including any terminating
! characters such as a trailing comma or blanks.
!
Implicit None

Logical ch_eqi
Character c
Double Precision dval
Integer ierror
Integer ihave
Integer isgn
Integer iterm
Integer jbot
Integer jsgn
Integer jtop
Integer length
Integer nchar
Integer ndig
Double Precision rbot
Double Precision rexp
Double Precision rtop
Character () s
Integer s_len_trim

nchar = s_len_trim(s)

ierror = 0
dval = 0.0D+00
length = -1
isgn = 1
rtop = 0
rbot = 1
jsgn = 1
jtop = 0
jbot = 1
ihave = 1
iterm = 0

10 Continue

length = length + 1

If (nchar<length+1) Then
Goto 20
End If

c = s(length+1:length+1)
!
! Blank character.
!
If (c==' ') Then

If (ihave==2) Then

Else If (ihave==6 .Or. ihave==7) Then
  iterm = 1
Else If (1<ihave) Then
  ihave = 11
End If

!
! Comma.
!
Else If (c==',' .Or. c==';') Then

If (ihave/=1) Then
  iterm = 1
  ihave = 12
  length = length + 1
End If

!
! Minus sign.
!
Else If (c=='-') Then

If (ihave==1) Then
  ihave = 2
  isgn = -1
Else If (ihave==6) Then
  ihave = 7
  jsgn = -1
Else
  iterm = 1
End If

!
! Plus sign.
!
Else If (c=='+') Then

If (ihave==1) Then
  ihave = 2
Else If (ihave==6) Then
  ihave = 7
Else
  iterm = 1
End If

!
! Decimal point.
!
Else If (c=='.') Then

If (ihave<4) Then
  ihave = 4
Else If (6<=ihave .And. ihave<=8) Then
  ihave = 9
Else
  iterm = 1
End If

!
! Scientific notation exponent marker.
!
Else If (ch_eqi(c,'E') .Or. ch_eqi(c,'D')) Then

If (ihave<6) Then
  ihave = 6
Else
  iterm = 1
End If

!
! Digit.
!
Else If (ihave<11 .And. lle('0',c) .And. lle(c,'9')) Then

If (ihave<=2) Then
  ihave = 3
Else If (ihave==4) Then
  ihave = 5
Else If (ihave==6 .Or. ihave==7) Then
  ihave = 8
Else If (ihave==9) Then
  ihave = 10
End If

Call ch_to_digit(c, ndig)

If (ihave==3) Then
  rtop = 10.0D+00*rtop + dble(ndig)
Else If (ihave==5) Then
  rtop = 10.0D+00*rtop + dble(ndig)
  rbot = 10.0D+00*rbot
Else If (ihave==8) Then
  jtop = 10*jtop + ndig
Else If (ihave==10) Then
  jtop = 10*jtop + ndig
  jbot = 10*jbot
End If

!
! Anything else is regarded as a terminator.
!
Else
iterm = 1
End If
!
! If we haven't seen a terminator, and we haven't examined the
! entire string, go get the next character.
!
If (iterm==1) Then
Goto 20
End If

Goto 10

20 Continue
!
! If we haven't seen a terminator, and we have examined the
! entire string, then we're done, and LENGTH is equal to NCHAR.
!
If (iterm/=1 .And. length+1==nchar) Then
length = nchar
End If
!
! Number seems to have terminated. Have we got a legal number?
! Not if we terminated in states 1, 2, 6 or 7.
!
If (ihave==1 .Or. ihave==2 .Or. ihave==6 .Or. ihave==7) Then
ierror = ihave
Write (*, '(a)') ' '
Write (*, '(a)') 'S_TO_R8 - Serious error!'
Write (*, '(a)') ' Illegal or nonnumeric input:'
Write (*, '(a,a)') ' ', s
Return
End If
!
! Number seems OK. Form it.
!
If (jtop==0) Then
rexp = 1.0D+00
Else
If (jbot==1) Then
rexp = 10.0D+00**(jsgn*jtop)
Else
rexp = 10.0D+00**(dble(jsgn*jtop)/dble(jbot))
End If
End If

dval = dble(isgn)rexprtop/rbot

Return
End Subroutine s_to_r8
Subroutine timestamp()

!*********************************************************************72
!
!c TIMESTAMP prints out the current YMDHMS date as a timestamp.
!
! Discussion:
!
! This FORTRAN version is made available for cases where the
! FORTRAN90 version cannot be used.
!
! Parameters:
!
! None
!
Implicit None

Character *(8) ampm
Integer d
Character *(8) date
Integer h
Integer m
Integer mm
Character *(9) month(12)
Integer n
Integer s
Character *(10) time
Integer y

Save month

Data month/'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', 'September',&
'October ', 'November ', 'December '/

Call date_and_time(date, time)

Read (date, '(i4,i2,i2)') y, m, d
Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm

If (h<12) Then
ampm = 'AM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Noon'
Else
ampm = 'PM'
End If
Else
h = h - 12
If (h<12) Then
ampm = 'PM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Midnight'
Else
ampm = 'AM'
End If
End If
End If

Write (*, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d, &
month(m), y, h, ':', n, ':', s, '.', mm, ampm

Return
End Subroutine timestamp
Subroutine transition(m, n, iterations, prob, seed, c1)

!*********************************************************************72
!
!c TRANSITION carries out a Monte Carlo simulation of a 3D Ising model.
!
! Parameters:
!
! Input, integer M, N, the number of rows and columns.
!
! Input, integer ITERATIONS, the number of iterations.
!
! Input, double precision PROB(1:5). PROB(I) represents the probability
! that the spin of a given cell will be reversed, given that it has I
! immediate neighbors (including itself) with spin the same as its own.
!
! Input/output, integer SEED, a seed for the random number
! generator.
!
! Input/output, integer C1(M,N). On input, the current
! state of the system. On output, the state of the system after
! the iterations.
!
Implicit None

Integer m
Integer n

Integer c1(m, n)
Integer c5(m, n)
Integer i
Integer iterations
Integer j
Double Precision prob(1:5)
Integer step
Double Precision r(m, n)
Integer seed
Double Precision thresh

step = 0
Call ising_2d_stats(step, m, n, c1)

Do step = 1, iterations
!
! C5 contains 1 through 5, the number of cells that agree with the center cell.
!
Call ising_2d_agree(m, n, c1, c5)

If (.False.) Then
  Call neighbor_2d_stats(step, m, n, c1, c5)
End If

!
! Determine the chances of flipping cell (I,J).
!
Call r8mat_uniform_01(m, n, seed, r)

Do j = 1, n
  Do i = 1, m
    If (r(i,j)<prob(c5(i,j))) Then
      c1(i, j) = -c1(i, j)
    End If
  End Do
End Do

Call ising_2d_stats(step, m, n, c1)

End Do

Return
End Subroutine transition

!c$$Run, Output:

7 September 2022 7:00:19.714 AM ISING_2D_SIMULATION FORTRAN95,
2013&18 version Monte Carlo simulation of a 2D Ising model.
The number of rows is M = 15
The number of columns is N = 20
The number of iterations taken is ITERATIONS = 35
The threshhold THRESH = 0.500000
The seed SEED = 123456789

The transition probability table, d on the number of
neighbors with the same spin.

1         2         3         4         5

0.9700    0.8600    0.5100    0.1600    0.0300

Created the gnuplot graphics file "ising_2d_initial.txt".

Step Positives Negatives

% # %

0     147   49.000     153   51.000
 1     147   49.000     153   51.000
 2     147   49.000     153   51.000
 3     147   49.000     153   51.000
 4     151   50.333     149   49.667
 5     152   50.667     148   49.333
 6     155   51.667     145   48.333
 7     153   51.000     147   49.000
 8     164   54.667     136   45.333
 9     163   54.333     137   45.667
10     160   53.333     140   46.667
11     167   55.667     133   44.333
12     165   55.000     135   45.000
13     169   56.333     131   43.667
14     162   54.000     138   46.000
15     165   55.000     135   45.000
16     167   55.667     133   44.333
17     172   57.333     128   42.667
18     167   55.667     133   44.333
19     160   53.333     140   46.667
20     156   52.000     144   48.000
21     158   52.667     142   47.333
22     163   54.333     137   45.667
23     155   51.667     145   48.333
24     161   53.667     139   46.333
25     163   54.333     137   45.667
26     143   47.667     157   52.333
27     144   48.000     156   52.000
28     137   45.667     163   54.333
29     140   46.667     160   53.333
30     137   45.667     163   54.333
31     136   45.333     164   54.667
32     140   46.667     160   53.333
33     144   48.000     156   52.000
34     144   48.000     156   52.000
35     144   48.000     156   52.000

Created the gnuplot graphics file "ising_2d_final.txt".

ISING_2D_SIMULATION
Normal end of execution.
April 2022
Program finished with exit code 0
Press ENTER to exit console.

set term png
set output "ising_2d_final.png"
set xrange [ 0 : 15 ]
set yrange [ 0 : 20 ]
set nokey
set "Final Configuration"
unset tics
set size ratio 1.333333
set rectangle from 14, 0 to 15, 1 fc rgb "blue"
set rectangle from 13, 0 to 14, 1 fc rgb "blue"
set rectangle from 12, 0 to 13, 1 fc rgb "blue"
......................................
set rectangle from 3, 19 to 4, 20 fc rgb "blue"
set rectangle from 2, 19 to 3, 20 fc rgb "blue"
set rectangle from 1, 19 to 2, 20 fc rgb "blue"
set rectangle from 0, 19 to 1, 20 fc rgb "red"
plot 1
quit

set term png
set output "ising_2d_initial.png"
set xrange [ 0 : 15 ]
set yrange [ 0 : 20 ]
set nokey
set "Initial Configuration"
unset tics
set size ratio 1.333333
set rectangle from 14, 0 to 15, 1 fc rgb "blue"
set rectangle from 13, 0 to 14, 1 fc rgb "red"
set rectangle from 12, 0 to 13, 1 fc rgb "red"
set rectangle from 11, 0 to 12, 1 fc rgb "red"
set rectangle from 10, 0 to 11, 1 fc rgb "blue"
set rectangle from 9, 0 to 10, 1 fc rgb "blue"
.......................................
set rectangle from 2, 19 to 3, 20 fc rgb "blue"
set rectangle from 1, 19 to 2, 20 fc rgb "red"
set rectangle from 0, 19 to 1, 20 fc rgb "blue"
plot 1
quit
Created the gnuplot graphics file "ising_2d_final.txt".
ISING_2D_SIMULATION Normal end of execution.
7 September 2022 7:00:19.716 AM

易辛二维度模拟正负电荷的2D阵列进化每个很会“翻转”与邻居一致,
程序运行验证演示数值如上
及图形处理命令，若按其执行应该可以会获得图形表达示意.

广州