/*---------------------------------------------------------------------------
illum3d.asm:   3-D Graphics Illumination
-----------------------------------------------------------------------------
Description:

Performs 3-Dimensional Graphics Illumination using the Phong model for
ambient lighting plus diffuse and specular relection.  Phong's light model
is:

    I = IaKaOd  + FaIp[ Kd Od cos(A) + W(A) (cos B)**n ]

where:
    Ia   = Ambient light intensity
    Ka   = ambient reflection constant
    Od   = object's diffuse color
    Fa   = lightsource attenuation (square-law)
    Ip   = Intensity of light source
    Kd   = Diffuse reflection constant
    A    = angle between lightsource and normal vectors
    B    = angle between reflection  and view vectors
    W(A) = specular relection function, dependent on angle A
    n   = specular reflection exponent

Now we use following facts to simplify the equation:

1)  cos(A) = L*N, the dot-product of the lightsource and surface-normal 
        direction vectors.
2)  cos(B) = R*V, the dot-product of the reflection and view direction
        vectors.
3)  W(A) can be approximated by a constant, the specular reflection constant.
4)  KaOd and KdOd can be lumped into a single constant.
5)  Lightsource attenuation can be ignored.

Thus this equation can be simplified to:

    I = IaKa + Ip(Kd(L*N) + Ks(R*V)**n)

where:
    L   = direction vector from object to point light source
    N   = surface normal vector of object
    R   = direction vector of reflection
    V   = direction vector from object to viewpoint
    Ip  = Intensity of light source
    Kd  = lumped Diffuse reflection constant
    Ia  = Ambient light intensity
    Ka  = lumped ambient reflection constant
    Ks  = specular reflection constant
    n   = specular reflection exponent

To save the exponent calculation, a lookup table can be used.  This table
will be referred to as the specular reflection table.  Each element in the
specular reflection table is formed by the following calculation:

    elem[i] = Ks(i/(num_elems-1))**n    for 0 <= i <= num_elems

The dot-product of the normalized vectors R and V will produce a number
between 0.0 and 1.0.  This number is then scaled by num_elems to form an
index into the lookup table.

-----------------------------------------------------------------------------
Program Characteristics:

Calling Values:
    REGISTER FILES
    f15 = num_elems in spctbl (in floating-point)

    DAG1 (Data Memory)
    i0  = V Vector List (Vx,Vy,Vz...)
    l0  = 0
    m0  = +1
    i1  = N Vector List (Nx,Ny,Nz...)
    l1  = 0
    i2  = Ip, Kd
    l2  = 2

    DAG2 (Program Memory)
    i8  = R Vector List (Rx,Ry,Rz...)
    l8  = 0
    m8  = +1
    i9  = L Vector List (Rx,Ry,Rz...)
    l9  = 0
    i10 = start of Phong spectral table
    l10 = num_elems
    i11 = Ia, Ka 
    l11 = 2

Registers Used:
    r0, r4, r8-r15

Return Value:       f0 = Intensity I

Computation Time:   16 cycles, straight-line
                    14 cycles, looped

-----------------------------------------------------------------------------
Author:     Jim Donahue, Analog Devices DSP Division
References: Foley, VanDam, et.al. "Computer Graphics, Principles and
            Practice", 2nd Edition (Addison-Wesley 1990), pp.729-730
Revised:    4-APR-91
---------------------------------------------------------------------------*/

.GLOBAL illum3d;

.SEGMENT /pm    pm_code;
illum3d:    /* calculate dot products R*V, L*N */
                                f0=dm(i0,m0),  f4=pm(i8,m8);
        f8=f0*f4,               f0=dm(i0,m0),  f4=pm(i8,m8);
        f12=f0*f4,              f0=dm(i0,m0),  f4=pm(i8,m8);
        f12=f0*f4,  f8=f8+f12,  f0=dm(i1,m0),  f4=pm(i9,m8);
        f9=f0*f4,   f8=f8+f12,  f0=dm(i1,m0),  f4=pm(i9,m8);
        f12=f0*f4,              f0=dm(i1,m0),  f4=pm(i9,m8);
        f12=f0*f4,  f9=f9+f12,  f13=dm(i2,m0), f10=pm(i11,m8); /* f8=R*V  */
                    f9=f9+f12,  f14=dm(i2,m0), f11=pm(i11,m8); /* f9=L*N  */
        f0=f15*f8;                  /* generate spctbl index              */
        r0=fix f0;                  /* fix index                          */
        f7=f10*f11, m10=r0;         /* f7=IaKa, m10 load with next instr. */
                                    /*      use causes 1 cycle hold       */
        f9=f9*f14, f0=pm(m10,i10);  /* f9=Kd(L*N), f0=Ks(R*V)**n          */
        rts(db), f0=f0+f9;          /*   =Kd(L*N) + Ks(R*V)**n            */
            f0=f0*f13;              /*   =Ip(Kd(L*N) + Ks(R*V)**n)        */
            f0=f0+f7;               /*   =IaKA + Ip(Kd(L*N) + Ks(R*V)**n) */
.ENDSEG;