/* Copyright (c) Mark J. Kilgard, 1998.  */

/* This program is freely distributable without licensing fees

   and is provided without guarantee or warrantee expressed or

   implied. This program is -not- in the public domain. */

#include <math.h>

#include "glsmapint.h"

/* (x,y,z) reflection vector --> (s,t) sphere map coordinates */

int

smapRvecToSt(float rvec[3], float st[2])

{

    double m, recipm;

    /* In Section 2.10.4 ("Generating texture coordinates")

       of the OpenGL 1.1 specification, you will find the

       GL_SPHERE_MAP equations:

          n' = normal after transformation to eye coordinates

          u  = unit vector from origin to vertex in eye coordinates

          (rx, ry, rz) = u - 2 * n' * transpose(n') * u

          m  = 2 * sqrt(rx^2 + ry^2 + (rz + 1)^2))

          s = rx/m + 0.5

          t = ry/m + 0.5

        The equation for calculating (rx, ry, rz) is the

        equation for calculating the reflection vector for

        a surface and observer.  The explanation and

        derivation for this equation is found in Roger's

        "Procedural Elements for Computer Graphics" 2nd ed.

        in Section 5-5 ("Determining the Reflection Vector").

        Note that Roger's convention has the Z axis in

        the opposite direction from the OpenGL convention. */

    m = 2 * sqrt(rvec[0]*rvec[0] +

                 rvec[1]*rvec[1] +

                 (rvec[2]+1)*(rvec[2]+1));

    if (m == 0.0) {

        /* Some point on the sphere map perimeter. */

        st[0] = 0.0;

        st[1] = 0.5;

        return 0;

    }

    recipm = 1.0/m;

    st[0] = rvec[0]*recipm + 0.5;

    st[1] = rvec[1]*recipm + 0.5;

    return 1;

}

/* (s,t) sphere map coordinate --> reflection verctor (x,y,z) */

void

smapStToRvec(float *st, float *rvec)

{

    double tmp1, tmp2;

    /* Using algebra to invert the sphere mapping equations

       shown above in smapRvecToSt, you get:

         rx = 2*sqrt(-4*s^2 + 4*s - 4*t^2 + 4*t - 1)*(2*s-1)

         ry = 2*sqrt(-4*s^2 + 4*s - 4*t^2 + 4*t - 1)*(2*t-1)

         rz = -8*s^2 + 8*s - 8*t^2 + 8*t - 3

       The C code below eliminates common subexpressions. */

    tmp1 = st[0]*(1-st[0]) + st[1]*(1-st[1]);

    tmp2 = 2 * sqrt(4*tmp1 - 1);

    rvec[0] = tmp2 * (2*st[0]-1);

    rvec[1] = tmp2 * (2*st[1]-1);

    rvec[2] = 8 * tmp1 - 3;

}