// Tex_sphere.cpp : Basic starting point for BHFS
//   a textured sphere viewed from inside that can warp texture at a pole.
//  Copyright Dewey Anderson 2001
//  You are free to use this code for educational purposes only.
 
#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glaux.h>
#include <math.h>   // for sin & cos
 

// Function prototypes
 
LONG WINAPI WndProc( HWND, UINT, WPARAM, LPARAM );
void DrawOpenGLScene( void );
HGLRC SetUpOpenGL( HWND hWnd );
 
// GR lensing th warping function
double warpth( double th, double warp );
 
// convert (phi_th) on the "black hole centered" sphere to (phi,th) on the texture image
// interpreted as a (phi,th) rectangle mapped to a sphere.  This is so that we can have the
// black hole located somewhere other than at the pole of the sky texture.
void bh_sphere_to_texture_sphere( double phi, double th, double *tex_phi, double *tex_th );
 
// convert (phi,th) coordinates on a textured sphere to (s,t) coordinates on texture
// image.  Possible correct for a kavraisky-like projection
void sphere_to_texture( double phi, double th, double *s, double *t );
 
 
 
// Globals
 
const double pi = 3.1415926535;
 
HDC hDC;    // Device context
 
// Animation parameters
double rot_angle;        // rotation angle of view
double delta_rot_angle;  // animation change per frame
double warp;             // animating black hole lensing warpage from pole
double delta_warp;       // animation change per frame
double warp_max;         // Upper limit to warp.
 
bool kavraisky_proj;  // true if texture image file is a kavraisky projection
                      // Otherwise we'll assume image file is rectangular (phi.th)
                      // This is so we can try different texture maps and see how they look
 
 
 
 
 
// Main program for Windows
// copied from OpenGl Programming for Windows book
//
int WINAPI WinMain(HINSTANCE hInstance,
                     HINSTANCE hPrevInstance,
                     LPSTR     lpCmdLine,
                     int       nCmdShow)
{
 static char szAppName[] = "OpenGL";
 static char szTitle[] = "BHFS Textured Celestial Sphere";
 WNDCLASS wc;
 MSG msg;
 HWND hWnd;
 
 wc.style = CS_HREDRAW | CS_VREDRAW;
 wc.lpfnWndProc = (WNDPROC)WndProc;
 wc.cbClsExtra = 0;
 wc.cbWndExtra = 0;
 wc.hInstance = hInstance;
 wc.hIcon = NULL;
 wc.hCursor = LoadCursor( NULL, IDC_ARROW);
 wc.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
 wc.lpszMenuName = NULL;
 wc.lpszClassName = szAppName;
 
 RegisterClass( &wc );
 
 hWnd = CreateWindow( szAppName, szTitle, WS_OVERLAPPEDWINDOW | WS_CLIPCHILDREN | WS_CLIPSIBLINGS,
                   CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL );
    if ( !hWnd )
  return 0;
 
 ShowWindow( hWnd, nCmdShow );
 UpdateWindow( hWnd );
 
 while ( GetMessage( &msg, NULL, 0, 0 )  )
 {
  TranslateMessage( &msg );
  DispatchMessage( &msg );
 }
 
 return( msg.wParam );
 
}   // WinMain()
 
 
 
// Window API message handling function.
// Mostly copied from OpenGl Programming for Windows book
//  I added code to WM_PAINT case to make it animate 
LONG WINAPI WndProc( HWND hWnd, UINT msg, WPARAM wParam, LPARAM lParam )
{
 HDC hDC;
 static HGLRC hRC;
 PAINTSTRUCT ps;
 GLsizei glnWidth, glnHeight;  // viewport dimensions
 GLdouble gldAspect;           // viewport aspect ratio
 
    
 switch (msg)
 {
   case WM_CREATE:
      hRC = SetUpOpenGL( hWnd );
  return 0;
 
  case WM_SIZE:
   hDC = GetDC( hWnd );
     wglMakeCurrent( hDC, hRC );
          // Extract width and height
       glnWidth = (GLsizei) LOWORD (lParam);
       glnHeight = (GLsizei) HIWORD (lParam);
       gldAspect = (GLdouble)glnWidth / (GLdouble)glnHeight;
          // Define projection matrix
       glMatrixMode( GL_PROJECTION );
       glLoadIdentity();
          gluPerspective( 45.0, gldAspect, 0.1, 1.2 );
       glViewport( 0, 0, glnWidth, glnHeight );
     wglMakeCurrent( NULL, NULL ); 
   ReleaseDC( hWnd, hDC );
  return 0;
 
  case WM_PAINT:
   hDC = BeginPaint( hWnd, &ps );
     wglMakeCurrent( hDC, hRC );
       DrawOpenGLScene();
     wglMakeCurrent( NULL, NULL );
   EndPaint( hWnd, &ps );
 
      // For animation, swap buffers
      SwapBuffers(hDC);
 
      // For animation, post a PAINT msg
      InvalidateRect(hWnd, 0, FALSE);
      PostMessage(GetParent(hWnd), WM_PAINT, 0, 0 );
   return 0;
 
  case WM_DESTROY:
   wglDeleteContext( hRC );
   PostQuitMessage( 0 );
  return 0;
 
 }  // msg switch
 
  return DefWindowProc( hWnd, msg, wParam, lParam );
 
}  // WndProc()
 

