Conduction in solids

globals [
  tick-advance-amount               ; how much we advance the tick counter this time through
  max-tick-advance-amount           ; the largest tick-advance-amount is allowed to be
  init-avg-speed init-avg-energy    ; initial averages
  avg-speed avg-energy              ; current average
  avg-temp-type-1
   avg-temp-type-2
   avg-temp-type-3

  avg-energy-gray
  particle-size
  toggle-red-state
  toggle-green-state
  min-particle-energy
  max-particle-energy
  gas-particles-to-add
  initial-gas-temperature

  #-3-wall-particles
  #-2-wall-particles
  #-1-wall-particles

  #-3-wall-total-te
  #-2-wall-total-te
  #-1-wall-total-te

  init-temp#1
  init-temp#2
  init-temp#3

  left-piece-label
  middle-piece-label
  right-piece-label
  mouse-interaction
  kelvin-celcius-offset
  a
  b
]

breed [ gas-particles particle ]
breed [ walls wall ]
breed [ flashes flash ]
breed [ erasers eraser ]
breed [ arrowheads arrowhead ]
breed [ boundaries boundary ]

erasers-own [ pressure? ]
flashes-own [ birthday ]

gas-particles-own [
  speed mass energy          ; gas-particles info
  last-collision
  color-type
]

walls-own [
  energy
  valve-1?
  valve-2?
  pressure?
  my-pxcor
  my-pycor
  my-heading
  distance-advance
  color-type
  particle-type
  init-yoffset
]

to setup
  clear-all
  reset-ticks
  set particle-size 1.0
  set max-tick-advance-amount 0.02

  set #-3-wall-particles 0
  set #-2-wall-particles 0
  set #-1-wall-particles 0
  set kelvin-celcius-offset 273
  set gas-particles-to-add 2
  set mouse-interaction "none"
  set a 200
  set b .75
 ; set initial-gas-temperature 500

  set init-temp#1 left.start.temp ;+ kelvin-celcius-offset
  set init-temp#2 middle.start.temp ;+ kelvin-celcius-offset
  set init-temp#3 right.start.temp; + kelvin-celcius-offset

  set-default-shape flashes "square2"
  set-default-shape walls "circle"
  set-default-shape erasers "eraser"
  set-default-shape arrowheads "default"

  set min-particle-energy 100
  set max-particle-energy 500  ;(.5 ) * ( max-dist-in-tick-advance-amount  / max-tick-advance-amount ) ^ 2

  create-erasers 1 [ set hidden? true set size 1.2 set color [255 255 255 100] ]

  make-box
  ask walls [transfer-energy-mytype]
  ask walls [recolorshade-walls]

  ask gas-particles [ apply-speed-visualization ]

  set init-avg-speed avg-speed
  set init-avg-energy avg-energy

  update-variables
  do-plots
end 

to go
   if end-simulation? [ stop ]
  ask walls [transfer-energy]
  ask walls [calculate-distance-advance-amount]

  ask walls [apply-vibration-motion-visualization]
  ask gas-particles [ apply-speed-visualization ]
  ask walls [recolorshade-walls]

  update-variables
  do-plots

  calculate-tick-advance-amount
  tick-advance tick-advance-amount

  display
end 

to calculate-distance-advance-amount
   set distance-advance (distance-advance + (tick-advance-amount * (energy) * 9 ))
end 

to-report end-simulation?
  let value-to-report false
  if auto-stop-simulation? = "at time = 5" and ticks >= 5    [set value-to-report true]
  if auto-stop-simulation? = "at time = 10" and ticks >= 10  [set value-to-report true]
  if auto-stop-simulation? = "at time = 15" and ticks >= 15  [set value-to-report true]
  if auto-stop-simulation? = "at time = 20" and ticks >= 20  [set value-to-report true]
 ; set value-to-report true
  report value-to-report
end 

