;; System dynamics model globals
globals [
  ;; constants
  penguin-death-rate
  ;; stock values
  BreedingGround
  Krill
  Penguin
  TotalTemperatureChange
  ;; size of each step, see SYSTEM-DYNAMICS-GO
  dt
]

;; Initializes the system dynamics model.
;; Call this in your model's SETUP procedure.
to system-dynamics-setup
  reset-ticks
  set dt 1.0
  ;; initialize constant values
  set penguin-death-rate 0.08748
  ;; initialize stock values
  set BreedingGround 57000
  set Krill 50000
  set Penguin 1000
  set TotalTemperatureChange 0
end

;; Step through the system dynamics model by performing next iteration of Euler's method.
;; Call this in your model's GO procedure.
to system-dynamics-go

  ;; compute variable and flow values once per step
  let local-penguin-birth-rate penguin-birth-rate
  let local-krill-birth-rate krill-birth-rate
  let local-krill-death-rate krill-death-rate
  let local-penguin-deaths penguin-deaths
  let local-penguin-births penguin-births
  let local-krill-deaths krill-deaths
  let local-krill-births krill-births
  let local-temperature-change-per-month temperature-change-per-month
  let local-BreedingGroundLoss BreedingGroundLoss

  ;; update stock values
  ;; use temporary variables so order of computation doesn't affect result.
  let new-BreedingGround max( list 0 ( BreedingGround - local-BreedingGroundLoss ) )
  let new-Krill max( list 0 ( Krill - local-krill-deaths + local-krill-births ) )
  let new-Penguin max( list 0 ( Penguin - local-penguin-deaths + local-penguin-births ) )
  let new-TotalTemperatureChange max( list 0 ( TotalTemperatureChange + local-temperature-change-per-month ) )
  set BreedingGround new-BreedingGround
  set Krill new-Krill
  set Penguin new-Penguin
  set TotalTemperatureChange new-TotalTemperatureChange

  tick-advance dt
end

;; Report value of flow
to-report penguin-deaths
  report ( Penguin * penguin-death-rate
  ) * dt
end

;; Report value of flow
to-report penguin-births
  report ( (Penguin * penguin-birth-rate * Krill) / 50000
  ) * dt
end

;; Report value of flow
to-report krill-deaths
  report ( (Penguin * Krill * krill-death-rate) / 1000
  ) * dt
end

;; Report value of flow
to-report krill-births
  report ( Krill * krill-birth-rate
  ) * dt
end

;; Report value of flow
to-report temperature-change-per-month
  report ( user-Annual-Temperature-Increase / 12
  ) * dt
end

;; Report value of flow
to-report BreedingGroundLoss
  report ( e ^ ((TotalTemperatureChange - 0.9495) / 0.2627) - e ^ ((TotalTemperatureChange - (user-Annual-Temperature-Increase / 12) - 0.9495) / 0.2627)
  ) * dt
end

;; Report value of variable
to-report penguin-birth-rate
  report 1 / 12 * BreedingGround / 57000
end

;; Report value of variable
to-report krill-birth-rate
  report 0.987 / 12 * BreedingGround / 57000
end

;; Report value of variable
to-report krill-death-rate
  report 0.987 / 12
end

;; Plot the current state of the system dynamics model's stocks
;; Call this procedure in your plot's update commands.
to system-dynamics-do-plot
  if plot-pen-exists? "BreedingGround" [
    set-current-plot-pen "BreedingGround"
    plotxy ticks BreedingGround
  ]
  if plot-pen-exists? "Krill" [
    set-current-plot-pen "Krill"
    plotxy ticks Krill
  ]
  if plot-pen-exists? "Penguin" [
    set-current-plot-pen "Penguin"
    plotxy ticks Penguin
  ]
  if plot-pen-exists? "TotalTemperatureChange" [
    set-current-plot-pen "TotalTemperatureChange"
    plotxy ticks TotalTemperatureChange
  ]
end
