Solar Oscillations Demo

One of my favourite programs from Sky and Telescope magazine's Astronomical Computing column was one which produced a visualisation of oscillations on the surface of the Sun. Although these oscillations were discovered in 1960 they gained more attention during the 1990s due to the very successful MDI instrument on board the SOHO spacecraft.

The program uses Legendre polynomials to illustrate the red- and blueshift of light due to pressure waves at the surface. The user specifies the degree L (typically up to 150) and order M (not exceeding L) of the oscillation mode. In reality the Sun contains many such wave modes superimposed.

Below is a mostly faithful port to C/C++ using bitmap output to ensure cross-platform compatibility. I'm using GCC version 5.3 and the code should compile on almost any PC or ARM-based computer (for example Raspberry Pi). The bitmap writing function assumes a little-endian CPU.

/*
 Solar Surface Oscillations
 SOLAROSC.BAS by John Kennewell
 */
#include <cmath>
#include <iostream>
#include <cstdio>

void writeBitmapFile(int x, int y, unsigned int *data);

using namespace std;

class LegendrePoly
{
public:
 static const double HALF_PI = 3.1415926535897936 / 2.0;
 static const int SAMPLES = 400;

 LegendrePoly(int _l, int _m)
  : l(_l), m(_m)
 {}

 ~LegendrePoly() {}

 double evaluate(const double cost)
 {
  int idx = (cost + 1.0) * SAMPLES / 2;
  return p[idx];
 }
 // Generate Legendre function
 void generate()
 {
  double max_val = 0;
  for(int xi = 0; xi < SAMPLES; ++xi)
  {
   double x = ((double)xi - SAMPLES/2) / (SAMPLES/2);
   double d = sqrt(1 - x*x);
   double p1 = exp(l * log(d)), p2 = 0; // p1 = (1 - x^2)^(l/2)

   for(int li = 1; li <= (2*l - 1); li += 2)
    p1 *= li;

   if ((m >= l) || p1 == 0) continue;
   double pm;
   for (int mi = l-1; mi>=m; --mi)
   {
    pm = 2 * (mi+1) * x * p1 / d - p2;
    pm /= ((l-mi) * (l+mi + 1));
    p2 = p1; p1 = pm;
   }
   p[xi] = p1;
   if(abs(p1) > max_val)
    max_val = abs(p1);
  }
  // normalise
  for(int xi = 0; xi < SAMPLES; ++xi)
   p[xi] /= max_val;
 }

private:
 const int l, m; // Degree and order of Legendre function
 double p[SAMPLES]; // polynomial coeffs
};


unsigned int colourLookup(double s)
{
 // rc = rd * INT((15 - i) * 4.5)
 // bc = bl * INT((i - 1) * 4.5)
 // gc = gn * INT((7 - ABS(8 - i)) * 4.5)
 unsigned int red =   (unsigned int)((1-s) * 255);
 unsigned int blue =  (unsigned int)(s * 255);
 unsigned int green =  (unsigned int)((0.6 - abs(s-0.5)) * 255);
 return (blue | (green << 8) | (red << 16));
}

int main(int argc, char **argv)
{
 int L, M;
 unsigned int data[400*400];

 cout << "Enter degree L and order M: ";
 cin >> L >> M;
 cout << "Solar Global Oscillations, Mode L =" << L << ", M =" << M << endl;

 LegendrePoly poly(L,M);
 poly.generate();

 // Plot solar hemisphere
 double s;
 for(int y=0; y<400; ++y)
 {
  double cost = (y-200.0) / 200;
  double sint = sqrt(1 - cost*cost);
  double pm = poly.evaluate(cost);
  for(int x=0; x<400; ++x)
  {
   double sinf = (x-200.0) / 200 / sint;
   if ((sinf > -1) || (sinf < 1))
   {
    double cosf = sqrt(1 - sinf*sinf);
    double phi = LegendrePoly::HALF_PI * sinf / abs(sinf);
    if (cosf != 0)
     phi = atan2(sinf, cosf);
    s = pm * cos(M * phi) * sint * cosf;
    data[400*y+x] = colourLookup(s/2+0.5);
   }
  }
 }
 writeBitmapFile(400, 400, data);
}


void writeBitmapFile(int x, int y, unsigned int *data)
{
 FILE *fp;
 fp = fopen("solar_osc.bmp", "w");

 /* BMP header */
 fwrite("BM", 2, 1, fp);
 int filesize = x*y*4 + 14 + 40; // 4 bytes per pixel
 fwrite(&filesize, sizeof(filesize), 1, fp);
 // reserved 
 fwrite("SOLO", 4, 1, fp);
 // offset to image data
 int offset = 40;
 fwrite(&offset, sizeof(offset), 1, fp);

 /* DIB header: BITMAPINFOHEADER */
 int hdrsize = 40; // header size in Bytes
 fwrite(&hdrsize, sizeof(hdrsize), 1, fp);
 fwrite(&x, 4, 1, fp);
 fwrite(&y, 4, 1, fp);

 short num_planes = 1, color_depth = 32;
 fwrite(&num_planes, 2, 1, fp);
 fwrite(&color_depth, 2, 1, fp);

 int comp_format = 0;
 fwrite(&comp_format, 4, 1, fp);
 int size = x*y*4; // 4 bytes per pixel
 fwrite(&size, sizeof(size), 1, fp);

 fwrite(&x, 4, 1, fp);
 fwrite(&y, 4, 1, fp);
 char padding[] = {0,0,0,0};
 fwrite(padding, 4, 2, fp);

 // image data
 fwrite(data, size, 1, fp);
 fclose(fp);
}

Links

• BASIC programs from Sky and Telescope magazine.
• Lecture course at UCL/MSSL.