# gis_analysis.tcl
 #     Use NAP for analysing the impact of a plume of smoke on a
 #     populated area
 #
 package require Tk
 package require nap
 namespace import NAP::*

 # linearx --
 #     Construct the x-coordinate of a rectilinear grid
 #
 # Arguments:
 #     nx            Number of cells in x-direction
 #     ny            Number of cells in y-direction
 #     xmin          Minimum x-coordinate
 #     xmax          Maximum x-coordinate
 #
 # Result:
 #     2D array with x-coordinates
 #
 # Note:
 #     We need to deal with the reference count of the returned
 #     NAO (NAP data object) as the local variable that holds it
 #     will cease to exist when the procedure returns
 #
 proc linearx {nx ny xmin xmax} {
     nap "xstep = (xmax-xmin)/(nx-1)"

     nap "x = xmin + xstep * (0 .. nx-1)"
     nap "x = reshape(x, {[$nx] [$ny]})"
     $x set count 2 -keep

     return $x
 }

 # lineary --
 #     Construct the y-coordinate of a rectilinear grid
 #
 # Arguments:
 #     nx            Number of cells in x-direction
 #     ny            Number of cells in y-direction
 #     ymin          Minimum x-coordinate
 #     ymax          Maximum x-coordinate
 #
 # Result:
 #     2D array with x-coordinates
 #
 proc lineary {nx ny ymin ymax} {

     nap "y = transpose(linearx($nx,$ny,$ymin,$ymax))"
     $y set count 2 -keep

     return $y
 }

 # uniform --
 #     Construct a grid with a unform value
 #
 # Arguments:
 #     nx            Number of cells in x-direction
 #     ny            Number of cells in y-direction
 #     value         Value to be used
 #
 # Result:
 #     2D array with a uniform value
 #
 proc uniform {nx ny value} {

     nap "v = reshape(($nx*$ny)#$value, {[$nx] [$ny]})"
     $v set count 2 -keep

     return $v
 }

 # fillCircle --
 #     Fill cells in an existing grid that lie within a given circle
 #
 # Arguments:
 #     grid          Grid to be partially filled
 #     x             X-coordinate of the grid points
 #     y             Y-coordinate of the grid points
 #     xcentre       X-coordinate of the centre of the circle
 #     ycentre       Y-coordinate of the centre of the circle
 #     radius        Radius of the circle
 #     value         Value to be used
 #
 # Result:
 #     2D array with a filled-up area
 #
 proc fillCircle {grid x y xcentre ycentre radius value} {

     nap "v = (hypot(($x-$xcentre),($y-$ycentre)) < $radius)? $value : $grid"
     $v set count 2 -keep

     return $v
 }

 # gaussian --
 #     Create a matrix containing a Gaussian plume
 #
 # Arguments:
 #     x             X-coordinate of the grid points
 #     y             Y-coordinate of the grid points
 #     xorig         X-coordinate of the origin of the plume
 #     yorig         Y-coordinate of the origin of the plume
 #     disp          (Transversal) dispersion rate
 #     wind          Wind velocity
 #
 # Result:
 #     2D array with a plume in positive x-direction (unscaled)
 #
 proc gaussian {x y xorig yorig disp wind} {

 #         1       -y**2/(Dx/u)
 #       -----   e
 #  sqrt(2piDx/u)
 #
 #   x given in km, hence a factor 1000

     nap "pi = 4.0*atan(1.0)"
     nap "v = (1.0/sqrt(2.0*pi*disp*x/wind)) * exp( - y*y/(disp*x*1000.0/wind) )"
     $v set count 2 -keep

     return $v
 }

 # classes --
 #     Fill in values for each class
 #
 # Arguments:
 #     matrix        Matrix with class numbers
 #     values        Values to fill in per class
 #
 # Result:
 #     2D array with values filled in
 #
 proc classes {matrix values} {

     set classno 0
     nap "result = matrix"
     foreach value [$values value] {
         nap "result = (matrix == classno)? value : result"
         incr classno
     }
     $result set count 2 -keep

     return $result
 }

 #
 # Set up the map with the population
 #
 set nx 100
 set ny 100

 nap "x = linearx($nx, $ny,   0.0, 50.0)"
 nap "y = lineary($nx, $ny, -25.0, 50.0)"
 nap "population = uniform($nx, $ny, 1)"  ;# Very low density

 #
 # Fill in the urban and suburban areas
 #
 nap "population = fillCircle($population, $x, $y, 50.0, -5.0, 30.0, 2)"
 nap "population = fillCircle($population, $x, $y, 40.0, 10.0, 10.0, 3)"

 #
 # The plume of smoke
 #
 nap "conc = 40.0 * gaussian($x, $y, 0.0, 0.0, 10.0e-3, 5.0)"

 #
 # Which cells have a "high" concentration?
 #
 nap "irritating = conc > 0.2"

 #
 # Compute the number of people that will suffer:
 # - Mark the cells with a high concentration (simply: value of
 #   irritating times the type of population density)
 # - Count the cells per category (via the unary tally operator #)
 # - Multiply each category with the number of people per cell
 # - Sum all categories
 #
 nap "influenced = #reshape(irritating*population,{[expr {$nx*$ny}]})"

 set area_per_cell [expr {50.0*50.0/($nx*$ny)}]

 nap "casualties = sum( influenced * {0 1 10 50} * $area_per_cell )"
 nap "casualties_per_category = influenced * {0 1 10 50} * $area_per_cell"

 catch {
     console show
 }
 puts "Number of casualties: [$casualties value]"
 puts "Per category:"
 puts "Category\tNumber\tArea:"
 set values [lrange [$casualties_per_category value] 1 end]
 set influenced_area [lrange [$influenced value] 1 end]

 foreach c {Low Medium High} v $values a $influenced_area {
     puts "$c\t\t$v\t[expr {$a*$area_per_cell}]"
 }

 #
 # Now some pictures
 # Note:
 # We have not defined any coordinates for the 2D datasets ...
 #
 plot_nao casualties_per_category -type bar -title "Casualties per category"
 plot_nao [nap "conc>0.2? conc : _"] -title "Concentration"
 plot_nao [nap "classes(population,{0 1 10 50})"] -title "Population density"