Once you translate the point of intersection to the origin you then have a normal vector and a direction vector of the incoming light ray. If you reverse the direction vector, normalize it and calculate <n,d> you will get cos(theta).
Now you need to calculate a new direction vector that is going out from the surface. To do this you must reflect the vector about the normal vector.
This is because the reversed direction vector (from the incoming vector of light) and the final direction vector are symmetrical with respect to the normal.
The rotated vector needs to be rotated 2*theta along the plane formed by the normal vector and the reversed direction vector. The plane has normal r x n (r is reversed incoming vector, n is normal) which is equal to n x i (incoming vector) where all vectors are unit length.
Then you rotate the direction vector about this axis of 2*theta and you're done. Also you don't need to calculate theta directly: you can use cos(theta) calculated from your inner product (and you can calculate sin(theta) from cross product) and substitute this into the rotation matrices directly.
So in summary:
1) Translate point of contact to the origin (subtract the value of the contact point)
2) Form vectors for normal of surface and reversed direction of incoming light
3) Calculate cos(theta) and sin(theta) using dot and cross products
4) Form rotation matrix that rotates reverse direction vector about normal formed by normal of surface and reversed direction vector
5) Obtain final direction vector
6) Translate this vector back to original co-ordinate system (add the contact point obtained from 1) and repeat process for further rays.