# MollweideProjection.frink

``` /** This program demonstrates the use of the country polygons in     Countries.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 g = new graphics g.stroke[0.001] g.font["SansSerif", "bold", 1 degree] // Iterate through all countries. for [code, country] = Country.getCountryList[] {    cc = new color[randomFloat[0,1], randomFloat[0,1], randomFloat[0,1], .8]    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, -sqrt, 2 sqrt, sqrt] g.show[] g.write["Mollweide.svg", 1000, undef] g.write["Mollweide.png", 1000, undef] g.write["Mollweide.html", 1000, 500] /** Convert a latitude, longitude, and optional center longitude (long0)     into x,y coordinates.     The x coordinate ranges from -2 sqrt to 2 sqrt     The y cooridnate ranges from -sqrt to sqrt]     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/pi (long - long0) cos[theta]    y = sqrt sin[theta]    return [x,y] } ```