README.md

Compile with

gcc -o earthprojs earthprojs.c  -O6 -Wall -std=c99 -pedantic -Wextra `pkg-config --cflags --libs gdk-pixbuf-2.0` -lm

after making TEXTURE_DIR point to a directory containing subdirectories level$n each containing files tx_$i_$j.jpg of size 1024×1024 with $i ranging from 0 inclusive to 2^($n+1) exclusive and $j ranging from 0 inclusive to 2^$n exclusive. (I think this is a fairly "standard" earth texture file format, I don't remember where I got it from.) TEXTURE_LEVEL indicates which texture level to use.

Center of projection is defined by CENTER_LAT (in radians) at compile time and argv[1] (in full turns) at runtime. This is not at all logical, but that's what I found useful to generate the images in https://www.youtube.com/watch?v=LKcTbIsWS9c

PROJECTION indicates the projection to use: 0 is stereographic, 1 is gnomonic, 2 is orthographic, 3 is Lambert azimuthal equal-area, 4 is azimuthal equidistant, 5 is Mercator, 6 is cylindrical and 7 is trivial (plate-carrée).

earthprojs.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <math.h>

#include <glib.h>
#include <gdk-pixbuf/gdk-pixbuf.h>

#ifndef SIZE
#define SIZE 512
#endif

#ifndef IMGWIDTH
#define IMGWIDTH SIZE
#endif

#ifndef IMGHEIGHT
#define IMGHEIGHT SIZE
#endif

#ifndef PROJECTION
#define PROJECTION 0
#endif

#ifndef ANTIALIAS
#define ANTIALIAS 0
#endif

#ifndef CENTER_LAT
#define CENTER_LAT 0.409092814
#endif

#ifndef M_PI
#define M_PI   3.14159265358979323846  /* pi */
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923  /* pi/2 */
#endif
#ifndef M_PI_4
#define M_PI_4 0.78539816339744830962  /* pi/4 */
#endif

#ifndef TEXTURE_BASEDIR
#define TEXTURE_BASEDIR "."
#endif

#ifndef TEXTURE_LEVEL
#define TEXTURE_LEVEL 1
#endif

#define C_PROCESSOR_SUCKS(x) #x
#define C_PROCESSOR_SUCKS_2(x) C_PROCESSOR_SUCKS(x)

#ifndef TEXTURE_NAME_TEMPLATE
#define TEXTURE_NAME_TEMPLATE TEXTURE_BASEDIR "/level" C_PROCESSOR_SUCKS_2(TEXTURE_LEVEL) "/tx_%d_%d.jpg"
#endif

GdkPixbuf *texture_images[(1<<(TEXTURE_LEVEL+1))][(1<<TEXTURE_LEVEL)];

void
get_texture_int (int full_i, int full_j, int *red, int *green, int *blue)
{
  g_assert ( full_i >=0 && full_i<(1<<(11+TEXTURE_LEVEL)) && full_j>=0 && full_j<(1<<(10+TEXTURE_LEVEL)) );
  int image_i = full_i >> 10;
  int image_j = full_j >> 10;
  int pixel_i = full_i & 1023;
  int pixel_j = full_j & 1023;
  GdkPixbuf *pixbuf;
  if ( ! (texture_images[image_i][image_j]) )
    {
      char fnbuf[512];
      snprintf (fnbuf, sizeof(fnbuf),
		TEXTURE_NAME_TEMPLATE, image_i, image_j);
      GError *err = NULL;
      pixbuf = gdk_pixbuf_new_from_file (fnbuf, &err);
      if ( pixbuf == NULL )
	{
	  g_printerr ("%s\n", err->message);
	  exit (EXIT_FAILURE);
	}
      g_assert (gdk_pixbuf_get_colorspace (pixbuf) == GDK_COLORSPACE_RGB);
      g_assert (gdk_pixbuf_get_bits_per_sample (pixbuf) == 8);
      g_assert (gdk_pixbuf_get_n_channels (pixbuf) == 3);
      g_assert (gdk_pixbuf_get_width (pixbuf) == 1024);
      g_assert (gdk_pixbuf_get_height (pixbuf) == 1024);
      texture_images[image_i][image_j] = pixbuf;
    }
  else
    pixbuf = texture_images[image_i][image_j];
  int rowstride = gdk_pixbuf_get_rowstride (pixbuf);
  guchar *pixels = gdk_pixbuf_get_pixels (pixbuf);
  *red = pixels[pixel_j * rowstride + pixel_i * 3];
  *green = pixels[pixel_j * rowstride + pixel_i * 3 + 1];
  *blue = pixels[pixel_j * rowstride + pixel_i * 3 + 2];
}

void
get_texture (double pcx, double pcy, int *red, int *green, int *blue)
{
  g_assert ( pcx>=0 && pcx<2 && pcy>=0 && pcy<1 );
  int float_i = pcx*(1<<(10+TEXTURE_LEVEL));
  int float_j = pcy*(1<<(10+TEXTURE_LEVEL));
  int full_i = (int)(float_i);
  int full_j = (int)(float_j);
  double frac_i = float_i - full_i;
  double frac_j = float_j - full_j;
  int red00, green00, blue00;
  int red01, green01, blue01;
  int red10, green10, blue10;
  int red11, green11, blue11;
  get_texture_int (full_i, full_j, &red00, &green00, &blue00);
  get_texture_int (full_i, full_j+(full_j < (1<<(10+TEXTURE_LEVEL))-1 ? 1 : 0),
		   &red01, &green01, &blue01);
  get_texture_int ((full_i+1)%(1<<(11+TEXTURE_LEVEL)), full_j,
		   &red10, &green10, &blue10);
  get_texture_int ((full_i+1)%(1<<(11+TEXTURE_LEVEL)),
		   full_j+(full_j < (1<<(10+TEXTURE_LEVEL))-1 ? 1 : 0),
		   &red11, &green11, &blue11);
  *red = ((1-frac_i)*(1-frac_j)*red00
	  + (1-frac_i)*frac_j*red01
	  + frac_i*(1-frac_j)*red10
	  + frac_i*frac_j*red11 + 0.5);
  *green = ((1-frac_i)*(1-frac_j)*green00
	  + (1-frac_i)*frac_j*green01
	  + frac_i*(1-frac_j)*green10
	  + frac_i*frac_j*green11 + 0.5);
  *blue = ((1-frac_i)*(1-frac_j)*blue00
	  + (1-frac_i)*frac_j*blue01
	  + frac_i*(1-frac_j)*blue10
	  + frac_i*frac_j*blue11 + 0.5);
}

int
main (int argc, char *argv[])
{
  g_type_init ();
  double time = 0.;
  if ( argc >= 2 )
    sscanf (argv[1], "%lf", &time);
  time -= floor(time+0.5);
  printf ("P3\n%d %d\n255\n", IMGWIDTH, IMGHEIGHT);
  double sc = (SIZE-1.)/2;
  double scx = (IMGWIDTH-1.)/2;
  double scy = (IMGHEIGHT-1.)/2;
  for ( int i=0 ; i<IMGHEIGHT ; i++ )
    {
      for ( int j=0 ; j<IMGWIDTH ; j++ )
	{
	  int tred = 0, tgreen = 0, tblue = 0;
#if ANTIALIAS
	  for ( int smp = 0 ; smp < 9 ; smp++ )
	    {
#endif
	      int red, green, blue;
	      double ifl, jfl;
#if ANTIALIAS
	      ifl = i+((double)((smp%3)-1))/3;
	      jfl = j+((double)((smp/3)-1))/3;
#else
	      ifl = i;  jfl = j;
#endif
	      double px = (jfl-scx)/sc;
	      double py = (scy-ifl)/sc;
#if PROJECTION==7 // plate carrée
	      if ( fabs(py) >= M_PI_2 )
		{
		  red = 0;  green = 0;  blue = 0;
		  goto next;
		}
	      double mlat = py;
	      double ux = cos(mlat) * sin(px);
	      double uy = sin(mlat);
	      double uz = cos(mlat) * cos(px);
#elif PROJECTION==6 // cylindrical
	      if ( fabs(py) >= 1 )
		{
		  red = 0;  green = 0;  blue = 0;
		  goto next;
		}
	      double mlat = asin(py);
	      double ux = cos(mlat) * sin(px);
	      double uy = sin(mlat);
	      double uz = cos(mlat) * cos(px);
#elif PROJECTION==5 // Mercator
	      double mlat = 2*atan(exp(py)) - M_PI_2;
	      double ux = cos(mlat) * sin(px);
	      double uy = sin(mlat);
	      double uz = cos(mlat) * cos(px);
#elif PROJECTION==4 // Azimuthal equidistant
	      if ( px*px + py*py >= 1 )
		{
		  red = 0;  green = 0;  blue = 0;
		  goto next;
		}
	      double v = sqrt(px*px + py*py);
	      double ux = (px/v) * sin(M_PI*v);
	      double uy = (py/v) * sin(M_PI*v);
	      double uz = cos(M_PI*v);
#elif PROJECTION==3 // Lambert azimuthal equal-area
	      if ( px*px + py*py >= 1 )
		{
		  red = 0;  green = 0;  blue = 0;
		  goto next;
		}
	      double v = sqrt(1. - (px*px + py*py));
	      double ux = 2*px*v;
	      double uy = 2*py*v;
	      double uz = 1. - 2*(px*px + py*py);
#elif PROJECTION==2 // orthographic
	      if ( px*px + py*py >= 1 )
		{
		  red = 0;  green = 0;  blue = 0;
		  goto next;
		}
	      double ux = px;
	      double uy = py;
	      double uz = sqrt(1. - (px*px + py*py));
#elif PROJECTION==1 // gnomonic
	      double v = sqrt(1. + (px*px + py*py));
	      double ux = px/v;
	      double uy = py/v;
	      double uz = 1/v;
#elif PROJECTION==0 // stereographic
	      double v = 1. + (px*px + py*py);
	      double ux = 2*px/v;
	      double uy = 2*py/v;
	      double uz = (1. - (px*px + py*py))/v;
#else
#error Unknown projection.
#endif
	      g_assert (fabs(1 - (ux*ux + uy*uy + uz*uz)) < 1.e-6);
	      double ttx = ux;
	      double tty = cos(CENTER_LAT)*uy + sin(CENTER_LAT)*uz;
	      double ttz = -sin(CENTER_LAT)*uy + cos(CENTER_LAT)*uz;
	      double tx = cos(time*2*M_PI)*ttx - sin(time*2*M_PI)*ttz;
	      double ty = tty;
	      double tz = sin(time*2*M_PI)*ttx + cos(time*2*M_PI)*ttz;
	      double lng = atan2(tx, tz);
	      double lat = atan2(ty, hypot(tx, tz));
	      double pcx = lng/M_PI + 1;
	      double pcy = 0.5 - lat/M_PI;
	      get_texture (pcx, pcy, &red, &green, &blue);
	    next:
	      tred += red;  tgreen += green;  tblue += blue;
#if ANTIALIAS
	    }
	  tred = (tred+4)/9;  tgreen = (tgreen+4)/9;  tblue = (tblue+4)/9;
#endif
	  printf ("%3d %3d %3d  ", tred, tgreen, tblue);
	}
      printf ("\n");
    }
}