globals [
  ticks-per-second ; Ticks to make one simulated second, currently fixed at 100
  pressure ; Pressure of ambient air (Pa)
  temp ; Temperature of ambient air (K)
  particle-mass ; Mass of each simulated "particle," which moves like a single molecule but reperesents many. Inversely proportional to <num-particles>
  particles-per-tick ; Number of particles that should enter the view each tick, representing oncoming air
  skin-patches ; List of patches that make up the "skin" of the airfoil (colored black)
  max-tick-delta ; Maximum number of ticks that can pass in one cycle, currently fixed at 1
  tick-delta ; Number of ticks that will pass in this cycle, always less than or equal to <max-tick-delta>
  impulse-avg-ticks ; Number of ticks over which to average impulses to calculate forces. Inversely proportional to <num-particles>
  lift-impulses ; List of y-components of impulses applied to the airfoil over the last <impulse-avg-ticks> ticks
  drag-impulses ; List of x-components of impulses applied to the airfoil over the last <impulse-avg-ticks> ticks
  impulse-expires ; List of times (in ticks) which correspond to elements of <lift-impulses> and <drag-impulses>
                  ; Each time shows when an impulse "expires" and is removed from the list, <impulse-avg-ticks> after it occurs
  cycles ; Number of cycles of the go loop since the last setup, different from <ticks>
]

patches-own [
  exact-x ; For members of <skin-patches>, exact x-coordinate of the point on the airfoil's skin
  exact-y ; For members of <skin-patches>, exact y-coordinate of the point on the airfoil's skin
  prev ; For elements of <skin-patches>, adjacent skin patch in counter-clockwise direction
  next ; For elements of <skin-patches>, adjacent skin patch in clockwise direction
  skin-angle ; For elements of <skin-patches>, angle of airfoil skin at this patch
  particle-speeds ; For patches outside airfoil, list of speeds of particles that passed thorugh the patch in the last <impulse-avg-ticks> ticks
  particle-headings ; For patches outside airfoil, list of headings of particles that passed thorugh the patch in the last <impulse-avg-ticks> ticks
  particle-speed-expires ; For patches outside airfoil, list of times (in ticks) which correspond to elements of <particle-speeds> and <particle-headings>
                         ; Each time shows when a speed and heading "expires" and is removed from the list, <impulse-avg-ticks> after it occurs
  avg-velocity-heading ; For patches outside airfoil, heading of average velocity through this patch over the last <impulse-avg-ticks> ticks
                       ; Equivalent to streamline angle at this patch
  patch-pressure ; For patches outside airfoil, pressure at that specific patch
                 ; The number is about 3x too high for some reason, but it's correct relative to other patches
  streamline-turtle ; For patches with a streamline, a reference to the turtle representing it
]

breed [ particles particle ]
breed [ streamlines streamline ]

particles-own [
  speed ; Distance this particle will move per tick
  mass
  last-collision ; Last thing this particle collided with, can be patch or particle
]

to setup
  clear-all
  set-default-shape particles "circle"
  set-default-shape streamlines "line"
  ask patches [ ; Reset patch color and initialize lists for patch variables
    set pcolor white
    set particle-speeds (list)
    set particle-headings (list)
    set particle-speed-expires (list)
  ]
  set ticks-per-second 100
  set pressure 101325 * ((1 - (2.25577 * altitude * (10 ^ -5))) ^ 5.25588) ; Pressure in atmosphere, from
                                                                           ; http://www.engineeringtoolbox.com/air-altitude-pressure-d_462.html
  set temp 288 - (altitude * 6.5 / 1000) ; Temperature in troposhere, from http://www.srh.noaa.gov/jetstream/atmos/layers.htm
  let molar-mass .02897 ; Kilograms of air per mole of molecules (over all gases), from http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html
  let total-area (count patches) / 400 ; Total area of view in m^2, 400 patches = 1 m^2
  let gas-constant 8.314 ; Gas constant for ideal gas law
  set particle-mass (molar-mass * pressure * total-area) / (gas-constant * temp * num-particles) ; Caluclated from ideal gas law
  let y-length (max-pycor - min-pycor + 1) / 20 ; Length (m) of y-axis of view
  set particles-per-tick (num-particles * airspeed * y-length) / (total-area * ticks-per-second) ; Particle density times volume of entering air per tick
  set skin-patches (list) ; Initialize lists and constant
  set max-tick-delta 1
  set impulse-avg-ticks 5000 / num-particles ; Higher <impulse-avg-ticks> smooths data outputs more, only necessary with low numbers of particles
  set lift-impulses (list)
  set drag-impulses (list)
  set impulse-expires (list)
  if airfoil-length > 0 [ ; If airfoil is supposed to exists, draw airfoil using NACA-4 parameters
    create-naca-4 (-10 * airfoil-length) (10 * airfoil-length) (max-camber / 100) (max-camber-position / 100) (max-thickness / 100) angle-of-attack
  ]
  create-streamlines (count patches with [pcolor = white and (pxcor + 100) mod 5 = 2 and (pycor + 100) mod 5 = 2]) [ ; Create streamlines on a 5x5 grid (excluding airfoil)
                                                                                                                     ; Using mod 5 = 2 to avoid edges
    set color black
    set size 5
    set hidden? true ; Streamlines are hidden until initial data collection is done
    set xcor .01 ; xcor = ycor = .01 indicates streamline is unclaimed
    set ycor .01
  ]
  ask patches with [pcolor = white and (pxcor + 100) mod 5 = 2 and (pycor + 100) mod 5 = 2] [ ; For each patch that should have a streamline:
    set streamline-turtle one-of streamlines with [xcor = .01 and ycor = .01] ; Pick an unclaimed turtle
    ask streamline-turtle [ ; Move it to location of patch, claiming it
      set xcor [pxcor] of myself
      set ycor [pycor] of myself
    ]
  ]
  reset-ticks
end

