Hyporheic Exchange Model: Mouse drop

breed [ particles particle ]

particles-own [ residence-time in-bank? ]
globals [Boundary constant k2 height time-counter time-dx HydrCond2 hm2 Boundary2 xcutoff U2 depth2]

to-report surface-function [ x ]
  report height * sin ( k2 * x ) + Slope * x
end 

to setup
  clear-all
  set-default-shape particles "circle"
  set constant 270 / lambda


  set Boundary BedDepth
  set k2 ( 360 / (constant * lambda ))
  set height ( BedformHeight * constant ) / 2
  set time-counter 0
  set time-dx 100
  set HydrCond2 constant * HydrCond
  set U2 ChannelVelocity * constant
  set depth2 depth * constant
  ifelse (( height * 2 / depth2 ) > .34)
  [set hm2 .28 * ( ( U2 ) ^ 2 ) / ( 2 * 9.8 * constant * 100  ) * (( height * 2 / depth2 ) / .34 ) ^ ( 3 / 2 ) ]
  [set hm2 .28 * ( ( U2 ) ^ 2 ) / ( 2 * 9.8 * constant * 100  ) * (( height * 2 / depth2 ) / .34 ) ^ ( 3 / 8 ) ]

  set Boundary2 constant * Boundary
  set xcutoff constant * lambda * floor( 270 / (constant * lambda))

  ask patches
  [
    ifelse ( pycor < surface-function (pxcor))
    [ set pcolor brown + 2 ] ; sand color
    [ set pcolor blue - 2 ]  ; water color

    if ( pycor < -1 * boundary2 )
    [ set pcolor black ]

    if ( pxcor > xcutoff )
    [ set pcolor black ]


    if ( pxcor < 0 )
    [ set pcolor black ]
  ]

  reset-ticks
end 

to go
  ask particles with [ in-bank? ]
  [
    move
  ]

  ifelse (time-counter = 0 and time1 > 0)
    [ change-color 65 ]
    [
      ifelse (time-counter < time1 * 60 and time-counter + time-dx >= time1 * 60)
      [ change-color 75 ]
      [ifelse (time-counter < time2 * 60 and time-counter + time-dx >= time2 * 60)
        [ change-color 83 ]
        [ifelse (time-counter < time3 * 60 and time-counter + time-dx >= time3 * 60)
          [ change-color 125 ]
          [ifelse (time-counter < time4 * 60 and time-counter + time-dx >= time4 * 60)
            [ change-color 115 ]
            [ifelse (time-counter < time5 * 60 and time-counter + time-dx >= time5 * 60)
              [ change-color 105 ]
              [ if (time-counter < time6 * 60 and time-counter + time-dx >= time6 * 60)
                [ change-color 95 ]
              ]]]]]]

  set time-counter time-counter + time-dx
  tick
end 

to timestep

  ifelse (time-counter < time1 * 60 )

  [
    repeat ((time1 * 60 - time-counter) / time-dx) [go]
  ]

  [ifelse (time-counter < time2 * 60 )
    [
      repeat ((time2 * 60 - time-counter) / time-dx) [go]
    ]
    [ifelse (time-counter < time3 * 60 )
      [
        repeat ((time3 * 60 - time-counter) / time-dx) [go]
      ]
      [ifelse (time-counter < time4 * 60 )
        [
          repeat ((time4 * 60 - time-counter) / time-dx) [go]
        ]
        [ifelse (time-counter < time5 * 60 )
          [
            repeat ((time5 * 60 - time-counter) / time-dx) [go]
          ]
          [ if (time-counter < time6 * 60 )
            [
              repeat ((time6 * 60 - time-counter) / time-dx) [go]
            ]
          ]]]]]
end 

to change-color [newcolor]
  ask particles with [in-bank?] [
    set color newcolor
    stamp
  ]
end 

to move

  let ycor2 ycor - surface-function( xcor )

  let exp1 exp ( k2 * (2 * pi / 360 ) * ycor2 )
  let exp2 exp ( k2 * (2 * pi / 360 ) * (-2 * boundary2 - ycor2 ) )
  let exp3 exp ( k2 * (2 * pi / 360 ) * (-2 * boundary2) )

  let expsum exp1 + exp2

  let yvel  ((-1 * (cos ( xcor * k2 )) * ( ( exp1 - exp2) * (k2 * 2 * pi / 360 )) * hm2)) / (1 + exp3) * HydrCond2 / porosity ;
  let xvel  ((k2 * (2 * pi / 360 ) * ( sin ( xcor * k2 )) * ( expsum )  * hm2 ) ) / (1 + exp3) * HydrCond2 / porosity;


  set xcor xcor +  time-dx * (xvel - HydrCond2 * Slope / (porosity * (1 + Slope ^ 2))) ; / constant

  if (ycor + yvel > max-pycor)
  [ show "big problem -- particle lost above screen" die ]
  if (ycor + yvel < min-pycor)
  [
    ;show "Oh no, I'm lost!"
    die
  ]
  set ycor ycor + time-dx * (yvel - (HydrCond2 * Slope ^ 2) / (porosity * (1 + Slope ^ 2))); + HydrCond2 / porosity );/ constant

  set residence-time residence-time + time-dx

  if xcor >= xcutoff
  [ pen-up
    set ycor ycor - Slope * xcutoff
    set xcor xcor - xcutoff

    if show-paths?
      [pen-down]
  ]

  if xcor <= 0
  [ pen-up
    set ycor ycor + Slope * xcutoff
    set xcor xcor + xcutoff

    if show-paths?
      [pen-down]
  ]

  set in-bank? (ycor <= surface-function xcor)
  if not (in-bank?)
  [
    ; we want them to stamp them right on the function, for better visuals
    set ycor surface-function xcor
    ht
  ]
end 

to mouse-drop
  if mouse-down?
  [
    ask patch mouse-xcor mouse-ycor [
      sprout-particles 1
      [
        set xcor mouse-xcor
        set ycor mouse-ycor
        if any? other particles in-radius 0.24 ; don't let them drop particles too close to each other
          [ die ]
        if (ycor > surface-function xcor) ; don't let them drop particles above the surface
          [ die ]
        set in-bank? true
        set residence-time 0
        set size 2
        set color 65
        if (show-paths?) [ pd ]
      ]
    ]
  ]
end 


;; COMMENTED OUT 10/15/2021 to allow running in NetLogoWeb

;to export-residence-times
;  let file user-new-file
;  if ( file != false )
;  [
;    if ( file-exists? file )
;    [ file-delete file ]
;    file-open file
;  ]
;  let res-time-list [ residence-time ] of particles with [ not in-bank? ]
;  foreach res-time-list
;  [ ?1 ->
;    file-type ?1
;    file-type "\r\n"
;  ]
;  file-close
;
;end