Three level spin system

globals
[
  n_zero_a           n_zero_a_update             ; number of patches in the middle energy state in spin group a
  n_one_a            n_one_a_update              ; number of patches in the lowest energy state in spin group a
  n_two_a            n_two_a_update              ; number of patches in the highest energy state in spin group a

  n_zero_b           n_zero_b_update             ; number of patches in the middle energy state in spin group b
  n_one_b            n_one_b_update              ; number of patches in the lowest energy state in spin group b
  n_two_b            n_two_b_update              ; number of patches in the highest energy state in spin group b

  N_total                                        ; total number of patches
  N_total_a          N_total_b                   ; total number of patches in each spin group
  N_total_a_update   N_total_b_update            ; verification variables for total population

  A                                              ; orientation of a selected patch
  neighbor1                                      ; neighboring patch of patch A
  neighbor2                                      ; neighboring patch of patch A
  B1                                             ; orientation of neighbor1
  B2                                             ; orientation of neighbor2
  tempa_hold                                     ; temperature of spin group a of the previous tick
  tempb_hold                                     ; temperature of spin group b of the previous tick
]

patches-own[orientation group friend?]           ; orientation is based on the energy level of the patch
                                                 ; group is which spin group a patch is in a/b (blue/pink)
                                                 ; friend? is used in selecting neighbors, used in "update1"

to setup
  ca                                             ; clears all previous data from previous runs
  if world_division = false                      ; calls a specific "setup" function based on the configuration of the switches
  [
    ifelse dense_setup = false [setup1] [setup3]
  ]
  if world_division = true
  [
    ifelse dense_setup = true [setup4] [setup2]
  ]
  ask patches [recolor set friend? false]        ; calls every patch to run through the recolor function and has default friend? = false

  set n_two_a_update n_two_a                     ; initial conditions for updated n values
  set n_one_a_update n_one_a                     ; updated n values are used to recalculate the thermal properties in real time as seen in monitors
  set n_zero_a_update n_zero_a

  set n_two_b_update n_two_b
  set n_one_b_update n_one_b
  set n_zero_b_update n_zero_b
end 

to setup1  ; world divison false and dense setup false
  ca                                                                                        ; clears all previous data from previous runs
  reset-ticks                                                                               ; resets the tick counter from any previous runs

  if world_division = false
  [
    set N_total (world-width * world-height)                                                ; counts the number of patches in the entire world
    set N_total_a (0.5 * world-width * world-height)                                        ; spin groups a/b are each 1/2 of the world split vertically
    set N_total_b (0.5 * world-width * world-height)

    set n_two_a floor((percent_n_two_a / 100) * (0.5 * world-width * world-height))         ; populates each state based on the observers inputed slider values
    set n_one_a floor((percent_n_one_a / 100) * (0.5 * world-width * world-height))
    set n_two_b floor((percent_n_two_b / 100) * (0.5 * world-width * world-height))
    set n_one_b floor((percent_n_one_b / 100) * (0.5 * world-width * world-height))

    set n_zero_a  N_total_a - (n_one_a + n_two_a)                                           ; by default, n zero is all patches not assigned from sliders (n one and n two)
    set n_zero_b  N_total_b - (n_one_b + n_two_b)

    ask patches                                                                             ; patches are populated with initial color, orientation and group #
     [
              ifelse pxcor > (world-width / 2 - 1)                                          ; spin group a is also refered to as spin group 1 in code BLUE
                                                                                            ; spin group b is also refered to as spin group 2 in code PINK
                     [set pcolor pink set orientation 0 set group 2 ]
                     [set pcolor blue set orientation 0 set group 1 ]                       ; default orientation is 0 for every patch
     ]
  ]

  if dense_setup = false
  [
    ask n-of n_two_a patches with [group = 1] [set orientation 1]                           ; all n two patches have orientation = 1
    ask n-of n_one_a patches with [group = 1 and orientation != 1] [set orientation -1]     ; all n one patches have orientation = -1
                                                                                            ; note : all n zero patches haave orientation = 0 by default
    ask n-of n_two_b patches with [group = 2] [set orientation 1]
    ask n-of n_one_b patches with [group = 2 and orientation != 1] [set orientation -1]
  ]
end 