// Initialize OpenGL
// Pixel format & DC and RC code copied from OpenGl Programming for Windows book
HGLRC SetUpOpenGL( HWND hWnd )
{
 static PIXELFORMATDESCRIPTOR pfd = { sizeof (PIXELFORMATDESCRIPTOR),
                                   1, 
           PFD_DRAW_TO_WINDOW | PFD_SUPPORT_OPENGL | PFD_DOUBLEBUFFER,
           PFD_TYPE_RGBA,
           24,
           0, 0, 0,
           0, 0, 0,
           0, 0, 
           0, 0, 0, 0, 0, 
           16,
           0, 0,
           PFD_MAIN_PLANE,
           0, 0, 0, 0 };
 int nMyPixelFormatID;
 HGLRC hRC;
 
 hDC = GetDC( hWnd );
 nMyPixelFormatID = ChoosePixelFormat( hDC, &pfd );
 
 SetPixelFormat( hDC, nMyPixelFormatID, &pfd );
 
 hRC = wglCreateContext( hDC );
  // Make this context current so we can set up texture map
 wglMakeCurrent( hDC, hRC );
 

  // Read in the celestial sphere image for texture map
  // For development purposes, there are several different files to experiment with
  // Using different files is done by changing which line pair is not commented out
 
  _AUX_RGBImageRec *image;   // holds texture image read from file
 
//  image = auxDIBImageLoadA( "Virgo Cluster.dib" );
//  kavraisky_proj = true;  // not really but best viewed as one
 
//  image = auxDIBImageLoadA( "NGC 2264.dib" );
//  kavraisky_proj = true;  // not really but best viewed as one
 
//  image = auxDIBImageLoadA( "COBEsky.dib" );
//  kavraisky_proj = true;
 
  image = auxDIBImageLoadA( "milkyway_lund.dib" );
  kavraisky_proj = true;
 
//  image = auxDIBImageLoadA( "LivingEarth.dib" );
//  kavraisky_proj = false;
 
 
 
  // Might need to resize the universe image here for use as texture map
  // I made my image dimensions all powers of 2 with an external image program
  // so I don't need to resize.
 
  // Make a texture out of the universe image
  glPixelStorei( GL_UNPACK_ALIGNMENT, 1 );
 
  glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT );
  glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT );
  glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR );
  glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST );
 
  glTexImage2D( GL_TEXTURE_2D, 0, GL_RGB, image->sizeX, image->sizeY, 0, GL_RGB, GL_UNSIGNED_BYTE, image->data );
 
  // unmake the DC current and release it
 wglMakeCurrent( NULL, NULL ); ReleaseDC( hWnd, hDC );
 
  // Initialize animation parameters
  rot_angle = 0.;  // init to 0
  delta_rot_angle = 0;
  warp = 0.;
  delta_warp = 0.02;
  warp_max = 0.3;
 
 return hRC;
 
}  // SetUpOpenGL()
 