to go
  create-new-particles
  calc-tick-delta ; See how many ticks will pass in this cycle
  ask particles [ move ] ; Move right after calculating <tick-delta> to make sure other changes do not invalidate it
  ask particles [
    if patch-here != patch-ahead (-1 * speed * tick-delta) [ ; If the turtle is on a different patch than it was last cycle:
      ask patch-here [
        set particle-speeds lput ([speed] of myself) particle-speeds ; Register that it passed through the new patch, for pressure and streamline purposes
        set particle-headings lput ([heading] of myself) particle-headings
        set particle-speed-expires lput (ticks + impulse-avg-ticks) particle-speed-expires
      ]
    ]
  ]
  ask particles [ check-for-collision ] ; After moving and registering movements, process collisions
  tick-advance tick-delta ; Move tick counter forward by <tick-delta>, usually a fraction of a tick
  if floor ticks > floor (ticks - tick-delta) [ ; If this is a new tick (this code should run once per tick, not cycle):
    ask patches [ update-overlays ] ; Prepare patch data and show new patch pressures and stremlines
    prepare-plot-data ; Ensure impulse lists are up-to-date
    update-plots
  ]
  set cycles cycles + 1
end

to trace-particles
  ifelse count (particles with [not hidden?]) >= num-trace-particles; If there are at least <num-trace-particles> visible:
  [
    ask n-of num-trace-particles (particles with [not hidden?]) [ pen-down ] ; Trace <num-trace-particles> visible particles
  ]
  [
    ask particles with [not hidden?] [ pen-down ] ; Otherwise, trace all visible particles
  ]
end

to create-naca-4 [startx endx m p t aoa]
  ; Airfoil formulas from http://people.clarkson.edu/~pmarzocc/AE429/The%20NACA%20airfoil%20series.pdf
  let per-patch 100 ; Number of points to calculate per patch, higher numbers take longer to setup but lower numbers might leave holes in the airfoil
  let x startx
  let bot-patches (list) ; Initialize list of patches on the bottom of the airfoil
  repeat (endx - startx) * per-patch + 1 [ ; For every point between the start and end patches:
    let xfrac (x - startx) / (endx - startx) ; Fraction of way to airfoil end, used in formulas
    let yfrac 0 ; y-coordinate of camber line scaled at same factor as xfrac
    let dydx 0 ; Derivative of camber line at this point
    ifelse xfrac <= p [
      set yfrac (m / (p ^ 2) * (2 * p * xfrac - xfrac ^ 2)) ; From formula paper
      set dydx (m / (p ^ 2) * (2 * p - 2 * xfrac))
    ]
    [
      set yfrac (m / ((1 - p) ^ 2) * (1 - 2 * p + 2 * p * xfrac - xfrac ^ 2))
      set dydx (m / ((1 - p) ^ 2) * (2 * p - 2 * xfrac))
    ]
    let y yfrac * (endx - startx) ; Scale <yfrac> back to normal coordinates
    ask cors-to-patch rotated-cors (list x y) (list 0 0) aoa [ if not member? self skin-patches [ set pcolor gray ] ] ; Draw camber line patch in gray
    let angle ((atan 1 dydx) - 90) mod 360 ; Angle of camber line based on derivative
    ; Thickness distribution formula from paper
    let thick (t / 0.2) * (0.2969 * (sqrt xfrac) - 0.1260 * xfrac - 0.3516 * (xfrac ^ 2) + 0.2843 * (xfrac ^ 3) - 0.1015 * (xfrac ^ 4)) * (endx - startx)
    let top-cors rotated-cors (cors-at-vector (list x y) thick angle) (list 0 0) aoa ; Skin coordinates are <thick> and <-thick> distance at <angle> from camber line patch
    let bot-cors rotated-cors (cors-at-vector (list x y) (-1 * thick) angle) (list 0 0) aoa
    let top-patch cors-to-patch top-cors ; Convert coordinates to patches
    let bot-patch cors-to-patch bot-cors
    if not member? top-patch skin-patches [ ; If top patch is not already on the skin:
      ask (top-patch) [
        set pcolor black
        set exact-x item 0 top-cors ; Add exact coordinates as patch properties
        set exact-y item 1 top-cors
      ]
      set skin-patches lput top-patch skin-patches ; Add patch to end of list of skin patches
    ]
    if not member? bot-patch skin-patches [ ; Same for bottom patch
      ask (bot-patch) [
        set pcolor black
        set exact-x item 0 bot-cors
        set exact-y item 1 bot-cors
        set bot-patches lput self bot-patches
      ]
      set skin-patches fput bot-patch skin-patches ; Patch is added to beginning instead of end to keep the list ordered clockwise
    ]
    set x x + (1 / per-patch) ; Move to next point
  ]
  ask patch-set skin-patches [ calc-skin-angle ] ; Calculate skin angle to be used in collisions with particles
  ask patch-set bot-patches [ fill-off-white ] ; Fill patches inside airfoil with a different white color so they are not considered "empty"
end

to-report rotated-cors [cors center angle] ; Helper to rotate coordinates around a center point
  let cx item 0 center
  let cy item 1 center
  let relx item 0 cors - cx
  let rely item 1 cors - cy
  ; Negative angle is used because of nonstandard coordinate system
  report list (cx + relx * (cos (-1 * angle)) - rely * (sin (-1 * angle))) (cy + relx * (sin (-1 * angle)) + rely * (cos (-1 * angle)))
end

to-report cors-at-vector [cors mag angle] ; Helper for vector addition
  ; Sine and cosine are reversed because NetLogo assigns angle 0 to +y
  report list ((item 0 cors) + mag * (sin angle)) ((item 1 cors) + mag * (cos angle))
end

to-report cors-to-patch [cors]
  report patch (item 0 cors) (item 1 cors)
end