to setup2  ; world divsion true dense setup false
  ca                                                                                        ; clears all previous data from previous runs
  reset-ticks                                                                               ; resets the tick counter from any previous runs

  if world_division = true
  [
    set N_total (world-width * world-height)                                                ; counts the number of patches in the entire world
    set N_total_a (0.75 * world-width * world-height)                                       ; spin group a will now occupy 3/4 of the total patches in the world
    set N_total_b (0.25 * world-width * world-height)                                       ; spin group b will now occupy 1/4 of the total patches in the world

    set n_two_a floor((percent_n_two_a / 100) * (0.75 * world-width * world-height))        ; populates each state based on the observers inputed slider values
    set n_one_a floor((percent_n_one_a / 100) * (0.75 * world-width * world-height))
    set n_two_b floor((percent_n_two_b / 100) * (0.25 * world-width * world-height))
    set n_one_b floor((percent_n_one_b / 100) * (0.25 * world-width * world-height))

    set n_zero_a  N_total_a - (n_one_a + n_two_a)                                           ; by default, n zero is all patches not assigned from sliders (n one and n two)
    set n_zero_b  N_total_b - (n_one_b + n_two_b)

    ask patches                                                                             ; patches are populated with initial color, orientation and group #
    [
         ifelse pxcor > (world-width / 2 - 1) and (pycor < (world-height / 2 - 1))          ; spin group a is also refered to as spin group 1 in code BLUE
                                                                                            ; spin group b is also refered to as spin group 2 in code PINK
             [set pcolor pink set orientation 0 set group 2 ]
             [set pcolor blue set orientation 0 set group 1 ]
    ]
  ]

  if dense_setup = false
  [
    ask n-of n_two_a patches with [group = 1] [set orientation 1]                           ; all n two patches have orientation = 1
    ask n-of n_one_a patches with [group = 1 and orientation != 1] [set orientation -1]     ; all n one patches have orientation = -1
                                                                                            ; note : all n zero patches haave orientation = 0 by default
    ask n-of n_two_b patches with [group = 2] [set orientation 1]
    ask n-of n_one_b patches with [group = 2 and orientation != 1] [set orientation -1]
  ]
end 

to setup3 ; world divsion false and dense setup true
  ca                                                                                        ; clears all previous data from previous runs
  reset-ticks                                                                               ; resets the tick counter from any previous runs


  if world_division = false
  [
    set N_total (world-width * world-height)                                                ; counts the number of patches in the entire world
    set N_total_a (0.5 * world-width * world-height)                                        ; spin groups a/b are each 1/2 of the world split vertically
    set N_total_b (0.5 * world-width * world-height)

    set n_two_a floor((percent_n_two_a / 100) * (0.5 * world-width * world-height))         ; populates each state based on the observers inputed slider values
    set n_one_a floor((percent_n_one_a / 100) * (0.5 * world-width * world-height))
    set n_two_b floor((percent_n_two_b / 100) * (0.5 * world-width * world-height))
    set n_one_b floor((percent_n_one_b / 100) * (0.5 * world-width * world-height))

    set n_zero_a  N_total_a - (n_one_a + n_two_a)                                           ; by default, n zero is all patches not assigned from sliders (n one and n two)
    set n_zero_b  N_total_b - (n_one_b + n_two_b)

    ask patches                                                                             ; patches are populated with initial color, orientation and group #
     [
              ifelse pxcor > (world-width / 2 - 1)                                          ; spin group a is also refered to as spin group 1 in code BLUE
                                                                                            ; spin group b is also refered to as spin group 2 in code PINK
                     [set pcolor pink set orientation 0 set group 2 ]
                     [set pcolor blue set orientation 0 set group 1 ]                       ; default orientation is 0 for every patch
     ]
  ]

  if dense_setup = true                                                                     ; populates a concentrated group of patches on the outer most area of the world
                                                                                            ; creating a "sliver" of excited patches
  [
    let none_xa floor ((percent_n_one_a / 100) * (world-width))                             ; x variables are holding variables used to populate each energy state for this specific configuration
    let ntwo_xa floor ((percent_n_two_a / 100) * (world-width))
    ask n-of n_two_a patches with [group = 1 and pxcor <= ntwo_xa]  [set orientation 1]
    ask n-of n_one_a patches with [group = 1 and orientation = 0 and pxcor <= none_xa] [set orientation -1]

    let none_xb ceiling (world-width - (percent_n_one_b / 100) * (world-width))
    let ntwo_xb ceiling (world-width - (percent_n_two_b / 100) * (world-width))
    ask n-of n_one_b patches with [group = 2 and pxcor >= none_xb - 1] [set orientation -1]
    ask n-of n_two_b patches with [group = 2 and orientation = 0 and pxcor >= ntwo_xb - 1] [set orientation 1]
    ask patches [recolor]
  ]
