Sorry your browser is not supported!

   // calculate directional vector of 3D ray
   // using the view frustum
   // the result is a relative direction of the ray
   // to the camera
   vec4 rayDirection = normalize( vec4(
      tan(gl_TexCoord[0].x*fovy),
      tan(gl_TexCoord[0].y*fovy),
      1.0,
      1.0
   ));
   
   // calculate ray direction in global space
   // by adding view direction to ray direction
   // Y-rotation
   float angle = acos(eyeDirection.z/length(eyeDirection.xyz));
   float temp = rayDirection.x*cos(angle) - rayDirection.z*sin(angle);
   rayDirection.z = rayDirection.z*cos(angle) + rayDirection.x*sin(angle);
   // Z-rotation
   angle = atan(eyeDirection.y/eyeDirection.x);
   rayDirection.x = temp*cos(angle) - rayDirection.y*sin(angle);
   rayDirection.y = temp*sin(angle) + rayDirection.y*cos(angle);

   -------------------------------------------------------------
   Math
   -------------------------------------------------------------
   
   *** 3D Ray ***
   The direction of the ray relative to the camera
   can be determined by using the FOV of the camera frustum
                   /x*FOV\
   rayDirection = | y*FOV |
                  |   0   |
                   \  0  /
                   
   The absolute direction of the ray is calculated by adding
   the global angle of the camera to the angle of
   the relative direction of the ray (calculated above).
      The camera's direction is represented by a normalized directional
   vector, and refers to the global normalized Z-axis vector e3=(0,0,1).
   In order to rotate the rayDirection vector by the same amount,
   the angle between e3 and the view vector of the camera needs to be
   calculated:
   
             dot(e3,eyeDirection)     eyeDirection.z
   angle =  ---------------------- = ----------------
             |e3| * |eyeDirection|    |eyeDirection|
             
   A Y-rotation matrix can now be applied to the ray direction vector.
   
                      /cos(a)  0 sin(a)\
   newRayDirection = |   0     1   0    | * rayDirection
                      \-sin(a) 0 cos(a)/
                      
   The final position is then calculated by applying a Z-rotation
   matrix:
   
                        /cos(a) -sin(a) 0\
   finalRayDirection = | sin(a) cos(a)  0 | * newRayDirection
                        \  0      0     1/

txz={{u/Sqrt[u^2 + v^2], v/Sqrt[u^2 + v^2], 0},
 {-v/Sqrt[u^2 + v^2], u/Sqrt[u^2 + v^2], 0},
 {0, 0, 1}};

tz={{w/Sqrt[u^2 + v^2 + w^2], 0, -Sqrt[u^2 + v^2]/Sqrt[u^2 + v^2 + w^2]},
 {0, 1, 0},
 {Sqrt[u^2 + v^2]/Sqrt[u^2 + v^2 + w^2], 0, w/Sqrt[u^2 + v^2 + w^2]}
 };

r={{Cos[q],-Sin[q],0},{Sin[q],Cos[q],0},{0,0,1}};

txz={{u/Sqrt[u^2 + v^2], v/Sqrt[u^2 + v^2], 0},{-v/Sqrt[u^2 + v^2], u/Sqrt[u^2 + v^2], 0},{0, 0, 1}};
tz={{w/Sqrt[u^2 + v^2 + w^2], 0, -Sqrt[u^2 + v^2]/Sqrt[u^2 + v^2 + w^2]},{0, 1, 0},{Sqrt[u^2 + v^2]/Sqrt[u^2 + v^2 + w^2], 0, w/Sqrt[u^2 + v^2 + w^2]}};
r={{Cos[q],-Sin[q],0},{Sin[q],Cos[q],0},{0,0,1}};

finalMatrix=Inverse[txz].Inverse[tz].r.tz.txz//FullSimplify;

output: 

{{(u^2 + (v^2 + w^2) Cos[q])/(u^2 + v^2 + w^2), (
  u v - u v Cos[q] - w Sqrt[u^2 + v^2 + w^2] Sin[q])/(
  u^2 + v^2 + w^2), (
  u w - u w Cos[q] + v Sqrt[u^2 + v^2 + w^2] Sin[q])/(
  u^2 + v^2 + w^2)}, {(
  u v - u v Cos[q] + w Sqrt[u^2 + v^2 + w^2] Sin[q])/(
  u^2 + v^2 + w^2), (v^2 + (u^2 + w^2) Cos[q])/(u^2 + v^2 + w^2), (
  v w - v w Cos[q] - u Sqrt[u^2 + v^2 + w^2] Sin[q])/(
  u^2 + v^2 + w^2)}, {-((u w (-1 + Cos[q]))/(u^2 + v^2 + w^2)) - (
   v Sin[q])/Sqrt[
   u^2 + v^2 + w^2], -((v w (-1 + Cos[q]))/(u^2 + v^2 + w^2)) + (
   u Sin[q])/Sqrt[u^2 + v^2 + w^2], (w^2 + (u^2 + v^2) Cos[q])/(
  u^2 + v^2 + w^2)}}