to calc-skin-angle ; Patch procedure to calculate skin angle to be used in collisions with particles
  let pos position self skin-patches ; Find position of this patch in the list
  ifelse pos = 0 or pos = (length skin-patches - 1) [ ; Find previous and next patches in a clockwise ordering
    ifelse pos = 0 [ ; If at the beginnning of the list, rollover <prev> to the end
      set prev item (length skin-patches - 1) skin-patches
      set next item 1 skin-patches
    ]
    [ ; If at the end of the list, rollover <next> to the beginning
      set prev item (length skin-patches - 2) skin-patches
      set next item 0 skin-patches
    ]
  ]
  [ ; Otherwise, assign <prev> and <next> normally
    set prev item (pos - 1) skin-patches
    set next item (pos + 1) skin-patches
  ]
  set skin-angle atan ([exact-x] of next - [exact-x] of prev) ([exact-y] of next - [exact-y] of prev) ; Calculate angle of tangent line between adjacent patches
end

to fill-off-white ; Patch procedure to fill patches inside airfoil with a different white color so they are not considered "empty"
  let i 1 ; Start at adjacent upward patch
  let to-fill-list (list) ; Keep track of all patches we intend to fill
  let to-fill patch-at 0 1
  while [to-fill != nobody and [pcolor] of to-fill != black] [ ; While we're looking at an empty patch:
    set to-fill-list lput to-fill to-fill-list ; Add it to the list we intend to fill
    set i i + 1
    set to-fill patch-at 0 i ; Look at the next upward patch
  ]
  if to-fill != nobody [ ; Only fill if we encounter the top of the airfoil, rather than the edge of the view
    ask patch-set to-fill-list [
      if pcolor = white [ set pcolor 39.99 ] ; Color 39.99 is extremely light brown, which looks like white but isn't
    ]
  ]
end

to create-new-particles
  create-particles min list (round (2 * particles-per-tick)) (num-particles - count particles) [ ; Create particles to replace up to <num-particles>
                                                                                                 ; with a maximum of twice the expected number
    let gas-constant 8.314
    let molar-mass .02897
    set speed (sqrt (3 * gas-constant * temp / molar-mass )) * 20 / ticks-per-second ; Root-mean-square speed of air based on environmental parameters
    set heading random 360 ; Root-mean-square speed has random direction
    add-speed-heading (airspeed * 20 / ticks-per-second) 90 ; Vector add airspeed to root-mean-square speed
    set mass particle-mass
    set last-collision nobody
    ifelse sin heading >= 0
    [ set xcor min-pxcor ] ; If it's going forward, start on the incoming edge
    [ set xcor max-pxcor ] ; Otherwise, start on the outgoing edge
    set ycor (min-pycor + random (max-pycor - min-pycor + 1)) ; Random y-position
    set color orange
    set hidden? (count particles with [not hidden?] >= num-visible-particles) ; If we need a new visible particle, make this one visible
  ]
end

to add-speed-heading [add-speed add-heading] ; Particle procedure for vector addition of velocities
  let current-speed-x (speed * sin heading) ; Break speeds into components
  let current-speed-y (speed * cos heading)
  let add-speed-x (add-speed * sin add-heading)
  let add-speed-y (add-speed * cos add-heading)
  let new-speed-x current-speed-x + add-speed-x ; Sum speeds in each direction
  let new-speed-y current-speed-y + add-speed-y
  set speed sqrt (new-speed-x ^ 2 + new-speed-y ^ 2)
  ifelse speed = 0 ; If the particle has no speed, heading is undefined, so just set it to 0
  [ set heading 0 ]
  [ set heading atan new-speed-x new-speed-y ]
end

to calc-tick-delta
  ; From GasLab Free Gas
  ; tick-delta is calculated in such way that even the fastest
  ; particle will jump at most 1 patch length in a tick. As
  ; particles jump (speed * tick-delta) at every tick, making
  ; tick length the inverse of the speed of the fastest particle
  ; (1/max speed) assures that. Having each particle advance at most
  ; one patch-length is necessary for them not to jump over each other
  ; without colliding.
  ifelse any? particles with [speed > 0]
  [ set tick-delta min list (1 / (ceiling max [speed] of particles)) max-tick-delta ]
  [ set tick-delta max-tick-delta ]
end

to move ; Particle procedure
  if patch-ahead (speed * tick-delta) != patch-here ; Reset <last-collision> if this particle is moving to a new patch
  [ set last-collision nobody ]
  jump (speed * tick-delta) ; Since one cycle may be less than a tick, multiply speed by the length of a cycle to calculate movement
end

to check-for-collision ; Particle procedure
  if [pcolor] of patch-here = gray [ die ] ; Get rid of particles that glitch into airfoil
  if patch-here != last-collision [
    if [pxcor] of patch-here = min-pxcor and sin heading < 0 ; If this particle is near a wall and moving toward it, collide with it
    [ one-particle-collide 90 ]
    if [pxcor] of patch-here = max-pxcor and sin heading > 0 ; If this particle collides with the outgoing edge, remove it
    [ die ]
    if [pycor] of patch-here = min-pycor and cos heading < 0
    [ one-particle-collide 0 ]
    if [pycor] of patch-here = max-pycor and cos heading > 0
    [ one-particle-collide 180 ]
    if member? patch-here skin-patches [ ; If this particle collides with the airfoil:
      impulse-airfoil (2 * mass * speed-on-vector (([skin-angle] of patch-here + 90) mod 360)) (([skin-angle] of patch-here + 90) mod 360) ; Apply impulse to airfoil
                                                                                                                                           ; based on expected change in
                                                                                                                                           ; momentum
      one-particle-collide (([skin-angle] of patch-here + 90) mod 360) ; Process effect of collision on particle
    ]
  ]
  
  ; From GasLab Free Gas
  ;; Here we impose a rule that collisions only take place when there
  ;; are exactly two particles per patch.
  if count other particles-here = 1
  [
    ;; the following conditions are imposed on collision candidates:
    ;;   1. they must have a lower who number than my own, because collision
    ;;      code is asymmetrical: it must always happen from the point of view
    ;;      of just one particle.
    ;;   2. they must not be the same particle that we last collided with on
    ;;      this patch, so that we have a chance to leave the patch after we've
    ;;      collided with someone.
    let candidate one-of other particles-here with
      [who < [who] of myself and myself != last-collision]
    ;; we also only collide if one of us has non-zero speed. It's useless
    ;; (and incorrect, actually) for two particles with zero speed to collide.
    if (candidate != nobody) and (speed > 0 or [speed] of candidate > 0)
    [
      two-particle-collide candidate ; Process collision
      set last-collision candidate ; Change <last-collision> to reflect collision
      ask candidate [ set last-collision myself ]
    ]
  ]