to update-variables

  if any? gas-particles [
    set avg-speed  mean [ speed ] of gas-particles
    set avg-energy( mean [ energy ] of gas-particles - kelvin-celcius-offset)
  ]
  let wall-particles-type-1 walls with [particle-type = 1]
  let wall-particles-type-2 walls with [particle-type = 2]
  let wall-particles-type-3 walls with [particle-type = 3]



  ;show count gas-particles-type-1
  ifelse any? wall-particles-type-1 [
     set #-1-wall-particles count wall-particles-type-1
     set avg-temp-type-1 ( (mean [ energy ] of wall-particles-type-1)  - kelvin-celcius-offset)
     set #-1-wall-total-te ((mean [ energy ] of wall-particles-type-1)  * #-1-wall-particles )
     set left-piece-label avg-temp-type-1
   ; show "yes"
  ]
  [set left-piece-label "N/A"]
  ifelse any? wall-particles-type-2 [
     set #-2-wall-particles count wall-particles-type-2
     set avg-temp-type-2 ( (mean [ energy ] of wall-particles-type-2)  - kelvin-celcius-offset)
     set #-2-wall-total-te ((mean [ energy ] of wall-particles-type-2)  * #-2-wall-particles )
     set middle-piece-label avg-temp-type-2
   ; show "yes"
  ]
  [set right-piece-label "N/A"]
    ifelse any? wall-particles-type-3 [
     set #-3-wall-particles count wall-particles-type-3
     set avg-temp-type-3 ( (mean [ energy ] of wall-particles-type-3)  - kelvin-celcius-offset)
     set #-3-wall-total-te ((mean [ energy ] of wall-particles-type-3) * #-3-wall-particles )
     set right-piece-label avg-temp-type-3
   ; show "yes"
  ]
  [set right-piece-label "N/A"]