end 

to setup4  ; world division true and dense setup true
  ca                                                                                        ; clears all previous data from previous runs
  reset-ticks                                                                               ; resets the tick counter from any previous runs

  if world_division = true
  [
    set N_total (world-width * world-height)                                                ; counts the number of patches in the entire world
    set N_total_a (0.75 * world-width * world-height)                                       ; spin group a will now occupy 3/4 of the total patches in the world
    set N_total_b (0.25 * world-width * world-height)                                       ; spin group b will now occupy 1/4 of the total patches in the world

    set n_two_a ceiling((percent_n_two_a / 100) * (0.75 * world-width * world-height))      ; populates each state based on the observers inputed slider values
    set n_one_a ceiling((percent_n_one_a / 100) * (0.75 * world-width * world-height))
    set n_two_b ceiling((percent_n_two_b / 100) * (0.25 * world-width * world-height))
    set n_one_b ceiling((percent_n_one_b / 100) * (0.25 * world-width * world-height))

    set n_zero_a  N_total_a - (n_one_a + n_two_a)                                           ; by default, n zero is all patches not assigned from sliders (n one and n two)
    set n_zero_b  N_total_b - (n_one_b + n_two_b)

    ask patches                                                                             ; patches are populated with initial color, orientation and group #
    [
         ifelse pxcor > (world-width / 2 - 1) and (pycor < (world-height / 2 - 1))          ; spin group a is also refered to as spin group 1 in code BLUE
                                                                                            ; spin group b is also refered to as spin group 2 in code PINK
             [set pcolor pink set orientation 0 set group 2 ]
             [set pcolor blue set orientation 0 set group 1 ]
    ]
  ]

  if dense_setup = true
  [
    let bw (world-width / 4)
    let bh (world-height / 4)
    let block_area floor (bw * bh)
    let a_blocks ceiling ( (n_two_a + n_one_a) / block_area)                                ; blocks are used to populate group a around group b in a circular fashion
    let b_blocks ceiling ( (n_two_b + n_one_b) / block_area)                                ; the number of blocks need is determinded by the sliders of percent n one and n two


    ; BLUE GROUP A
    if a_blocks = 1
    [
      ask n-of n_two_a patches with [pxcor < bw and pycor <= bh and group = 1] [set orientation 1]
      ask n-of n_one_a patches with [pxcor < bw and pycor <= bh and group = 1 and orientation = 0] [set orientation -1]
    ]
    if a_blocks = 2
    [
      ask n-of n_two_a patches with
      [group = 1 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)) or (pxcor < bw and pycor <= (bh * 3)))]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)))]
      [set orientation -1]
    ]
    if a_blocks = 3
    [
      ask n-of n_two_a patches with
      [group = 1 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)) or (pxcor < bw and pycor <= (bh * 3)))]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)) or (pxcor < bw and pycor <= (bh * 3)))]
      [set orientation -1]
    ]
    if a_blocks = 4
    [
      ask n-of n_two_a patches with
      [group = 1 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)) or (pxcor < bw and pycor <= (bh * 3)) or (pxcor < bw and pycor <= (bh * 4)))]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ((pxcor < bw and pycor <= bh) or (pxcor < bw and pycor <= (bh * 2)) or (pxcor < bw and pycor <= (bh * 3)) or (pxcor < bw and pycor <= (bh * 4)))]
      [set orientation -1]
    ]
    if a_blocks = 5
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 2) and pycor > (bh * 3)) )]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 2) and pycor > (bh * 3))) ]
      [set orientation -1]
    ]
    if a_blocks = 6
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 3) and pycor >= (bh * 3)) )]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 3) and pycor >= (bh * 3))) ]
      [set orientation -1]
    ]
    if a_blocks = 7
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) )]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3))) ]
      [set orientation -1]
    ]
    if a_blocks = 8
    [
      ask n-of n_two_a patches with
      [group = 1 and ((pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= bh) ) ]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ((pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= bh) ) ]
      [set orientation -1]
    ]
    if a_blocks = 9
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 2)) ) ]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 2)) ) ]
      [set orientation -1]
    ]
    if a_blocks = 10
    [
      ask n-of n_two_a patches with
      [group = 1 and ((pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 3)) ) ]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 3)) ) ]
      [set orientation -1]
    ]
     if a_blocks = 11
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 3)) or (pxcor < (3 * bw) and pycor >= (2 * bh) ) ) ]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor < (bh * 3)) or (pxcor < (3 * bw) and pycor >= (2 * bh) )  ) ]
      [set orientation -1]
    ]
    if a_blocks = 12
    [
      ask n-of n_two_a patches with
      [group = 1 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 3)) or (pxcor < (4 * bw) and pycor > (2 * bh) ) ) ]
      [set orientation 1]
      ask n-of n_one_a patches with
      [group = 1 and orientation = 0 and ( (pxcor < bw and pycor <= (bh * 4)) or (pxcor < (bw * 4) and pycor >= (bh * 3)) or (pxcor < (bw * 2) and pycor <= (bh * 3)) or (pxcor < (4 * bw) and pycor > (2 * bh) )  ) ]
      [set orientation -1]
    ]


    ;PINK GROUP B
    if b_blocks = 1
    [
      ask n-of n_two_b patches with [group = 2 and pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)]
      [set orientation 1]
      ask n-of n_one_b patches with [group = 2 and orientation = 0 and pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)]
      [set orientation -1]
    ]
    if b_blocks = 2
    [
      ask n-of n_two_b patches with [group = 2 and  ( (pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1) )  )]
      [set orientation 1]
      ask n-of n_one_b patches with [group = 2 and orientation = 0 and ( (pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1) )  )]
      [set orientation -1]
    ]
    if b_blocks = 3
    [
      ask n-of n_two_b patches with [group = 2 and  ( (pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1) )  or (pxcor > ((3 * bw) - 1) and pycor <= (2 * bh)))]
      [set orientation 1]
      ask n-of n_one_b patches with [group = 2 and orientation = 0 and ( (pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1) ) or (pxcor > ((3 * bw) - 1) and pycor <= (2 * bh)) )]
      [set orientation -1]
    ]
    if b_blocks = 4
    [
      ask n-of n_two_b patches with [group = 2 and ((pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((3 * bw) - 1) and pycor <= (2 * bh)) or (pxcor >= ((2 * bw)) and pycor <= (2 * bh)))]
      [set orientation 1]
      ask n-of n_one_b patches with [group = 2 and orientation = 0 and ((pxcor > ((3 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((2 * bw) - 1) and pycor <= (bh - 1)) or (pxcor > ((3 * bw) - 1) and pycor <= (2 * bh)) or (pxcor >= ((2 * bw)) and pycor <= (2 * bh)))]
      [set orientation -1]
    ]
    ask patches [recolor]
  ]
end 

to go
  repeat 500 [ask one-of patches [update1] ]                                                   ; calls each patches to go through the update1 function in which will cause the patches to start flipping
  tick                                                                                         ; allows the model to coninuously run until the "stopping" condition is met, as seen bellow
  ask patches [recolor]                                                                        ; recolors the patches based on their orientation

  if world_division = false                                                                    ; spin groups have equal number of patches populated               ( 1/2 )
  [
      ask patches
      [
          ifelse pxcor > (world-width / 2 - 1)
              [ set group 2] [ set group 1]                                                    ; sets the spin group for each patch as 1(blue) or 2(pink) based on the patches location in the world
      ]
  ]
  if world_division = true                                                                     ; spin group 1(blue) has more patches than spin group 2(pink)    ( 3/4 vs 1/4 )
  [
      ask patches
      [
          ifelse (pxcor > (world-width / 2 - 1)) and (pycor <= (world-height / 2 - 1))
              [ set group 2] [ set group 1]                                                    ; sets the spin group for each patch as 1(blue) or 2(pink) based on the patches location in the world
      ]
  ]


  ; STOPPING CONDITION
  ; When thermal equalibrium is reached then the model will stop, this occurs when the temperatures of each spin groupd are equal to eachother

  let temp_diff abs(temperature_a_update - temperature_b_update)                              ; tests the differance between the two temperatures as an absolute value
  let tempstop 10                                                                             ; local dynamic variable used when the model has run for an excess amount of time,
                                                                                              ; therefore "tempstop" is the range is which the two temps can be within eachother that will cause the model to stop
  ifelse ticks < 5000
  [   if temp_diff <= 1 [stop]         ]                                                      ; before 5000 ticks the two temps must come between one degree of eachother
  [   if temp_diff <= tempstop [stop]  ]                                                      ; after 5000 ticks, the two temps must come between "tempstop" degrees of one another
  if ticks > 50000 [stop]                                                                     ; if the model has run on for more than 50,000 ticks then it will also stop here



  set N_total_a (n_zero_a + n_one_a + n_two_a)
  set N_total_b (n_zero_b + n_one_b + n_two_b)
  ask patches[ hold ]                                                                         ; calls each patch to run through the temperature holding condition, used to troubleshoot runtime errors
  update-plots                                                                                ; updates all plots in real time as the model is running so that the observer can register what the simulation is doing
end 


; FLIPPING ORIENTATION
; A specific patch can be infulenced based it surrounding, neighboring patches
; Beacause of this the energy levels (described here as "orientation") can flip from one patch to another
; In some cases there are multiple probabilities for flipping, when this occurs a random generator is used not to cause bias
; There are 14 possibilites once it is determinded that the ttal energy level of the selected patch is different than the total energy of the neighboring patches
; Broken down into six cases with subcases under each

to update1
  let sumNeighbors ( sum[orientation] of neighbors )                                          ; adds up the energy levels of all surrounding patches of a specific patch
  set A orientation                                                                           ; this variable becomes the orientation of the selected patch, this is used as a "holding" variable so that the new orientation is not overriding

  if (sumNeighbors > 0 and A <= 0) or                                                         ; checks if the neighboring sum is different then that of the selected patch
     (sumNeighbors = 0 and A != 0) or                                                         ; if in fact different, then it will flip
     (sumNeighbors < 0 and A >= 0)

  [
      ask n-of N_total patches [ set friend? false ]                                          ; choses a random neighbor and will flip based on that neighbors orientation
      set neighbor1 one-of neighbors                                                          ; there are times when two neighboring flips are possible when jumping from the lowest to the highest energy state or highest to lowest
      ask neighbor1 [set friend? true ]
      set neighbor2 one-of neighbors with [friend? = false]
      set B1 [orientation] of neighbor1                                                       ; orientation of the first neighbor
      set B2 [orientation] of neighbor2                                                       ; orientation of the second neighbor

      let num random-float 1                                                                  ; random generator

    ;  CASE #1
    if A = 1 and B1 = -1
    [
      ifelse num > .5
      [
        set orientation 0                               ; case 1a
        ask neighbor1 [ set orientation 0]
      ]
      [
        set orientation -1                              ; case 1b
        ask neighbor1 [set orientation 1]
      ]
    ]
    ;  CASE #2
    ifelse A = 1 and B1 = 0 and B2 = 0
    [
      set orientation -1                                 ; case 2a
      ask neighbor1 [set orientation 1]
      ask neighbor2 [set orientation 1]
    ]
    [
      if A = 1 and B1 = 0
      [
      ifelse num > .5
      [
          ifelse B2 = 0
          [
            set orientation -1                           ; case 2b
            ask neighbor1 [set orientation 1]
            ask neighbor2 [set orientation 1]
          ]
          [
            set orientation 0                            ; case 2c
            ask neighbor1 [set orientation 1]
          ]
      ]
      [
        set orientation 0                                ; case 2d
        ask neighbor1 [set orientation 1]
      ]
      ]
    ]
    ; CASE #3
    if A = 0 and B1 = 1
    [
      set orientation 1                                  ; case 3a
      ask neighbor1 [set orientation 0]
    ]
    ; CASE #4
    if A = 0 and B1 = 0
    [
      ifelse num > .5
      [
        set orientation 1                                ; case 4a
        ask neighbor1 [set orientation -1]
      ]
      [
        set orientation -1                               ; case 4b
        ask neighbor1 [set orientation 1]
      ]
    ]
    ; CASE #5
      if A = -1 and B1 = 1
      [
              ifelse num > .5
              [
                  set orientation 0                      ; case 5a
                  ask neighbor1 [set orientation 0]
              ]
              [
                  set orientation 1                      ; case 5b
                  ask neighbor1 [set orientation -1]
              ]
      ]

    ; CASE # 6
    if A = -1 and B1 = 0
    [
      ifelse num > .5 and B2 = 0
      [
        set orientation 1                                ; case 6a
        ask neighbor1 [set orientation -1]
        ask neighbor2 [set orientation -1]
      ]
      [
        set orientation 0                                ; case 6b
        ask neighbor1 [set orientation -1]
      ]
    ]
    ; CASE #7
    if A = 0 and B1 = -1
    [
      set orientation -1                                 ; case 7a
      ask neighbor1 [set orientation 0]
    ]
  ]

  set n_zero_a_update count patches with [group = 1 and orientation = 0 ]    ; resets a new updated variable for the ever changing number of patches in each state based on orientation set in code above
  set n_one_a_update  count patches with [group = 1 and orientation = -1]
  set n_two_a_update  count patches with [group = 1 and orientation = 1]

  set n_zero_b_update count patches with [group = 2 and orientation = 0 ]    ; resets a new updated variable for the ever changing number of patches in each state based on orientation set in code above
  set n_one_b_update  count patches with [group = 2 and orientation = -1]
  set n_two_b_update  count patches with [group = 2 and orientation = 1]


  set N_total_a_update (n_zero_a_update + n_one_a_update + n_two_a_update)  ; verifies that the total count of patches is correct
  set N_total_b_update (n_zero_b_update + n_one_b_update + n_two_b_update)
end 

to recolor                                                                  ; coloring is dependant on each patches orientation
  ifelse group = 1
  [
    if orientation = 1  [set pcolor blue - 2]                               ; darkest blue = highest energy state of  ( + epsilon )
    if orientation = -1 [set pcolor blue + 2]                               ; lightest blue = lowest energy state of  ( - epsilon )
    if orientation = 0  [set pcolor blue]                                   ; middle blue = middle energy state of    ( 0 )
  ]
  [
    if orientation = 1  [set pcolor pink - 2]                               ; darkest pink = highest energy state of  ( + epsilon )
    if orientation = -1 [set pcolor pink + 2]                               ; lightest pink = lowest energy state of  ( - epsilon )
    if orientation = 0  [set pcolor pink]                                   ; middle pink = middle energy state of    ( 0 )
  ]
end 


;;;;; Thermaldynamic Equations ;;;;;;;



; SPIN GROUP A EQUATIONS BLUE

to-report entropy_a         ; derived using the sterling approximation and from the multiplicity
    report (N_total_a * ln(N_total_a) - N_total_a) - (n_two_a * ln(n_two_a) - n_two_a) - (n_one_a * ln(n_one_a) - n_one_a) - (n_zero_a * ln(n_zero_a) - n_zero_a)
end 

to-report entropy_a_update  ; new entropy as time goes on will always be increasing as seen in graphs
report   (N_total_a * ln(N_total_a) - N_total_a) - (n_two_a_update * ln(n_two_a_update) - n_two_a_update) - (n_one_a_update * ln(n_one_a_update) - n_one_a_update) - (n_zero_a_update * ln(n_zero_a_update) - n_zero_a_update)
end 

to-report U_a              ; total energy is always constant
  report (n_two_a - n_one_a) * epsilon
end 

to-report U_a_update       ; total energy is always constant
  report (n_two_a_update - n_one_a_update) * epsilon
end 

to-report temperature_a ; derived using mathematica from entropy
  ifelse (n_one_a = n_two_a) or (n_one_a = n_zero_a) or (n_two_a = n_zero_a) [report 0]     ; if 2 states have equal number of patches populated then 0 will be reported as to not receive ive an error
  [
     report
  ( epsilon /
    ln(
      ( n_one_a - n_two_a + ( (n_two_a - n_one_a) ^ 2 - 4 * (-2 * n_one_a - n_zero_a) * (2 * n_two_a + n_zero_a) ) ^ .5)
                        /( 2 * (2 * n_two_a + n_zero_a) )
      )
  )
  ]
end 

to-report temperature_a_update ; derived using mathematica with new populations
  ifelse (n_one_a_update = n_two_a_update) or (n_one_a_update = n_zero_a_update) or (n_two_a_update = n_zero_a_update) [report tempa_hold]          ; if 2 states have equal number of patches populated then the previous run's temp will be reported
  [
        report
  ( epsilon /
    ln(
      ( n_one_a_update - n_two_a_update + ( (n_two_a_update - n_one_a_update) ^ 2 - 4 * (-2 * n_one_a_update - n_zero_a_update) * (2 * n_two_a_update + n_zero_a_update) ) ^ .5)
                        /( 2 * (2 * n_two_a_update + n_zero_a_update) )
      )
  )
  ]
end 

to-report free_energy_a  ; formula from from 2 state system
  report ( U_a - (temperature_a * entropy_a) )
end 

to-report free_energy_a_update  ; from 2 state   ; will change over time because entropy is changing
  report ( U_a_update - (temperature_a_update * entropy_a_update) )
end 
;to-report specific_heat_a  ; derived using mathematica from total energy ( U ) and temperature
;  report (
;
;
;            (  2 * (  ((epsilon * e ^ (-1 * epsilon / temperature_a)) / (temperature_a ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_a)) / (temperature_a ^ 2)))   * temperature_a * N_total  )
;                      / (  1 + e ^ (-1 * epsilon / temperature_a) + e ^ (epsilon / temperature_a)  )
;
;    +  (
;            ( (    (((epsilon ^ 2) * e ^ (-1 * epsilon / temperature_a)) / (temperature_a ^ 4))
;                 + (((epsilon ^ 2) * e ^ (epsilon / temperature_a)) / (temperature_a ^ 4))
;                 - ((epsilon       * e ^ (-1 * epsilon / temperature_a)) / (temperature_a ^ 3))
;                 + ((epsilon       * e ^ (epsilon / temperature_a)) / (temperature_a ^ 3))
;             ) * (temperature_a ^ 2) * (N_total)  )
;         / (  1 + e ^ (-1 * epsilon / temperature_a) + e ^ (epsilon / temperature_a)    )
;       )
;   -      (
;            ((  ((epsilon * e ^ (-1 * epsilon / temperature_a)) / (temperature_a ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_a)) / (temperature_a ^ 2))) )
;                     / ( (  1 + e ^ (-1 * epsilon / temperature_a) + e ^ (epsilon / temperature_a)    ) ^ 2 )
;
;          )
;
;         )
;end
;to-report specific_heat_a_update  ; derived using mathematica from total energy ( U ) and temperature
;
;  report (
;
;
;            (  2 * (  ((epsilon * e ^ (-1 * epsilon / temperature_a_update)) / (temperature_a_update ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_a_update)) / (temperature_a_update ^ 2)))   * temperature_a_update * N_total  )
;                      / (  1 + e ^ (-1 * epsilon / temperature_a_update) + e ^ (epsilon / temperature_a_update)  )
;
;    +  (
;            ( (    (((epsilon ^ 2) * e ^ (-1 * epsilon / temperature_a_update)) / (temperature_a_update ^ 4))
;                 + (((epsilon ^ 2) * e ^ (epsilon / temperature_a_update)) / (temperature_a_update ^ 4))
;                 - ((epsilon       * e ^ (-1 * epsilon / temperature_a_update)) / (temperature_a_update ^ 3))
;                 + ((epsilon       * e ^ (epsilon / temperature_a_update)) / (temperature_a_update ^ 3))
;             ) * (temperature_a_update ^ 2) * (N_total)  )
;         / (  1 + e ^ (-1 * epsilon / temperature_a_update) + e ^ (epsilon / temperature_a_update)    )
;       )
;   -      (
;            ((  ((epsilon * e ^ (-1 * epsilon / temperature_a_update)) / (temperature_a_update ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_a_update)) / (temperature_a_update ^ 2))) )
;                     / ( (  1 + e ^ (-1 * epsilon / temperature_a_update) + e ^ (epsilon / temperature_a_update)    ) ^ 2 )
;
;          )
;
;         )
;end