end

to one-particle-collide [collide-vector] ; Particle procedure for particle colliding with static object
  let v-par speed-on-vector collide-vector ; Break particle velocity in components along collision vector
  let v-perp speed-perp-vector collide-vector
  if speed > 0
  [ set heading collide-vector - (atan v-perp (-1 * v-par)) ] ; Reverse parallel velocity when calculating new heading
                                                              ; No need to change speed, since energy and therefore speed in conserved
  set last-collision patch-here ; Change <last-collision> to reflect collision
end

; From GasLab Free Gas
;; implements a collision with another particle.
;;
;; THIS IS THE HEART OF THE PARTICLE SIMULATION, AND YOU ARE STRONGLY ADVISED
;; NOT TO CHANGE IT UNLESS YOU REALLY UNDERSTAND WHAT YOU'RE DOING!
;;
;; The two particles colliding are self and other-particle, and while the
;; collision is performed from the point of view of self, both particles are
;; modified to reflect its effects. This is somewhat complicated, so I'll
;; give a general outline here:
;;   1. Do initial setup, and determine the heading between particle centers
;;      (call it theta).
;;   2. Convert the representation of the velocity of each particle from
;;      speed/heading to a theta-based vector whose first component is the
;;      particle's speed along theta, and whose second component is the speed
;;      perpendicular to theta.
;;   3. Modify the velocity vectors to reflect the effects of the collision.
;;      This involves:
;;        a. computing the velocity of the center of mass of the whole system
;;           along direction theta
;;        b. updating the along-theta components of the two velocity vectors.
;;   4. Convert from the theta-based vector representation of velocity back to
;;      the usual speed/heading representation for each particle.
;;   5. Perform final cleanup and update derived quantities.
to two-particle-collide [other-particle] ; Particle procedure
  ;;; PHASE 1: initial setup

  ;; for convenience, grab some quantities from other-particle
  let mass2 [mass] of other-particle
  let speed2 [speed] of other-particle
  let heading2 [heading] of other-particle

  ;; since particles are modeled as zero-size points, theta isn't meaningfully
  ;; defined. we can assign it randomly without affecting the model's outcome.
  let theta (random-float 360)



  ;;; PHASE 2: convert velocities to theta-based vector representation

  ;; now convert my velocity from speed/heading representation to components
  ;; along theta and perpendicular to theta
  let v1t (speed * cos (theta - heading))
  let v1l (speed * sin (theta - heading))

  ;; do the same for other-particle
  let v2t (speed2 * cos (theta - heading2))
  let v2l (speed2 * sin (theta - heading2))



  ;;; PHASE 3: manipulate vectors to implement collision

  ;; compute the velocity of the system's center of mass along theta
  let vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2) )

  ;; now compute the new velocity for each particle along direction theta.
  ;; velocity perpendicular to theta is unaffected by a collision along theta,
  ;; so the next two lines actually implement the collision itself, in the
  ;; sense that the effects of the collision are exactly the following changes
  ;; in particle velocity.
  set v1t (2 * vcm - v1t)
  set v2t (2 * vcm - v2t)



  ;;; PHASE 4: convert back to normal speed/heading

  ;; now convert my velocity vector into my new speed and heading
  set speed sqrt ((v1t ^ 2) + (v1l ^ 2))
  ;; if the magnitude of the velocity vector is 0, atan is undefined. but
  ;; speed will be 0, so heading is irrelevant anyway. therefore, in that
  ;; case we'll just leave it unmodified.
  if v1l != 0 or v1t != 0
    [ set heading (theta - (atan v1l v1t)) ]

  ;; and do the same for other-particle
  ask other-particle [
    set speed sqrt ((v2t ^ 2) + (v2l ^ 2))
    if v2l != 0 or v2t != 0
      [ set heading (theta - (atan v2l v2t)) ]
  ]
end

to-report speed-on-vector [angle] ; Particle procedure, helper for vector math
  report speed * cos (angle - heading)
end

to-report speed-perp-vector [angle] ; Particle procedure, helper for vector math
  report speed * sin (angle - heading)
end

to impulse-airfoil [impulse angle] ; Apply an impulse to be calculated into airfoil forces
  set lift-impulses lput (impulse * cos angle) lift-impulses ; Add components of impulse into corresponding lists
  set drag-impulses lput (impulse * sin angle) drag-impulses
  set impulse-expires lput (ticks + impulse-avg-ticks) impulse-expires ; Set expire time for impulses
end