// Main function to actually draw the scene
// 
void DrawOpenGLScene()
{
 
 glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
  glEnable( GL_TEXTURE_2D );
  glTexEnvi( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL );
 

 glMatrixMode( GL_MODELVIEW );
  glLoadIdentity();
  glRotated( 180., 0.0, 1.0, 0.0 );  // so that default view is looking in the +Z direction
                                     // This is because I place the black hole in the th=0
                                     // direction which is at Z=+1.  This is behind the default
                                     // view so we turn around
  glRotated( rot_angle, 0.0, 1.0, 0.0 );  // animated view "panning around sphere"
 
  // This is a sphere made of quads bordered by lines of constant colatitude (th)
  // and constant longitude, phi.  The borders are th1, th2, phi1, phi2.
  // The quads are mapped to a texture map representing the celestial sphere.
  // To illustrate the gravitational lensing effect, we take the nominal th values
  // and "warp" them out by some function, wth = f(th).  f(th) should really represent
  // gravitational lensing.  For now we fake it.
 
  // th details
  int num_th = 16;               // number of colatitude bands
  double th1 = 0.0;              // starting colatitude
  double delta_th = pi/ num_th;  // loop step size.  theta goes from 0 to pi
 
  // phi details
  int num_phi = 32;                       // number of longitude "slices"
  double delta_phi = 2 * pi / num_phi;    // loop step size.  phi goes from 0 to 2*pi
 
  // Make the actual quads comprising this sphere
 glBegin( GL_QUADS );
 
    // Outer loop is each colatitude, th
    for ( int jj=0; jj<num_th; jj++ )
    {
      double th2 = th1 + delta_th;  // larger colatitude edge of band
 
      // compute the warped angular positions of th1 and th2
      // REPLACE WITH GR CALCULATION
      double wth1 = warpth(th1, warp);
      double wth2 = warpth(th2, warp);
 
      // precalculate some terms used inside the inner loop
      double r1 = sin(wth1);     // r is the radial distance from axis to surface
      double r2 = sin(wth2);
      double coswth1 = cos(wth1);
      double coswth2 = cos(wth2);
 
      double phi1 = 0;                        
 
      // inner loop for each phi slice in a band
      for ( int ii=0; ii<num_phi; ii++ )
      {
        double phi2 = phi1 + delta_phi;    // ccw border of slice
 
        // convert the (phi,th) coordinates on the black hole pole sphere to
        // (phi,th) coordinates on the textured celestial sphere
        // This enables placing the black hole somewhere other than at the top of the texture
        // Polygon side may no longer be "coordinate aligned" so we have to go from just
        // knowing 2 values each for phi and th and knowing 4 separate pairs
        // Polygon has vertices a, b, c & d
        double phia, tha, phib, thb, phic, thc, phid, thd;
        bh_sphere_to_texture_sphere( phi1, th1, &phia, &tha );
        bh_sphere_to_texture_sphere( phi1, th2, &phib, &thb );
        bh_sphere_to_texture_sphere( phi2, th2, &phic, &thc );
        bh_sphere_to_texture_sphere( phi2, th1, &phid, &thd );
 
        // get (s,t) texture coordinates for each point
        double sa, ta, sb, tb, sc, tc, sd, td;
        sphere_to_texture( phia, tha, &sa, &ta );
        sphere_to_texture( phib, thb, &sb, &tb );
        sphere_to_texture( phic, thc, &sc, &tc );
        sphere_to_texture( phid, thd, &sd, &td );
 
        // just to avoid duplicate computations
        double cos1 = cos(phi1);  
        double sin1 = sin(phi1);
        double cos2 = cos(phi2);
        double sin2 = sin(phi2);
 
        // Now define the actual vertices along with their texture coords
        glTexCoord2d( sa, ta );
        glVertex3d( r1*cos1, r1*sin1, coswth1 );
 
        glTexCoord2d( sb, tb );
        glVertex3d( r2*cos1, r2*sin1, coswth2);
 
        glTexCoord2d( sc, tc );
        glVertex3d( r2*cos2, r2*sin2, coswth2);
 
        glTexCoord2d( sd, td );
        glVertex3d( r1*cos2, r1*sin2, coswth1);
 
        // increment to next phi
        phi1 = phi2;           
 
      }  // for each phi slice
 
      // step out to next colatitude
      th1 = th2;           
 
    }  // for each ring
 
  glEnd();
 
  glFlush();    // draw everything in the queue
 

  // update animated parameters.
  // The basic animation is to stare at the pole as it warps larger like we're
  // getting closer to a black hole.  Then stop and pan around the sky behind and 
  // back to the black hole.  Then "back away" as warp shrinks.  Repeat.
 
  // increment the animation parameters
  rot_angle += delta_rot_angle;
  warp += delta_warp;
 
  // See if we've reached a point where we need to change deltas for different motion
 
  // Is warp big enough that we should start pan?
  if ( warp > warp_max )
  {                         // when warp is max
    delta_warp = 0.0;       // stop changing warp
    delta_rot_angle = 3.0;  // start rotating 
  }
 
  // Have we panned all the way around?
  if ( rot_angle >= 360 )
  {                         // when rotated around
    delta_rot_angle = 0.0;  // stop rotating, 
    delta_warp = -0.02;    // reverse warp to shrinking
    rot_angle = 0.0;       // reset just to keep things repeatable
  }
 
  // Has warp shrunk all the way down?
  if ( warp < 0 )
  {
    delta_warp = 0.02;  // reverse to growing
    warp = 0.;          // reset just to keep things repeatable
  }
 
 
 
}  //  DrawOpenGLScene()
 
 
 