; GROUP B EQUATIONS PINK

to-report entropy_b         ; derived using the sterling approximation and from the multiplicity
    report (N_total_b * ln(N_total_b) - N_total_b) - (n_two_b * ln(n_two_b) - n_two_b) - (n_one_b * ln(n_one_b) - n_one_b) - (n_zero_b * ln(n_zero_b) - n_zero_b)
end 

to-report entropy_b_update  ; new entropy_b as time goes on will always be increasing as seen in graphs
report   (N_total_b * ln(N_total_b) - N_total_b) - (n_two_b_update * ln(n_two_b_update) - n_two_b_update) - (n_one_b_update * ln(n_one_b_update) - n_one_b_update) - (n_zero_b_update * ln(n_zero_b_update) - n_zero_b_update)
end 

to-report temperature_b ; derived using mathematica from entropy_b
  ifelse (n_one_b = n_two_b) or (n_one_b = n_zero_b) or (n_two_b = n_zero_b) [report 0]
  [
    report
  ( epsilon /
    ln(
      ( n_one_b - n_two_b + ( (n_two_b - n_one_b) ^ 2 - 4 * (-2 * n_one_b - n_zero_b) * (2 * n_two_b + n_zero_b) ) ^ .5)
                        /( 2 * (2 * n_two_b + n_zero_b) )
      )
  )
  ]
