/** This program demonstrates the use of the country polygons in
    Country.frink to draw a map of the world.  It can be readily altered to
    draw the map in your favorite projection.  This demonstrates the
    Mollweide projection:

    https://mathworld.wolfram.com/MollweideProjection.html
    https://en.wikipedia.org/wiki/Mollweide_projection
*/

use Country.frink

/** This is a simple data struture that holds the data for a country that we
    read from a file. */
class CountryData
{
   var iso2
   var value
   var name

   new[iso, displayName, val] :=
   {
      iso2 = iso
      name = displayName
      value = val
   }
}

filename = "countries.txt"

// This is a dictionary where the key is the 2-letter ISO code and the
// value is a CountryData object.
countryDict = new dict

// Read the data from the file and make a dictionary out of it.
largestValue = 0
for line = lines[filenameToURL[filename]]
{
   // Split line on semicolons
   [iso2, name, val] = split[%r/;/, line]
   val = eval[val]   // Turn string into a number
   if isUnit[val]
   {
      countryDict@iso2 = new CountryData[iso2, name, val]

      // Keep track of largest value for scaling
      if val > largestValue
         largestValue = val
   }
}


g = new graphics
g.stroke[0.001]
g.font["SansSerif", "bold", 1 degree]

// Iterate through all countries' map data.
for [code, country] = Country.getCountryList[]
{
   cd = countryDict@code   // Fetch the CountryData from the dictionary
   if cd != undef
   {
      // Autoscale to largestValue
      red   = 1
      green = 1-(cd.value / largestValue)
      blue  = green
      cc = new color[red, green, blue]
   }
   
   first = true
   for poly = country.borders  // Iterate through polygons in a country.
   {
      p = new filledPolygon   // This polygon is the filled country
      po = new polygon        // This is the outline of the country
      for [long, lat] = poly  // Iterate through points in polygon
      {
         [x,y] = latLongToXYMollweide[lat degree, long degree]
         p.addPoint[x, -y]
         po.addPoint[x, -y]
      }

      // Draw filled countries
      g.color[cc]
      g.add[p]

      // The polygons in Country.frink are sorted so the first polygon is the
      // largest. Just label the largest.
      if first
      {
         [cx, cy] = p.getCentroid[]
         g.color[0,0,0]
         g.text[code, cx, cy]
      }

      // Draw country outlines
      g.color[0.2, 0.2, 0.2, 0.8]
      g.add[po]

      first = false
   }
}

// Optional:  Draw the ellipse that makes the boundary
// You could do this first but with fillEllipse... to draw water
//g.color[0,0,0,.3]
//g.drawEllipseSides[-2 sqrt[2], -sqrt[2], 2 sqrt[2], sqrt[2]]

g.show[]
g.write["heatmap.svg", 1000, undef]
g.write["heatmap.png", 1000, undef]
g.write["heatmap.html", 1000, 500]


/** Convert a latitude, longitude, and optional center longitude (long0)
    into x,y coordinates.

    The x coordinate ranges from -2 sqrt[2] to 2 sqrt[2]
    The y cooridnate ranges from -sqrt[2] to sqrt[2]]

    returns
     [x ,y]
*/
latLongToXYMollweide[lat, long, long0 = 0 degrees West] :=
{
   // The Mollweide projection uses an "auxiliary angle" theta where theta is
   // given by
   //   2 theta + sin(2 theta) = pi sin[lat]
   //
   // Where theta is found iteratively by iterating:
   // theta = theta-(2 theta + sin[2 theta] - pi sin[lat])/(2 + 2 cos[2 theta])
   //
   // or, by introducing a different supplementary angle theta1 where
   // theta = 1/2 theta1
   // The equations can be iterated as
   // theta1 = theta1 - (theta1 + sin[theta1] - pi sin[lat]) / (1 + cos[theta1])
   // with an initial guess
   // theta1 = 2 arcsin[2 lat / pi]
   theta1 = 2 arcsin[2 lat / pi]
   slat = sin[lat]

   do
   {
      ct = cos[theta1]
      if ct == -1
         break
      delta = - (theta1 + sin[theta1] - pi slat) / (1 + ct)
      theta1 = theta1 + delta
   } while abs[delta] > 1e-7  // error of 0.63 m on the earth
   
   theta = 1/2 theta1

   x = 2 sqrt[2]/pi (long - long0) cos[theta]
   y = sqrt[2] sin[theta]

   return [x,y]
}