end 



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;GAS MOLECULES MOVEMENT;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to bounce  ; gas-particles procedure
  ; get the coordinates of the patch we'll be on if we go forward 1
  let bounce-patch nobody
  let bounce-patches nobody
  let hit-angle 0
  let this-patch patch-here
  let new-px 0
  let new-py 0
  let visible-wall nobody
  let particle-energy energy

  set bounce-patch  min-one-of walls in-cone ((sqrt (2)) / 2) 180 with [ myself != this-patch ] [ distance myself ]

  if bounce-patch != nobody [
    let wall-energy [energy] of bounce-patch
    let energy-difference (particle-energy - wall-energy)
    set new-px [ pxcor ] of bounce-patch
    set new-py [ pycor ] of bounce-patch
    set visible-wall walls-on bounce-patch

    if any? visible-wall with [ not hidden? ]  [
      set hit-angle towards bounce-patch
      ifelse (hit-angle <= 135 and hit-angle >= 45) or (hit-angle <= 315 and hit-angle >= 225) [
        set heading (- heading)
      ][
        set heading (180 - heading)
      ]
      set energy (energy - (energy-difference / 4))
      set speed speed-from-energy
      ask bounce-patch [set energy (energy + (energy-difference / 4)) ]
      if true [; show-wall-hits? [
        ask patch new-px new-py [
          sprout 1 [
            set breed flashes
            set color gray - 2
            set birthday ticks
          ]
        ]
      ]
    ]
  ]
end 

to rewind-to-bounce  ; gas-particles procedure
  ; attempts to deal with particle penetration by rewinding the particle path back to a point
  ; where it is about to hit a wall
  ; the particle path is reversed 49% of the previous tick-advance-amount it made,
  ; then particle collision with the wall is detected again.
  ; and the particle bounces off the wall using the remaining 51% of the tick-advance-amount.
  ; this use of slightly more of the tick-advance-amount for forward motion off the wall, helps
  ; insure the particle doesn't get stuck inside the wall on the bounce.

  let bounce-patch nobody
  let bounce-patches nobody
  let hit-angle 0
  let this-patch nobody
  let new-px 0
  let new-py 0
  let visible-wall nobody

  bk (speed) * tick-advance-amount * .49
  set this-patch  patch-here

  set bounce-patch  min-one-of walls in-cone ((sqrt (2)) / 2) 180 with [ self != this-patch ] [ distance myself ]

  if bounce-patch != nobody [

    set new-px [pxcor] of bounce-patch
    set new-py [pycor] of bounce-patch
    set visible-wall walls-on bounce-patch

    if any? visible-wall with [not hidden?] [
      set hit-angle towards bounce-patch

      ifelse (hit-angle <= 135 and hit-angle >= 45) or (hit-angle <= 315 and hit-angle >= 225) [
        set heading (- heading)
      ][
        set heading (180 - heading)
      ]

      if true [;[show-wall-hits? [
        ask patch new-px new-py [
          sprout 1 [
            set breed flashes
            set color gray - 2
            set birthday ticks
          ]
        ]
      ]
    ]
  ]
  fd (speed) * tick-advance-amount * .51
end 

to move  ; gas-particles procedure
  if patch-ahead (speed * tick-advance-amount) != patch-here [ set last-collision nobody ]
  jump (speed * tick-advance-amount)
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;GAS MOLECULES COLLISIONS;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;from GasLab

to calculate-tick-advance-amount
  set tick-advance-amount max-tick-advance-amount
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  mouse interaction procedures
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to mouse-action
 let snap-xcor 0
 let snap-ycor 0
 let orig-xcor 0
 let orig-ycor 0
 let eraser-window-walls nobody
 let eraser-window-gas-particles nobody

  ifelse mouse-down? [
    set orig-xcor mouse-xcor
    set orig-ycor mouse-ycor
    set snap-xcor round orig-xcor
    set snap-ycor round orig-ycor

    ask patches with [ pxcor = snap-xcor and pycor = snap-ycor ] [
      set eraser-window-walls  walls-here
     ; show eraser-window-walls
      set eraser-window-gas-particles gas-particles-here



      if mouse-interaction = "speed up gas-particles"  [

        ask erasers [
          set hidden? false
          set shape "spray paint"
          setxy orig-xcor orig-ycor
        ]
        ask eraser-window-gas-particles [
          set energy (energy * 1.1)
          set energy limited-particle-energy
          set speed speed-from-energy
        ;  apply-vibration-color-visualization
        ]
        ask eraser-window-walls [
          set energy (energy * 1.1)
          set energy limited-particle-energy
         ; apply-vibration-color-visualization
        ]



      ]

      if mouse-interaction = "slow down gas-particles" [
        ask erasers [
          set hidden? false
          set shape "spray paint"
          setxy orig-xcor orig-ycor
        ]
        ask eraser-window-gas-particles [
          set energy (energy / 1.1)
          set energy limited-particle-energy
          set speed speed-from-energy
          apply-speed-visualization
        ]
        ask eraser-window-walls [
          set energy (energy / 1.1)
          set energy limited-particle-energy
        ;  apply-vibration-color-visualization
        ]
      ]
    ]
    ask gas-particles with [ any? walls-here ] [ remove-from-walls ] ; deal with any walls drawn on top of gas-particles
  ][
    ask erasers [ set hidden? true ]
  ]
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; thermal energy transfer
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to transfer-energy
  let neighboring-walls walls-on neighbors4
  let my-energy energy
  let total-energy-from-neighbors 0

  if any? neighboring-walls [

  ask neighboring-walls [
    let energy-from-me (((energy) * (0.25)))
    set energy (energy - energy-from-me)
    set total-energy-from-neighbors total-energy-from-neighbors + energy-from-me
  ]

  set energy energy + total-energy-from-neighbors
  ]
end 

to transfer-energy-mytype
  let my-type particle-type
  let neighboring-walls-all walls-on neighbors4
  let neighboring-walls neighboring-walls-all with [particle-type = my-type]
  let my-energy energy
  let total-energy-from-neighbors 0

  if any? neighboring-walls [

  ask neighboring-walls [
    let energy-from-me (((energy) * (0.25)))
    set energy (energy - energy-from-me)
    set total-energy-from-neighbors total-energy-from-neighbors + energy-from-me
  ]

  set energy energy + total-energy-from-neighbors
  ]
end 

to check-slide-up
  if particle-type = 1 [set ycor slide-left-solid-up + init-yoffset]
  if particle-type = 2 [set ycor slide-middle-solid-up + init-yoffset]
  if particle-type = 3 [set ycor slide-right-solid-up + init-yoffset]
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; particle speed and flash visualization procedures
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to apply-vibration-motion-visualization
  set xcor my-pxcor
  check-slide-up
 ; set my-heading my-heading + (energy / 100) * (random-float 5  - random-float 5) ; (energy / 100)  - random-float (energy / 100)
;  set heading my-heading
  let displacement (0.05 * ((energy / a) ^ b) *  sin (distance-advance) )
  fd displacement
end 

to apply-speed-visualization
  recolorshade-gas-particles
 ; if visualize-speed? = "arrows" [ scale-arrowheads ]
 ; if visualize-speed? = "different shades" [ recolorshade-gas-particles ]
 ; if visualize-speed? = "none" [ recolornone ]
end 

to color-particle-and-link
  let this-link my-out-links
  let this-color-type color-type
  set color this-color-type
  ask this-link [ set color this-color-type ]
end 

to scale-arrowheads
  let this-xcor xcor
  let this-ycor ycor
  let this-speed speed
  let this-heading heading
  let this-arrowhead out-link-neighbors
  let this-link my-out-links
  ask this-link [ set hidden? false ]
  ask this-arrowhead [
    set xcor this-xcor
    set ycor this-ycor
    set heading this-heading
    fd .5 + this-speed / 3
  ]
end 

to-report my-color-type
  let this-color-type color-type
  report this-color-type
end 

to recolorshade-gas-particles
  let this-link my-out-links
  let this-arrowhead out-link-neighbors

  let local-color-type my-color-type
  let energy-difference (energy)
  let r-value ((4.9) * (1 - 0.999 ^ (energy)))


  let this-color ((my-color-type + 4.99) - r-value)

  set color this-color
  ask this-link [ set color this-color ]
  ask this-arrowhead [ set color this-color ]
end 

to recolornone
  let this-link my-out-links
  let this-arrowhead out-link-neighbors
  let this-color color-type
  set color color-type
  ask this-link [set color this-color]
  ask this-arrowhead [ set color this-color ]
end 

to recolorshade-walls
  let this-link my-out-links
  let this-arrowhead out-link-neighbors

  let local-color-type my-color-type
  let energy-difference (energy)
  let r-value ((8.9) * (1 - 0.999 ^ (energy)))

  let this-color ((my-color-type + 4.99) - r-value)

  set color this-color
  ask this-link [ set color this-color ]
  ask this-arrowhead [ set color this-color ]
end 



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  initialization procedures
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to make-box
  let vertical-line min-pycor
  let vertical-line-xcor 0.5
  let vertical-line-offset 8

  create-boundaries 1 [
    set xcor (vertical-line-xcor)  set shape "new-line" set color [255 255 255 120] set size 17
  ]
  create-boundaries 1 [
    set xcor (vertical-line-xcor + left-solid-width) set shape "new-line" set color [255 255 255 120] set size 17
  ]
 create-boundaries 1 [
    set xcor (vertical-line-xcor + left-solid-width + middle-solid-width) set shape "new-line" set color [255 255 255 120] set size 17
  ]

 create-boundaries 1 [
    set xcor (vertical-line-xcor + left-solid-width + middle-solid-width + right-solid-width) set shape "new-line" set color [255 255 255 120] set size 17
  ]


  ask patches with [ ((pycor >= (min-pycor + slide-left-solid-up)) and (pycor < (min-pycor + slide-left-solid-up + left-solid-length )) )
      and (pxcor > min-pxcor ) and (pxcor <= (left-solid-width)) ] [
      sprout 1 [set breed walls set color-type red set size 1.0 set particle-type 1 set shape "rounded-square" initialize-this-wall set init-yoffset( pycor - slide-left-solid-up)]
  ]

  ask patches with [ ((pycor >= (min-pycor + slide-middle-solid-up )) and (pycor < (min-pycor + slide-middle-solid-up + middle-solid-length )) )
      and (pxcor > (min-pxcor + left-solid-width )  and (pxcor <=  (min-pxcor + left-solid-width  + middle-solid-width))) ] [
      sprout 1 [set breed walls set color-type red set particle-type 2 set size 1.0 set shape "rounded-square" initialize-this-wall set init-yoffset( pycor - slide-middle-solid-up) ]
  ]

  ask patches with [ ((pycor >= (min-pycor + slide-right-solid-up )) and (pycor < (min-pycor + slide-right-solid-up + right-solid-length )) )
      and (pxcor > (min-pxcor + left-solid-width + middle-solid-width))  and (pxcor <= (min-pxcor  + middle-solid-width + left-solid-width + right-solid-width)) ] [
      sprout 1 [set breed walls set color-type red set size 1.0  set particle-type 3 set shape "rounded-square" initialize-this-wall set init-yoffset( pycor - slide-right-solid-up)]
  ]
end 

to initialize-this-wall
  set valve-1? false
  set valve-2? false
  set pressure? false
  if particle-type =  1 [set energy left.start.temp + kelvin-celcius-offset]
  if particle-type =  2 [set energy middle.start.temp  + kelvin-celcius-offset]
  if particle-type =  3 [set energy right.start.temp  + kelvin-celcius-offset]
  set my-pxcor pxcor
  set my-pycor pycor
  set my-heading random 360
  set distance-advance 0
end 

to setup-gas-particles  ; gas-particles procedure
  set shape "circle"
  ;set size particle-size
  set energy initial-gas-temperature
  set color-type 5
  set color color-type
  set mass (10)  ; atomic mass
  hatch 1 [
    set breed arrowheads
    set hidden? true
    create-link-from myself [ tie ]
  ]
  set speed speed-from-energy
  set last-collision nobody
end 


; Place gas-particles at random, but they must not be placed on top of wall atoms.
; This procedure takes into account the fact that wall molecules could have two possible arrangements,
; i.e. high-surface area ot low-surface area.

to random-position ;; gas-particles procedure
  let open-patches nobody
  let open-patch nobody
  set open-patches patches with [not any? turtles-here and pxcor != max-pxcor and pxcor != min-pxcor and pycor != min-pycor and pycor != max-pycor]
  set open-patch one-of open-patches

  ; Reuven added the following "if" so that we can get through setup without a runtime error.
  if open-patch = nobody [
    user-message "No open patches found.  Exiting."
    stop
  ]

  setxy ([ pxcor ] of open-patch) ([ pycor ] of open-patch)
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; wall penetration error handling procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; if gas-particles actually end up within the wall

to remove-from-walls
  let this-wall walls-here with [ not hidden? ]

  if count this-wall != 0 [
    let available-patches patches with [ not any? walls-here ]
    let closest-patch nobody
    if (any? available-patches) [
      set closest-patch min-one-of available-patches [ distance myself ]
      set heading towards closest-patch
      setxy ([ pxcor ] of closest-patch)  ([ pycor ] of closest-patch)
    ]
  ]
end 


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;GRAPHS;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;REPORTERS;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report speed-from-energy
  report sqrt (2 * energy / mass)
end 

to-report energy-from-speed
  report (mass * speed * speed / 2)
end 

to-report limited-particle-energy
  let limited-energy energy
  if limited-energy > max-particle-energy [ set limited-energy max-particle-energy ]
  if limited-energy < min-particle-energy [ set limited-energy min-particle-energy ]
  report limited-energy
end 

to-report vibration-displacement
end 

to do-plots
  set-current-plot "Temperature of solids"
  if #-1-wall-particles > 0 [
    set-current-plot-pen "left"
    plotxy ticks avg-temp-type-1
  ]
    if #-2-wall-particles > 0 [
    set-current-plot-pen "middle"
    plotxy ticks avg-temp-type-2
  ]
    if #-3-wall-particles > 0 [
    set-current-plot-pen "right"
    plotxy ticks avg-temp-type-3
  ]
end 



; Copyright 2006 Uri Wilensky.
; See Info tab for full copyright and license.