end 

to-report temperature_b_update ; derived using mathematica with new populations
  ifelse (n_one_b_update = n_two_b_update) or (n_one_b_update = n_zero_b_update) or (n_two_b_update = n_zero_b_update) [report tempb_hold]
  [
        report
  ( epsilon /
    ln(
      ( n_one_b_update - n_two_b_update + ( (n_two_b_update - n_one_b_update) ^ 2 - 4 * (-2 * n_one_b_update - n_zero_b_update) * (2 * n_two_b_update + n_zero_b_update) ) ^ .5)
                        /( 2 * (2 * n_two_b_update + n_zero_b_update) )
      )
  )
  ]
end 

to-report U_b               ; total energy is always constant
  report ( n_two_b - n_one_b ) * epsilon
end 

to-report U_b_update        ; verifies that the total energy is stayig constant based on the new populations
  report ( n_two_b_update - n_one_b_update ) * epsilon
end 

to-report free_energy_b  ; formula from from 2 state system
  report ( U_b - (temperature_b * entropy_b) )
end 

to-report free_energy_b_update  ; from 2 state   ; will change over time because entropy is changing
  report ( U_b_update - (temperature_b_update * entropy_b_update) )
end 
;to-report specific_heat_b  ; derived using mathematica from total energy ( U ) and temperature
;  report (
;
;
;            (  2 * (  ((epsilon * e ^ (-1 * epsilon / temperature_b)) / (temperature_b ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_b)) / (temperature_b ^ 2)))   * temperature_b * N_total_b  )
;                      / (  1 + e ^ (-1 * epsilon / temperature_b) + e ^ (epsilon / temperature_b)  )
;
;    +  (
;            ( (    (((epsilon ^ 2) * e ^ (-1 * epsilon / temperature_b)) / (temperature_b ^ 4))
;                 + (((epsilon ^ 2) * e ^ (epsilon / temperature_b)) / (temperature_b ^ 4))
;                 - ((epsilon       * e ^ (-1 * epsilon / temperature_b)) / (temperature_b ^ 3))
;                 + ((epsilon       * e ^ (epsilon / temperature_b)) / (temperature_b ^ 3))
;             ) * (temperature_b ^ 2) * (N_total_b)  )
;         / (  1 + e ^ (-1 * epsilon / temperature_b) + e ^ (epsilon / temperature_b)    )
;       )
;   -      (
;            ((  ((epsilon * e ^ (-1 * epsilon / temperature_b)) / (temperature_b ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_b)) / (temperature_b ^ 2))) )
;                     / ( (  1 + e ^ (-1 * epsilon / temperature_b) + e ^ (epsilon / temperature_b)    ) ^ 2 )
;
;          )
;
;         )
;end
;to-report specific_heat_b_update  ; derived using mathematica from total energy ( U ) and temperature
;  report (
;
;
;            (  2 * (  ((epsilon * e ^ (-1 * epsilon / temperature_b_update)) / (temperature_b_update ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_b_update)) / (temperature_b_update ^ 2)))   * temperature_b_update * N_total_b  )
;                      / (  1 + e ^ (-1 * epsilon / temperature_b_update) + e ^ (epsilon / temperature_b_update)  )
;
;    +  (
;            ( (    (((epsilon ^ 2) * e ^ (-1 * epsilon / temperature_b_update)) / (temperature_b_update ^ 4))
;                 + (((epsilon ^ 2) * e ^ (epsilon / temperature_b_update)) / (temperature_b_update ^ 4))
;                 - ((epsilon       * e ^ (-1 * epsilon / temperature_b_update)) / (temperature_b_update ^ 3))
;                 + ((epsilon       * e ^ (epsilon / temperature_b_update)) / (temperature_b_update ^ 3))
;             ) * (temperature_b_update ^ 2) * (N_total_b)  )
;         / (  1 + e ^ (-1 * epsilon / temperature_b_update) + e ^ (epsilon / temperature_b_update)    )
;       )
;   -      (
;            ((  ((epsilon * e ^ (-1 * epsilon / temperature_b_update)) / (temperature_b_update ^ 2))   -   ((epsilon * e ^ (epsilon / temperature_b_update)) / (temperature_b_update ^ 2))) )
;                     / ( (  1 + e ^ (-1 * epsilon / temperature_b_update) + e ^ (epsilon / temperature_b_update)    ) ^ 2 )
;
;          )
;
;         )
;
;end

to hold                                                 ; in case of error when values of n_two n_one or n_zero equal eachother
  set tempa_hold temperature_a_update                   ; sets temphold as the previous run's temp value
  set tempb_hold temperature_b_update
end 

to-report total_energy
  report U_a_update + U_b_update
end 

to-report final_total_entropy
  report entropy_a_update + entropy_b_update
end 

to-report initial_total_entropy
  report entropy_a + entropy_b
end