to update-overlays ; Patch procedure to show patch pressure and streamlines
  while [(length particle-speed-expires) > 0 and ticks > (first particle-speed-expires)] [ ; While there are some expired impulses on this patch:
    set particle-speeds but-first particle-speeds ; Remove them from the impulse lists
    set particle-headings but-first particle-headings
    set particle-speed-expires but-first particle-speed-expires ; Remove their expiration times from the expiration list
  ]
  if length particle-speeds > 0 and (pcolor = white or shade-of? sky pcolor) and ticks > impulse-avg-ticks [ ; If some particles have passed through this patch, this patch
                                                                                                             ; is empty, and we are ready to display data:
    let x-speed-sum 0 ; Create variables to store total velocities in each direction
    let y-speed-sum 0
    foreach n-values (length particle-speeds) [?] [ ; For every index in the particle speed and heading lists...
      set x-speed-sum (x-speed-sum + (item ? particle-speeds * (sin item ? particle-headings))) ; Add components of speed and heading at that index to the sums
      set y-speed-sum (y-speed-sum + (item ? particle-speeds * (cos item ? particle-headings)))
    ]
    set avg-velocity-heading atan x-speed-sum y-speed-sum ; Heading of average velocity is the same as heading of summed velocity
    let perp-impulse 0 ; Sum of hypothetical impulses perpendicular to <avg-velocity-heading>
    foreach n-values (length particle-speeds) [?] [ ; For every index in the particle speed and heading lists...
      let perp-speed ((item ? particle-speeds) * sin (avg-velocity-heading - (item ? particle-headings))) ; Find particle speed perpendicular to <avg-velocity-heading>
      if perp-speed > 0 [
        set perp-impulse (perp-impulse + (perp-speed * particle-mass)) ; Add impulse that would be applied to a solid surface along <avg-velocity-heading> to sum
      ]
    ]
    set patch-pressure (perp-impulse / (impulse-avg-ticks / ticks-per-second * .0025)) ; Pressure is impulse per unit time per unit area
    ifelse show-pressure? [
      set pcolor scale-color sky patch-pressure (8 * pressure) 0 ; Scale <patch-pressure> into color value, where higher is darker
    ]
    [ set pcolor white ]
    if (pxcor + 100) mod 5 = 2 and (pycor + 100) mod 5 = 2 and ticks > impulse-avg-ticks + 1 [ ; If this patch has a streamline:
      ifelse show-streamlines? [
        ask streamline-turtle [
          set hidden? false
          set heading mean [avg-velocity-heading] of (neighbors with [pcolor = white or shade-of? sky pcolor]) ; Average <avg-velocity-heading> from neighboring patches
                                                                                                               ; for smoothing
        ]
      ]
      [ ask streamline-turtle [ set hidden? true ] ]
    ]
  ]
end

to prepare-plot-data ; Remove expired impulses from airfoil impulse lists before updating force graphs
  while [length impulse-expires > 0 and ticks > (first impulse-expires)] [ ; While there are some expired impulses on the airfoil:
    set lift-impulses but-first lift-impulses ; Remove them from the impulse lists
    set drag-impulses but-first drag-impulses
    set impulse-expires but-first impulse-expires ; Remove their expiration times from the expiration list
  ]
end
@#$#@#$#@
GRAPHICS-WINDOW
9
11
473
496
50
50
4.5
1
10
1
1
1
0
0
0
1
-50
50
-50
50
1
1
1
ticks
30.0

BUTTON
486
12
552
45
NIL
setup
NIL
1
T
OBSERVER
NIL
S
NIL
NIL
1

BUTTON
570
12
633
45
NIL
go
T
1
T
OBSERVER
NIL
G
NIL
NIL
0

SLIDER
486
80
658
113
altitude
altitude
0
15000
9000
500
1
meters
HORIZONTAL

TEXTBOX
507
59
657
77
Environmental options
11
0.0
1

SLIDER
486
121
658
154
airspeed
airspeed
0
400
220
10
1
m/s
HORIZONTAL

TEXTBOX
530
206
625
224
Airfoil options
11
0.0
1

SLIDER
486
267
689
300
max-camber
max-camber
0
9
2
0.5
1
% of length
HORIZONTAL

SLIDER
486
226
663
259
airfoil-length
airfoil-length
0
4.5
2
.1
1
meters
HORIZONTAL

SLIDER
486
308
747
341
max-camber-position
max-camber-position
0
99
40
1
1
% of length
HORIZONTAL

SLIDER
486
349
702
382
max-thickness
max-thickness
0
99
15
1
1
% of length
HORIZONTAL

SLIDER
486
391
689
424
angle-of-attack
angle-of-attack
-50
50
20
1
1
degrees
HORIZONTAL

SLIDER
486
163
691
196
num-particles
num-particles
0
1000
250
10
1
particles
HORIZONTAL

TEXTBOX
797
14
947
32
Visualization options
11
0.0
1

SLIDER
763
34
950
67
num-visible-particles
num-visible-particles
0
500
100
10
1
NIL
HORIZONTAL

SLIDER
763
76
942
109
num-trace-particles
num-trace-particles
0
100
1
1
1
NIL
HORIZONTAL

BUTTON
894
117
1000
150
clear-traces
clear-drawing
NIL
1
T
OBSERVER
NIL
C
NIL
NIL
0

PLOT
763
201
963
425
Lift of airfoil
Time (ticks)
Lift per meter of width (N)
0.0
10.0
0.0
10.0
true
false
"" "if length lift-impulses = 0 or ticks < impulse-avg-ticks [stop]"
PENS
"default" 1.0 0 -16777216 true "" "plot (sum lift-impulses) / (impulse-avg-ticks / ticks-per-second)"

PLOT
975
201
1175
425
Drag of airfoil
Time (ticks)
Drag per meter of width (N)
0.0
10.0
0.0
10.0
true
false
"" "if length drag-impulses = 0 or ticks < impulse-avg-ticks [stop]"
PENS
"default" 1.0 0 -16777216 true "" "plot (sum drag-impulses) / (impulse-avg-ticks / ticks-per-second)"

TEXTBOX
763
430
953
457
Note: 100 ticks = 1 second
11
0.0
1

BUTTON
763
117
886
150
NIL
trace-particles
NIL
1
T
OBSERVER
NIL
T
NIL
NIL
0

SWITCH
763
158
920
191
show-pressure?
show-pressure?
0
1
-1000

SWITCH
929
158
1102
191
show-streamlines?
show-streamlines?
0
1
-1000

TEXTBOX
975
431
1148
487
Some displays are disabled for (5000 / num-particles) ticks to allow initial data collection
11
0.0
1

@#$#@#$#@
## WHAT IS IT?

