Path: utzoo!utgpu!utstat!jarvis.csri.toronto.edu!mailrus!nrl-cmf!ames!pasteur!miro.Berkeley.EDU!ph
From: ph@miro.Berkeley.EDU (Paul Heckbert)
Newsgroups: comp.graphics
Subject: Re: Ray tracing refraction
Summary: C code
Message-ID: <9389@pasteur.Berkeley.EDU>
Date: 2 Feb 89 04:54:49 GMT
References: <65@sdcc10.ucsd.EDU> <0XryqWy00Uo1875Ud-@andrew.cmu.edu> <25779@sgi.SGI.COM>
Sender: news@pasteur.Berkeley.EDU
Reply-To: ph@miro.Berkeley.EDU (Paul Heckbert)
Organization: University of California at Berkeley
Lines: 71

Here's some C code to compute the refracted ray direction.

Aside to Thant:
    Actually, Snell's law is n1*sin(theta1)=n2*sin(theta2);
    you were using the reciprocals of the indices of refraction.

Below is an excerpt of some notes I wrote for the Intro to Ray Tracing
SIGGRAPH tutorial.  (These notes are coming out soon as a book from Academic
Press, by Glassner, Arvo, Cook, Haines, Hanrahan, and Heckbert.
Included in the book is a derivation of the following formulas
from Snell's Law, which I would have included here except it's written in
eqn and troff and uses paste-up figures).

---------------------------------

Below is C code for SpecularDirection and TransmissionDirection,
routines which compute secondary ray directions.
The following formulas generate unit output vectors if given unit input vectors.

/*
 * SpecularDirection: compute specular direction R from incident direction
 * I and normal N.
 * All vectors unit.
 */
SpecularDirection(I, N, R)
Point I, N, R;
{
    VecAddS(-2.*VecDot(I, N), N, I, R);
}

/*
 * TransmissionDirection: compute transmission direction T from incident
 * direction I, normal N, going from medium with refractive index n1 to
 * medium with refractive index n2, with refraction governed
 * by Snell's law:  n1*sin(theta1) = n2*sin(theta2).
 * If there is total internal reflection, return 0, else set T and return 1.
 * All vectors unit.
 */
TransmissionDirection(n1, n2, I, N, T)
double n1, n2;
Point I, N, T;
{
    double eta, c1, cs2;

    eta = n1/n2;			/* relative index of refraction */
    c1 = -VecDot(I, N);			/* cos(theta1) */
    cs2 = 1.-eta*eta*(1.-c1*c1);	/* cos^2(theta2) */
    if (cs2<0.) return 0;		/* total internal reflection */
    VecComb(eta, I, eta*c1-sqrt(cs2), N, T);
    return 1;
}

where
    double VecDot(A, B)		dot product: returns A.B
    VecComb(a, A, b, B, C)	linear combination: C = aA+bB
    VecAddS(a, A, B, C)		add scalar multiple: C = aA+B

    typedef double Point[3];	/* xyz point data type */
    Point A, B, C;
    double a, b;

-------------

Whitted gives formulas for the refracted ray direction in his classic
paper on ray tracing: "An Improved Illumination Model for Shaded Display",
CACM, June 1980, but the formulas above compute faster.  It's a fun exercise
in trig and vector algebra to prove that Whitted's formulas are equivalent.

Paul Heckbert, CS grad student
508-7 Evans Hall, UC Berkeley		UUCP: ucbvax!miro.berkeley.edu!ph
Berkeley, CA 94720			ARPA: ph@miro.berkeley.edu