// Relativistic gravitational lense warping, "warp theta"
//  This is just a stand-in and not the real GR warp function
//    hyperbola asymptotic to wth=th, then normalized so what warpth(pi) = pi
double
warpth( double th, double warp )
{
  double warpedth = sqrt(warp*warp + th*th);   
  double warpedpi = sqrt(warp*warp + pi*pi);
 
  double normalized_warp = warpedth * pi/warpedpi ;  // mormalize so that w(pi) = pi
 
  return normalized_warp;
 
}   // warpth
 

// convert spherical polar coordinates on the sphere with the "black hole pole"
// to spherical polar coordinates on a textured celestial sphere.
// Imagine the polygon sphere as a clear plastic sphere over a textured celestial sphere
// The black hole is located somewhere on that undistorted celestial sphere
// Rotate the clear sphere around to make its north pole "point toward" (overlap) the
// black hole's position.  Convert (phi,th) coordinate positions from the clear sphere
// to (tex_phi, tex_th) coordinate positions on the textured sphere.
// Back in the main program, you can imagine the textured sphere being "unwrapped" to a 
// rectangle to get the (s,t) coordinates on the texture. 
void
bh_sphere_to_texture_sphere( double phi, double th, double *tex_phi, double *tex_th )
{
  // convert to Cartesian coordinates
  double r = sin(th);
  double x = r*cos(phi);
  double y = r*sin(phi);
  double z = cos(th);
 
  // rotate the point alpha degrees around the X axis.
  // This should be replaced with a more general "pointing" algorithm that lets you 
  // "point" anywhere.
  double alpha = 1.4;  // a little less than pi/2
  double xp = x;
  double yp = y*cos(alpha) - z*sin(alpha);
  double zp = y*sin(alpha) + z*cos(alpha);
 
  // convert back to spherical coords
  *tex_th = acos(zp);
  *tex_phi = ( yp >= 0 )  ?  acos(xp/sin(*tex_th)) : 2*pi - acos(xp/sin(*tex_th));
 
}   // bh_sphere_to_texture_sphere()
 

// convert (phi,th) coordinates on a textured sphere to (s,t) coordinates on texture
// image.  Possible correct for a kavraisky-like projection
void
sphere_to_texture( double phi, double th, double *s, double *t )
{
  // WARNING: This algorithm doesn't handle polygons that span over the edge correctly
 
  // scale th, 0->pi to be 0->1.
  // th runs from 0 at the "pole" to pi/2 at the equator to pi at the bottom
  // It seems more natural to think of th=0 as being at the top of the texture, so invert it
  *t = 1.0 - th / pi;
 
  // Actual s coordinates depend on whether texture image is an kavraisky-like projection
  // or just an orthogonal phi,th image
  if ( kavraisky_proj )
  {
    // Assume a mapping of the sphere onto an ellipse in the texture rectangle
    // Since the rectanglular texture is not just a (phi,th) map, we compress s based on t
    // I don't know that this is the exact undoing of the kavraisky-like projection
    // used to generate the actual image from sky data.
    *s = ( phi / (2*pi) - 0.5 ) * sin(th) + 0.5;
  }
 
  else  // Texture is just a rectangular ph,th image
    *s = phi / (2*pi) ;
 
}   // sphere_to_texture