This model simulates in two dimensions the motion of air around an airfoil shape, such as an airplane wing. It uses the motion of individual particles according to the Ideal Gas Law and elastic collisions, which are similar to the behavior of molecules in the real world.

This model can be used to analyze the effects of different airfoil shapes or environmental conditions, such as altitude or airspeed, on the forces on the airfoil and the turbulence patterns it creates. It also shows how the general principles of fluid dynamics come from the motion of individual particles.

## HOW IT WORKS

Particles begin on the left side of the screen, with a velocity consisting of the vector sum of the airspeed, directed to the right, and a random velocity determined by the temperature. They move in straight lines until the collide with walls, the airfoil, or other particles.

When two particles collide, they exchange velocity along a random vector, and thereby change direction. They bounce off walls and the airfoil like balls reflecting off a solid surface, but when they collide with the airfoil, the impulse they apply to it is tracked in order to to calculate the overall lift and drag. When particles reach the right edge of the view, they are removed from the simulation.

The motion of particles through each patch is also tracked to provide information for the overlays, as described under Visualization Options below.

## HOW TO USE IT

To run the model, press the `setup` button, then press the `go` button. Press the `go` button again to pause. After changing any options other than the Visualization Options, press the `setup` button to reset the model.

#### Environmental options

- `altitude` The altitude this airfoil is flying at, which affects the ambient pressure and temperature. Normal cruising altitude is approximately 9000 meters.

- `airspeed` The speed the airfoil is flying at, which is simulated as the average velocity of the air entering the view. Normal cruising speed is approximately 220 meters per second.

- `num-particles` Number of particles the model simulates at a time. This should not affect overall results, but higher values will produce more stable output at the expense of slower updates.

#### Airfoil options

These options specify the airfoil using shapes in the NACA-4 series, see the last link in the Credits and References.

- `airfoil-length` The distance from the leading edge of the airfoil to its trailing edge.

- `max-camber` The maximum deviation from a straight line in the airfoil's camber line (shown in gray), specified as a percent of `airfoil-length`.

- `max-camber-position` The position of the camber line's maximum deviation, specified as a percent of `airfoil-length`.

- `max-thickness` Maximum distance the airfoil skin can be from the camber line, specified as a percent of `airfoil-length`.

- `angle-of-attack` Rotation of the entire airfoil, relative to the orientation where the two ends of the camber line lie on a horizontal line.

#### Visualization options

- `num-visible-particles` How many particles appear in the view, must be less than or equal to `num-particles`.

- `num-trace-particles` How many particles to trace when `trace-particles` is clicked.

- `trace-particles` Ask `num-trace-particles` visible particles to enable their pens, meaning they will leave orange lines showing their paths until they leave the view.

- `clear-traces` Clear all lines from traced particles.

- `show-pressure?` Whether to show the overlay of pressure at each patch, which is enabled after `5000 / num-particles` ticks. Patches are colored shades of blue, where darker shades are higher pressures.

- `show-streamlines?` Whether to show the overlay of streamlines, which is enabled after `5000 / num-particles` ticks. The angle of each streamline shows the direction of the average velocity of particles near it.

#### Plots

- `Lift of airfoil` Shows the lift that would be applied to one meter of the displayed airfoil shape if its was in three-dimensional space. Values will fluctuate, but should approach a long-term average. Enabled after `5000 / num-particles` ticks.

- `Drag of airfoil` Shows the drag that would be applied to one meter of the displayed airfoil shape if its was in three-dimensional space. Values will fluctuate, but should approach a long-term average. Enabled after `5000 / num-particles` ticks.


## THINGS TO NOTICE

Notice how the streamlines have an upward slope behind the trailing edge, showing the "uplift" that is a significant cause of turbulence behind wings.

Trace a few particles. Notice how their motion does not closely reflect the overall motion of the air as shown by the streamlines.

Notice how the airfoil changes the pressure and streamlines for a quite large area around it.

Notice that due to the limited number of particles we can simulate, the forces on the airfoil fluctuate significantly, while in the real world they would be much steadier.

Notice how variations in the lift and drag plots tend to roughly track each other, possibly because of variations in the overall number of particles hitting the airfoil at any particular time.

## THINGS TO TRY

Change the altitude slider. What happens to the forces on the airfoil? Why do airplanes fly at high altitude?

Change the airspeed slider. What happens to the forces on the airfoil? How could this be used to determine the speed a particular airplane needs to take off?

Change the airfoil's angle of attack. What tradeoffs might this decision bring?

Try designing your own airfoils with the Airfoil Options. What seems like a good design? How might this vary for different kinds of airplanes?

## EXTENDING THE MODEL

This model doesn't properly deal with supersonic effects. Why not? Implement the speed of sound. How does this change afffect supersonic flight?

Add a more detailed view of air motion than the five-patch streamlines offer. What new effects are visible? Is there turbulence near the airfoil's trailing edge?

Make this model into an evolutionary design where a generation of random airfoils is tested, and the ones which perform best are more likely to pass their parameters to the next generation. How do you define which airfoils perform best? How does this affect the results after many generations?

## NETLOGO FEATURES

This model uses the `tick-advance` feature to divide time in fractions of ticks, ensuring that particles can never glitch through barriers by moving more than one patch per processing cycle, no matter their speed. This technique comes from the GasLab models.

As NetLogo does not have a data structure similar to a map, corresponding lists were used instead, where the first element of one list is paired with the first element of the other, the second is paired with the second, and so on. This structure can be seen with the airfoil impulse lists and their expiration times, and with the patch speed/heading lists and their expiration times.

## RELATED MODELS

The GasLab example models under the Chemistry and Physics section of the Models Library provided some of the particle collision code for this model and explore in more detail various aspects of gas behavior.

## CREDITS AND REFERENCES

This model was created as a final project for John Pickle's course "Computer Modeling in Science" at Concord Academy. His advice and suggestions are gratefully acknowledged.

