MapRasterSeries

const float PI = 3.141592653589793;
const float EQUAL_EARTH_TOP = ${
  EQUAL_EARTH_TOP
};
const float NATURAL_EARTH_TOP = ${
  NATURAL_EARTH_TOP
};

// Equal Earth: the latitude of the row at raw y,
and the longitude a unit of raw x covers there
vec2 equalEarthRow(float y) {
  
const float A1 = 1.340264,
  A2 = -0.081106,
  A3 = 0.000893,
  A4 = 0.003796,
  M = 0.8660254037844386;
float l = y;
for (int i = 0; i < 8; i++) {
    
float l2 = l * l;
float l6 = l2 * l2 * l2;
l -= (l * (A1 + A2 * l2 + l6 * (A3 + A4 * l2)) - y) / (A1 + 3.0 * A2 * l2 + l6 * (7.0 * A3 + 9.0 * A4 * l2));

  }
float l2 = l * l;
float l6 = l2 * l2 * l2;
return vec2(asin(clamp(sin(l) / M,
  -1.0,
  1.0)),
  M * (A1 + 3.0 * A2 * l2 + l6 * (7.0 * A3 + 9.0 * A4 * l2)) / cos(l));

}

// Natural Earth,
likewise
vec2 naturalEarthRow(float y) {
  
float phi = y;
for (int i = 0; i < 12; i++) {
    
float p2 = phi * phi;
float p4 = p2 * p2;
phi -= (phi * (1.007226 + p2 * (0.015085 + p4 * (-0.044475 + 0.028874 * p2 - 0.005916 * p4))) - y) / (1.007226 + p2 * (0.045255 + p4 * (-0.311325 + 0.259866 * p2 - 0.065076 * p4)));

  }
float p2 = phi * phi;
return vec2(phi,
  1.0 / (0.8707 + p2 * (-0.131979 + p2 * (-0.013791 + p2 * p2 * p2 * (0.003971 - 0.001529 * p2)))));

}

void main() {
  
vec2 css = vec2(gl_FragCoord.x,
  u_height - gl_FragCoord.y) / u_resolution;
vec2 p = vec2(css.x - u_translate.x,
  u_translate.y - css.y) / u_scale;
vec3 v;
float alpha;
if (u_kind == 0) {
    
float r = length(p);
float edge = 0.5 / (u_scale * u_resolution);
alpha = clamp((1.0 + edge - r) / (2.0 * edge),
    0.0,
    1.0);
v = vec3(sqrt(max(0.0,
    1.0 - r * r)),
    p.x,
    p.y);

  }
else if (u_kind >= 3) {
    
// Rows off the map are clamped onto it,
    so the iteration stays finite
float limit = u_kind == 3 ? EQUAL_EARTH_TOP : NATURAL_EARTH_TOP;
float y = clamp(p.y,
    -limit,
    limit);
vec2 row = u_kind == 3 ? equalEarthRow(y) : naturalEarthRow(y);
float phi = clamp(row.x,
    -PI / 2.0,
    PI / 2.0);
float lambda = p.x * row.y;
// The curved sides are anti-aliased half a pixel either side
float edge = max(fwidth(lambda),
    1e-6);
alpha = abs(p.y) <= limit ? clamp((PI - abs(lambda)) / edge + 0.5,
    0.0,
    1.0) : 0.0;
lambda = clamp(lambda,
    -PI,
    PI);
v = vec3(cos(phi) * cos(lambda),
    cos(phi) * sin(lambda),
    sin(phi));

  }
else {
    
float phi = u_kind == 2 ? 2.0 * atan(exp(p.y)) - PI / 2.0 : p.y;
alpha = abs(p.x) <= PI && abs(phi) <= PI / 2.0 ? 1.0 : 0.0;
v = vec3(cos(phi) * cos(p.x),
    cos(phi) * sin(p.x),
    sin(phi));

  }
v = u_rotation * v;
vec2 uv = vec2(atan(v.y,
  v.x) / (2.0 * PI) + 0.5,
  0.5 - asin(clamp(v.z,
  -1.0,
  1.0)) / PI);

// Gradients without the jump where longitude wraps round
vec2 wrapped = vec2(fract(uv.x + 0.5),
  uv.y);
vec2 dx = dFdx(uv);
vec2 dy = dFdy(uv);
vec2 dxWrapped = dFdx(wrapped);
vec2 dyWrapped = dFdy(wrapped);
if (abs(dxWrapped.x) + abs(dyWrapped.x) < abs(dx.x) + abs(dy.x)) {
    
dx = dxWrapped;
dy = dyWrapped;

  }
vec4 day = textureGrad(u_image,
  uv,
  dx,
  dy);

// The sine of the sun's height there: day above the twilight band,
  night below
if (u_dayNight == 1) {
    
float t = smoothstep(-u_twilight,
    u_twilight,
    dot(normalize(v),
    u_sun));
day = mix(textureGrad(u_night,
    uv,
    dx,
    dy),
    day,
    t);

  }
color = day * alpha;