The core particle motion code comes from the GasLab Free Gas example model, by Uri Wilensky. Algorithms used in this model come from:

http://www.engineeringtoolbox.com/air-altitude-pressure-d_462.html
http://www.srh.noaa.gov/jetstream/atmos/layers.htm
http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html
http://people.clarkson.edu/~pmarzocc/AE429/The%20NACA%20airfoil%20series.pdf
@#$#@#$#@
default
true
0
Polygon -7500403 true true 150 5 40 250 150 205 260 250

airplane
true
0
Polygon -7500403 true true 150 0 135 15 120 60 120 105 15 165 15 195 120 180 135 240 105 270 120 285 150 270 180 285 210 270 165 240 180 180 285 195 285 165 180 105 180 60 165 15

arrow
true
0
Polygon -7500403 true true 150 0 0 150 105 150 105 293 195 293 195 150 300 150

box
false
0
Polygon -7500403 true true 150 285 285 225 285 75 150 135
Polygon -7500403 true true 150 135 15 75 150 15 285 75
Polygon -7500403 true true 15 75 15 225 150 285 150 135
Line -16777216 false 150 285 150 135
Line -16777216 false 150 135 15 75
Line -16777216 false 150 135 285 75

bug
true
0
Circle -7500403 true true 96 182 108
Circle -7500403 true true 110 127 80
Circle -7500403 true true 110 75 80
Line -7500403 true 150 100 80 30
Line -7500403 true 150 100 220 30

butterfly
true
0
Polygon -7500403 true true 150 165 209 199 225 225 225 255 195 270 165 255 150 240
Polygon -7500403 true true 150 165 89 198 75 225 75 255 105 270 135 255 150 240
Polygon -7500403 true true 139 148 100 105 55 90 25 90 10 105 10 135 25 180 40 195 85 194 139 163
Polygon -7500403 true true 162 150 200 105 245 90 275 90 290 105 290 135 275 180 260 195 215 195 162 165
Polygon -16777216 true false 150 255 135 225 120 150 135 120 150 105 165 120 180 150 165 225
Circle -16777216 true false 135 90 30
Line -16777216 false 150 105 195 60
Line -16777216 false 150 105 105 60

car
false
0
Polygon -7500403 true true 300 180 279 164 261 144 240 135 226 132 213 106 203 84 185 63 159 50 135 50 75 60 0 150 0 165 0 225 300 225 300 180
Circle -16777216 true false 180 180 90
Circle -16777216 true false 30 180 90
Polygon -16777216 true false 162 80 132 78 134 135 209 135 194 105 189 96 180 89
Circle -7500403 true true 47 195 58
Circle -7500403 true true 195 195 58

circle
false
0
Circle -7500403 true true 0 0 300

circle 2
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240

cow
false
0
Polygon -7500403 true true 200 193 197 249 179 249 177 196 166 187 140 189 93 191 78 179 72 211 49 209 48 181 37 149 25 120 25 89 45 72 103 84 179 75 198 76 252 64 272 81 293 103 285 121 255 121 242 118 224 167
Polygon -7500403 true true 73 210 86 251 62 249 48 208
Polygon -7500403 true true 25 114 16 195 9 204 23 213 25 200 39 123

cylinder
false
0
Circle -7500403 true true 0 0 300

dot
false
0
Circle -7500403 true true 90 90 120

face happy
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 255 90 239 62 213 47 191 67 179 90 203 109 218 150 225 192 218 210 203 227 181 251 194 236 217 212 240

face neutral
false
0
Circle -7500403 true true 8 7 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Rectangle -16777216 true false 60 195 240 225

face sad
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 168 90 184 62 210 47 232 67 244 90 220 109 205 150 198 192 205 210 220 227 242 251 229 236 206 212 183

fish
false
0
Polygon -1 true false 44 131 21 87 15 86 0 120 15 150 0 180 13 214 20 212 45 166
Polygon -1 true false 135 195 119 235 95 218 76 210 46 204 60 165
Polygon -1 true false 75 45 83 77 71 103 86 114 166 78 135 60
Polygon -7500403 true true 30 136 151 77 226 81 280 119 292 146 292 160 287 170 270 195 195 210 151 212 30 166
Circle -16777216 true false 215 106 30

flag
false
0
Rectangle -7500403 true true 60 15 75 300
Polygon -7500403 true true 90 150 270 90 90 30
Line -7500403 true 75 135 90 135
Line -7500403 true 75 45 90 45

flower
false
0
Polygon -10899396 true false 135 120 165 165 180 210 180 240 150 300 165 300 195 240 195 195 165 135
Circle -7500403 true true 85 132 38
Circle -7500403 true true 130 147 38
Circle -7500403 true true 192 85 38
Circle -7500403 true true 85 40 38
Circle -7500403 true true 177 40 38
Circle -7500403 true true 177 132 38
Circle -7500403 true true 70 85 38
Circle -7500403 true true 130 25 38
Circle -7500403 true true 96 51 108
Circle -16777216 true false 113 68 74
Polygon -10899396 true false 189 233 219 188 249 173 279 188 234 218
Polygon -10899396 true false 180 255 150 210 105 210 75 240 135 240

house
false
0
Rectangle -7500403 true true 45 120 255 285
Rectangle -16777216 true false 120 210 180 285
Polygon -7500403 true true 15 120 150 15 285 120
Line -16777216 false 30 120 270 120

leaf
false
0
Polygon -7500403 true true 150 210 135 195 120 210 60 210 30 195 60 180 60 165 15 135 30 120 15 105 40 104 45 90 60 90 90 105 105 120 120 120 105 60 120 60 135 30 150 15 165 30 180 60 195 60 180 120 195 120 210 105 240 90 255 90 263 104 285 105 270 120 285 135 240 165 240 180 270 195 240 210 180 210 165 195
Polygon -7500403 true true 135 195 135 240 120 255 105 255 105 285 135 285 165 240 165 195

line
true
0
Line -7500403 true 150 0 150 300

line half
true
0
Line -7500403 true 150 0 150 150

pentagon
false
0
Polygon -7500403 true true 150 15 15 120 60 285 240 285 285 120

person
false
0
Circle -7500403 true true 110 5 80
Polygon -7500403 true true 105 90 120 195 90 285 105 300 135 300 150 225 165 300 195 300 210 285 180 195 195 90
Rectangle -7500403 true true 127 79 172 94
Polygon -7500403 true true 195 90 240 150 225 180 165 105
Polygon -7500403 true true 105 90 60 150 75 180 135 105

plant
false
0
Rectangle -7500403 true true 135 90 165 300
Polygon -7500403 true true 135 255 90 210 45 195 75 255 135 285
Polygon -7500403 true true 165 255 210 210 255 195 225 255 165 285
Polygon -7500403 true true 135 180 90 135 45 120 75 180 135 210
Polygon -7500403 true true 165 180 165 210 225 180 255 120 210 135
Polygon -7500403 true true 135 105 90 60 45 45 75 105 135 135
Polygon -7500403 true true 165 105 165 135 225 105 255 45 210 60
Polygon -7500403 true true 135 90 120 45 150 15 180 45 165 90

sheep
false
15
Circle -1 true true 203 65 88
Circle -1 true true 70 65 162
Circle -1 true true 150 105 120
Polygon -7500403 true false 218 120 240 165 255 165 278 120
Circle -7500403 true false 214 72 67
Rectangle -1 true true 164 223 179 298
Polygon -1 true true 45 285 30 285 30 240 15 195 45 210
Circle -1 true true 3 83 150
Rectangle -1 true true 65 221 80 296
Polygon -1 true true 195 285 210 285 210 240 240 210 195 210
Polygon -7500403 true false 276 85 285 105 302 99 294 83
Polygon -7500403 true false 219 85 210 105 193 99 201 83

square
false
0
Rectangle -7500403 true true 30 30 270 270

square 2
false
0
Rectangle -7500403 true true 30 30 270 270
Rectangle -16777216 true false 60 60 240 240

star
false
0
Polygon -7500403 true true 151 1 185 108 298 108 207 175 242 282 151 216 59 282 94 175 3 108 116 108

target
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
Circle -7500403 true true 60 60 180
Circle -16777216 true false 90 90 120
Circle -7500403 true true 120 120 60

tree
false
0
Circle -7500403 true true 118 3 94
Rectangle -6459832 true false 120 195 180 300
Circle -7500403 true true 65 21 108
Circle -7500403 true true 116 41 127
Circle -7500403 true true 45 90 120
Circle -7500403 true true 104 74 152

triangle
false
0
Polygon -7500403 true true 150 30 15 255 285 255

triangle 2
false
0
Polygon -7500403 true true 150 30 15 255 285 255
Polygon -16777216 true false 151 99 225 223 75 224

truck
false
0
Rectangle -7500403 true true 4 45 195 187
Polygon -7500403 true true 296 193 296 150 259 134 244 104 208 104 207 194
Rectangle -1 true false 195 60 195 105
Polygon -16777216 true false 238 112 252 141 219 141 218 112
Circle -16777216 true false 234 174 42
Rectangle -7500403 true true 181 185 214 194
Circle -16777216 true false 144 174 42
Circle -16777216 true false 24 174 42
Circle -7500403 false true 24 174 42
Circle -7500403 false true 144 174 42
Circle -7500403 false true 234 174 42

turtle
true
0
Polygon -10899396 true false 215 204 240 233 246 254 228 266 215 252 193 210
Polygon -10899396 true false 195 90 225 75 245 75 260 89 269 108 261 124 240 105 225 105 210 105
Polygon -10899396 true false 105 90 75 75 55 75 40 89 31 108 39 124 60 105 75 105 90 105
Polygon -10899396 true false 132 85 134 64 107 51 108 17 150 2 192 18 192 52 169 65 172 87
Polygon -10899396 true false 85 204 60 233 54 254 72 266 85 252 107 210
Polygon -7500403 true true 119 75 179 75 209 101 224 135 220 225 175 261 128 261 81 224 74 135 88 99

wheel
false
0
Circle -7500403 true true 3 3 294
Circle -16777216 true false 30 30 240
Line -7500403 true 150 285 150 15
Line -7500403 true 15 150 285 150
Circle -7500403 true true 120 120 60
Line -7500403 true 216 40 79 269
Line -7500403 true 40 84 269 221
Line -7500403 true 40 216 269 79
Line -7500403 true 84 40 221 269

wolf
false
0
Polygon -16777216 true false 253 133 245 131 245 133
Polygon -7500403 true true 2 194 13 197 30 191 38 193 38 205 20 226 20 257 27 265 38 266 40 260 31 253 31 230 60 206 68 198 75 209 66 228 65 243 82 261 84 268 100 267 103 261 77 239 79 231 100 207 98 196 119 201 143 202 160 195 166 210 172 213 173 238 167 251 160 248 154 265 169 264 178 247 186 240 198 260 200 271 217 271 219 262 207 258 195 230 192 198 210 184 227 164 242 144 259 145 284 151 277 141 293 140 299 134 297 127 273 119 270 105
Polygon -7500403 true true -1 195 14 180 36 166 40 153 53 140 82 131 134 133 159 126 188 115 227 108 236 102 238 98 268 86 269 92 281 87 269 103 269 113

x
false
0
Polygon -7500403 true true 270 75 225 30 30 225 75 270
Polygon -7500403 true true 30 75 75 30 270 225 225 270

@#$#@#$#@
NetLogo 5.2.0
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
default
0.0
-0.2 0 0.0 1.0
0.0 1 1.0 0.0
0.2 0 0.0 1.0
link direction
true
0
Line -7500403 true 150 150 90 180
Line -7500403 true 150 150 210 180

@#$#@#$#@
0
@#$#@#$#@