Thread Subject:

How to rotate a (3-D) line in a plane by an angle of theta about the point of origin ?

Hello All,
I have rays traveling in different directions, hitting spherical bubbles and moving on. However, upon hitting the bubbles, the direction of rays changes by an angle theta with respect to the normal at the point of intersection (ray with bubble). Despite the direction change, the new ray should be in the same plane as the original ray which as well passes through the center of the bubble.
I have the equation of the original ray as I have the point (P1) and direction of emission, I calculate the point of intersection (P2) with the sphere. I know the center of the bubble (Pc), hence I have the equation of the plane as well . I also write an equation for the normal to the sphere at the point of intersection (Pc-P2), which I would like to rotate with an angle theta that I calculate elsewhere about point (P2).
I do not know how to ask Matlab to calculate the equation of the new ray? Specially since this should happen in a loop, I can not do the calculation on pen and paper for once on pen and paper and be over with it
On pen and paper, I make use of orthogonality of the normal to the plane of rays and the new ray and the angle theta (dot product) between the normal and the new ray and guessing a value for x for an unknown point, I calculate one point on the new line (ray).
There should be a way to do it simpler and I was wondering if anyone knows something that could help me.
Thanks in advance.

Regards,
Mahsa.

"Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> Hello All,
> I have rays traveling in different directions, hitting spherical bubbles and moving on. However, upon hitting the bubbles, the direction of rays changes by an angle theta with respect to the normal at the point of intersection (ray with bubble). Despite the direction change, the new ray should be in the same plane as the original ray which as well passes through the center of the bubble.
> I have the equation of the original ray as I have the point (P1) and direction of emission, I calculate the point of intersection (P2) with the sphere. I know the center of the bubble (Pc), hence I have the equation of the plane as well . I also write an equation for the normal to the sphere at the point of intersection (Pc-P2), which I
would like to rotate with an angle theta that I calculate elsewhere about point (P2).
============================

This might help.

http://www.mathworks.com/matlabcentral/fileexchange/30864-3d-rotation-about-shifted-axis

You just need to know the normal vector to the plane of rotation, which I can't quite glean from your description. It might be given by cross(P1-P2,Pc-P2)

"Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> Hello All,
> I have rays traveling in different directions, hitting spherical bubbles and moving on. However, upon hitting the bubbles, the direction of rays changes by an angle theta with respect to the normal at the point of intersection (ray with bubble). Despite the direction change, the new ray should be in the same plane as the original ray which as well passes through the center of the bubble.
> I have the equation of the original ray as I have the point (P1) and direction of emission, I calculate the point of intersection (P2) with the sphere. I know the center of the bubble (Pc), hence I have the equation of the plane as well . I also write an equation for the normal to the sphere at the point of intersection (Pc-P2), which I would like to rotate with an angle theta that I calculate elsewhere about point (P2).
> I do not know how to ask Matlab to calculate the equation of the new ray? Specially since this should happen in a loop, I can not do the calculation on pen and paper for once on pen and paper and be over with it
> On pen and paper, I make use of orthogonality of the normal to the plane of rays and the new ray and the angle theta (dot product) between the normal and the new ray and guessing a value for x for an unknown point, I calculate one point on the new line (ray).
> There should be a way to do it simpler and I was wondering if anyone knows something that could help me.
> Thanks in advance.
>
> Regards,
> Mahsa.
- - - - - - - - -
  Are you talking about a reflection of the ray off the surface of your bubble, or a refraction in which the ray enters the bubble at an altered angle determined by the refractive index of the bubble's interior?

  Assuming this is a simple reflection and therefore your incoming and outgoing rays make equal angles with respect to the bubble's surface normal and lie on opposite sides of it, you can compute the outgoing direction as follows. Let P be a vector pointing along the direction of the incoming ray and let N be a vector pointing along the surface normal direction. Then the direction of the reflected ray will be along the vector:

 Q = P - 2*dot(P,N)/dot(N,N)*N;

In other words the component of P along a direction orthogonal to the normal is preserved while the component along the normal itself is reversed.

  If we are talking about a refraction, let P again be pointing in the direction of the incoming ray and let N be a surface normal which points inward toward the bubble's interior. If I interpret your statement correctly, you want the direction of the ray to bend through an angle theta toward the normal. Is that correct? If so, you can compute the refracted ray direction along vector Q as:

 Q = cross(P,N);
 Q = cos(theta)*P + sin(theta)*cross(Q/norm(Q),P);

  Of course these Q's are only vectors pointing along the altered rays. To establish the equations of their lines, you must make the lines pass through the point of surface intersection.

Roger Stafford

> "Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> > ...... the direction of rays changes by an angle theta with respect to the normal ....
- - - - - - - - -
  By the way, if your question was about refraction, when you said, "the direction of rays changes by an angle theta with respect to the normal", by Snell's Law of refraction that angle theta would depend on the angle of incidence of the ray, and therefore you would need to calculate that angle from the direction of the incoming ray with respect to the surface normal.

Roger Stafford

"Roger Stafford" wrote in message <ir8rqu$hkl$1@newscl01ah.mathworks.com>...
> > "Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> > > ...... the direction of rays changes by an angle theta with respect to the normal ....
> - - - - - - - - -
> By the way, if your question was about refraction, when you said, "the direction of rays changes by an angle theta with respect to the normal", by Snell's Law of refraction that angle theta would depend on the angle of incidence of the ray, and therefore you would need to calculate that angle from the direction of the incoming ray with respect to the surface normal.
>
> Roger Stafford

Dear Mr. Stafford,
Thanks a lot for the advice. I was indeed talking about simple reflection and refraction in ray tracing. I have bubbles in a liquid containing particles. The very tiny particles have a phase function indicating the new direction of flight which I have implemented. Therefore, I have the direction of incidence, I try to solve the equation of line(ray) with sphere to calculate the point of incidence. Then I can have the equation of normal to the surface of sphere that passes through the center of bubble and the incidence point. Then simple with a dot product I calculate the angle of incidence(theta1). It is a Monte Carlo scheme, so the ray is either reflected or refracted, and I need to define the new direction of flight and this direction (as you know) must be in the same plane as the incidence ray and the normal to the sphere surface. With Snell's law I calculate the angle of refraction
with respect to normal (theta2). Then basically the new ray will make theta2 with normal (depending on how we have defined the direction of normal, inward or outward). So, you understood correctly what I was talking bout, but I don't want the ray to bent toward the normal, I simply want theta2 to be the angle they make with each other. May I kindly ask you to elaborate the second equation for Q in the refraction case.
Thanks a million.
Regards,
Mahsa.



Thanks again. I appreciate

"Matt J" wrote in message <ir67mj$he2$1@newscl01ah.mathworks.com>...
> "Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> > Hello All,
> > I have rays traveling in different directions, hitting spherical bubbles and moving on. However, upon hitting the bubbles, the direction of rays changes by an angle theta with respect to the normal at the point of intersection (ray with bubble). Despite the direction change, the new ray should be in the same plane as the original ray which as well passes through the center of the bubble.
> > I have the equation of the original ray as I have the point (P1) and direction of emission, I calculate the point of intersection (P2) with the sphere. I know the center of the bubble (Pc), hence I have the equation of the plane as well . I also write an equation for the normal to the sphere at the point of intersection (Pc-P2), which I
> would like to rotate with an angle theta that I calculate elsewhere about point (P2).
> ============================
>
> This might help.
>
> http://www.mathworks.com/matlabcentral/fileexchange/30864-3d-rotation-about-shifted-axis

Thanks a lot. I am looking at the link. The normal is obtained as you mentioned.
>
> You just need to know the normal vector to the plane of rotation, which I can't quite glean from your description. It might be given by cross(P1-P2,Pc-P2)

"Mahsa " <newsreader@mathworks.com> wrote in message <ir8vbe$ppj$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ir8rqu$hkl$1@newscl01ah.mathworks.com>...
> > > "Mahsa " <newsreader@mathworks.com> wrote in message <ir62ul$31b$1@newscl01ah.mathworks.com>...
> > > > ...... the direction of rays changes by an angle theta with respect to the normal ....
> > - - - - - - - - -
> > By the way, if your question was about refraction, when you said, "the direction of rays changes by an angle theta with respect to the normal", by Snell's Law of refraction that angle theta would depend on the angle of incidence of the ray, and therefore you would need to calculate that angle from the direction of the incoming ray with respect to the surface normal.
> >
> > Roger Stafford
>
> Dear Mr. Stafford,
> Thanks a lot for the advice. I was indeed talking about simple reflection and refraction in ray tracing. I have bubbles in a liquid containing particles. The very tiny particles have a phase function indicating the new direction of flight which I have implemented. Therefore, I have the direction of incidence, I try to solve the equation of line(ray) with sphere to calculate the point of incidence. Then I can have the equation of normal to the surface of sphere that passes through the center of bubble and the incidence point. Then simple with a dot product I calculate the angle of incidence(theta1). It is a Monte Carlo scheme, so the ray is either reflected or refracted, and I need to define the new direction of flight and this direction (as you know) must be in the same plane as the incidence ray and the normal to the sphere surface. With Snell's law I calculate the angle of refraction

> ....... Then basically the new ray will make theta2 with normal .......
- - - - - - - - - - -
  Yes I misunderstood what angle you meant by theta. As I understand it now, if we again regard the normal N as pointing inward from the surface, the angle theta is to be measured from this normal toward, and in the same plane as, the incident ray P. In that case the refracted vector Q would be:

 Q = cross(P,N);
 Q = (cos(theta)*N + sin(theta)*cross(N,Q/norm(Q)))/norm(N);

which is a very different formula.

Roger Stafford

"Roger Stafford" wrote in message <ir95nq$be6$1@newscl01ah.mathworks.com>...
> Yes I misunderstood what angle you meant by theta. As I understand it now, if we again regard the normal N as pointing inward from the surface, the angle theta is to be measured from this normal toward, and in the same plane as, the incident ray P. In that case the refracted vector Q would be:
>
> Q = cross(P,N);
> Q = (cos(theta)*N + sin(theta)*cross(N,Q/norm(Q)))/norm(N);
>
> which is a very different formula.
>
> Roger Stafford
- - - - - - - -
  In case you wish to bypass the calculation of theta, you can compute Q directly. Let r be the ratio of the outer index of refraction to the inner index, let P be the vector of the incoming ray, and N be the inward pointing normal direction. Then do this:

 t1 = r*dot(P,N);
 t2 = dot(P,P)*(1-r^2);
 Q = r*P + t2/(sqrt(t1^2+t2*dot(N,N))+t1)*N;

  If the quantity within the square root is negative, you have an r greater than one, (that is the velocity increases,) and an impossible combination of the P and N vectors.

Roger Stafford

"Roger Stafford" wrote in message <ir9mcc$kmm$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <ir95nq$be6$1@newscl01ah.mathworks.com>...
> > Yes I misunderstood what angle you meant by theta. As I understand it now, if we again regard the normal N as pointing inward from the surface, the angle theta is to be measured from this normal toward, and in the same plane as, the incident ray P. In that case the refracted vector Q would be:
> >
> > Q = cross(P,N);
> > Q = (cos(theta)*N + sin(theta)*cross(N,Q/norm(Q)))/norm(N);
> >
> > which is a very different formula.
> >
> > Roger Stafford
> - - - - - - - -
> In case you wish to bypass the calculation of theta, you can compute Q directly. Let r be the ratio of the outer index of refraction to the inner index, let P be the vector of the incoming ray, and N be the inward pointing normal direction. Then do this:
>
> t1 = r*dot(P,N);
> t2 = dot(P,P)*(1-r^2);
> Q = r*P + t2/(sqrt(t1^2+t2*dot(N,N))+t1)*N;
>
> If the quantity within the square root is negative, you have an r greater than one, (that is the velocity increases,) and an impossible combination of the P and N vectors.
>
> Roger Stafford

Thanks for the tip on the vector notation for the ray and normal. I was a bit confused because I was working with points and line equations rather than vector notation. It is much easier this way, thanks.
It was an interesting suggestion to skip the calculation of theta, however, the formulation is not straight forward for me. I do have increasing speed in my simulation (water+solid to air bubble), so my r is greater than one. If we assume that the vectors are normalized, could we define the new vector without having to separately calculate theta of refraction (although it is not a big deal)? thanks.

"Mahsa " <newsreader@mathworks.com> wrote in message <irbt3g$8el$1@newscl01ah.mathworks.com>...
> Thanks for the tip on the vector notation for the ray and normal. I was a bit confused because I was working with points and line equations rather than vector notation. It is much easier this way, thanks.
> It was an interesting suggestion to skip the calculation of theta, however, the formulation is not straight forward for me. I do have increasing speed in my simulation (water+solid to air bubble), so my r is greater than one. If we assume that the vectors are normalized, could we define the new vector without having to separately calculate theta of refraction (although it is not a big deal)? thanks.
- - - - - - - - - -
  That code I gave you will in fact work for r greater than one. Also you do not have to normalize P and N for it to work properly. It is designed to always give the norm of Q equal to that of P.

  However to get meaningful results you do need to restrict its use to cases where the quantity within the square root,

 t1^2+t2*dot(N,N),

is non-negative. When this is zero that is the extreme case where Q would be bent over so far as to be orthogonal to the surface normal - in other words the outgoing ray would travel tangent to the surface. Beyond that, no refraction can occur - it becomes total reflection. This nicely corresponds to a complex result for the square root in the code.

  In your statement, "the formulation is not straight forward for me", are you asking how the expression was derived? It is not difficult. Let the vector Q be expressed as a linear combination of P and N:

 Q = a*P + b*N

which forces Q to lie in the same plane as P and N. Then impose two conditions on Q: 1) that the ratio of the orthogonal distance of the endpoint of Q from the line of N to the orthogonal distance of P to that line be r (Snell's Law,) and 2) that the norm of Q be equal to the norm of P. The first requirement implies that coefficient a must equal r since the ratio of the lengths of vectors cross(Q,N) and cross(P,N), being the ratios of the sines of the theta angles, must be equal to r. The second requirement gives rise to a quadratic equation in b with two solutions in which only one of them is valid for the refraction problem, and is the one used in the code I gave you. (You will note that I rationalized the numerator of this quadratic solution to provide better computation accuracy when r is close to one.)

  It is equivalent to the geometry problem of drawing a circle in the plane of N and P with radius norm(P). Vectors N and P are based at the circle's center and the endpoint of P would therefore be a point on the circle. We are to move in or out from the endpoint of vector P along the line of the vector by the ratio r, and from the resulting point draw a line parallel to the line of N until we meet the circle. This will actually occur for two points, but only the nearer point to the endpoint of N is the solution we seek. The "total reflection" case noted above corresponds to moving so far out on the line of P with an r greater than one that a line parallel to N would miss the circle altogether.

Roger Stafford

"Roger Stafford" wrote in message <irc5p0$t0t$1@newscl01ah.mathworks.com>...
> "Mahsa " <newsreader@mathworks.com> wrote in message <irbt3g$8el$1@newscl01ah.mathworks.com>...
> > Thanks for the tip on the vector notation for the ray and normal. I was a bit confused because I was working with points and line equations rather than vector notation. It is much easier this way, thanks.
> > It was an interesting suggestion to skip the calculation of theta, however, the formulation is not straight forward for me. I do have increasing speed in my simulation (water+solid to air bubble), so my r is greater than one. If we assume that the vectors are normalized, could we define the new vector without having to separately calculate theta of refraction (although it is not a big deal)? thanks.
> - - - - - - - - - -
> That code I gave you will in fact work for r greater than one. Also you do not have to normalize P and N for it to work properly. It is designed to always give the norm of Q equal to that of P.
>
> However to get meaningful results you do need to restrict its use to cases where the quantity within the square root,
>
> t1^2+t2*dot(N,N),
>
> is non-negative. When this is zero that is the extreme case where Q would be bent over so far as to be orthogonal to the surface normal - in other words the outgoing ray would travel tangent to the surface. Beyond that, no refraction can occur - it becomes total reflection. This nicely corresponds to a complex result for the square root in the code.
>
> In your statement, "the formulation is not straight forward for me", are you asking how the expression was derived? It is not difficult. Let the vector Q be expressed as a linear combination of P and N:
>
> Q = a*P + b*N
>
> which forces Q to lie in the same plane as P and N. Then impose two conditions on Q: 1) that the ratio of the orthogonal distance of the endpoint of Q from the line of N to the orthogonal distance of P to that line be r (Snell's Law,) and 2) that the norm of Q be equal to the norm of P. The first requirement implies that coefficient a must equal r since the ratio of the lengths of vectors cross(Q,N) and cross(P,N), being the ratios of the sines of the theta angles, must be equal to r. The second requirement gives rise to a quadratic equation in b with two solutions in which only one of them is valid for the refraction problem, and is the one used in the code I gave you. (You will note that I rationalized the numerator of this quadratic solution to provide better computation accuracy when r is close to one.)
>
> It is equivalent to the geometry problem of drawing a circle in the plane of N and P with radius norm(P). Vectors N and P are based at the circle's center and the endpoint of P would therefore be a point on the circle. We are to move in or out from the endpoint of vector P along the line of the vector by the ratio r, and from the resulting point draw a line parallel to the line of N until we meet the circle. This will actually occur for two points, but only the nearer point to the endpoint of N is the solution we seek. The "total reflection" case noted above corresponds to moving so far out on the line of P with an r greater than one that a line parallel to N would miss the circle altogether.
>
> Roger Stafford

Dear Roger,

Thanks for the elaboration. I can only use this formula for when the ray is traveling inside the bubble and reaches the interface of bubble and liquid to enter the liquid(r<1).
I put the two criteria, but I get a different formula, I don't know maybe am doing something different:
Q=aP+bN
a=r
t1=r*dot(P,N) & t2=(1-r^2)*dot(P,P)
norm(Q)=norm(P) ==> dot(Q,Q)=dot(P,P)
dot(Q,Q)=r^2*dot(P,P)+b^2*dot(N,N)+2*b*r*dot(P,N) ==>
(1-r^2)*dot(P,P)=b^2*dot(N,N)+2*b*r*dot(P,N)
t2=b^2*dot(N,N)+2*b*t1
t2/dot(N,N)=(b^2+2*b*t1/dot(N,N))
b^2+2*b*t1/dot(N,N)+(t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
(b+t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
(b+t1/dot(N,N))^2=t2/dot(N,N)+(t1/dot(N,N))^2=(t2*dot(N,N)+t1^2)/(dot(N,N))^2
==>b+t1/dot(N,N)=1/dot(N,N)*sqrt(t1^2+t2*dot(N,N))
b=1/dot(N,N)*(-t1+sqrt(t1^2+t2*dot(N,N))
which is a bit different than your formula of: b=t2/(sqrt(t1^2+t2*dot(N,N))+t1)
I do not need to rationalize the denominator, as have no sqrt in there, Am I making a mistake somewhere?
Thanks,
Regards,
Mahsa.

"Mahsa " <newsreader@mathworks.com> wrote in message <irdsal$jjq$1@newscl01ah.mathworks.com>...
> Dear Roger,
>
> Thanks for the elaboration. I can only use this formula for when the ray is traveling inside the bubble and reaches the interface of bubble and liquid to enter the liquid(r<1).
> I put the two criteria, but I get a different formula, I don't know maybe am doing something different:
> Q=aP+bN
> a=r
> t1=r*dot(P,N) & t2=(1-r^2)*dot(P,P)
> norm(Q)=norm(P) ==> dot(Q,Q)=dot(P,P)
> dot(Q,Q)=r^2*dot(P,P)+b^2*dot(N,N)+2*b*r*dot(P,N) ==>
> (1-r^2)*dot(P,P)=b^2*dot(N,N)+2*b*r*dot(P,N)
> t2=b^2*dot(N,N)+2*b*t1
> t2/dot(N,N)=(b^2+2*b*t1/dot(N,N))
> b^2+2*b*t1/dot(N,N)+(t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
> (b+t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
> (b+t1/dot(N,N))^2=t2/dot(N,N)+(t1/dot(N,N))^2=(t2*dot(N,N)+t1^2)/(dot(N,N))^2
> ==>b+t1/dot(N,N)=1/dot(N,N)*sqrt(t1^2+t2*dot(N,N))
> b=1/dot(N,N)*(-t1+sqrt(t1^2+t2*dot(N,N))
> which is a bit different than your formula of: b=t2/(sqrt(t1^2+t2*dot(N,N))+t1)
> I do not need to rationalize the denominator, as have no sqrt in there, Am I making a mistake somewhere?
> Thanks,
> Regards,
> Mahsa.
- - - - - - - - -
  You haven't made any mistakes, Mahsa, but yours and my expressions for b are identically equal. Multiply and divide your expression for b by

 sqrt(t1^2+t2*dot(N,N))+t1,

which wouldn't change its value, and you will get my expression for b. To see that, first do the multiplication:

 1/dot(N,N)*(-t1+sqrt(t1^2+t2*dot(N,N))*(sqrt(t1^2+t2*dot(N,N))+t1) =
 1/dot(N,N)*(sqrt(t1^2+t2*dot(N,N))-t1)*(sqrt(t1^2+t2*dot(N,N))+t1) =
 1/dot(N,N)*(t1^2+t2*dot(N,N)-t1^2) = 1/dot(N,N)*t2*dot(N,N) = t2

Now divide by the same quantity and get

 t2/(sqrt(t1^2+t2*dot(N,N))+t1)

which is my expression for b. They are exactly the same quantity.

  To settle any lingering doubts on this score try entering random vectors in the two expressions to see that the results are always equal (except for round off errors.)

  The above is a standard method of removing (rationalizing) the square root from the numerator (not denominator!) of expressions. The advantage of doing so in this case is that when r is very close to 1, your expression for b requires subtracting two sizable quantities which are nearly equal, which can created large round off errors, whereas in the other expression the round off error tends to be smaller.

  I stand by my statement that the same formula is valid for cases when r is greater than 1 as well those in which r is less than or equal to 1. Of course as I said before, when r is greater than 1 you need to cancel the refractive ray and use only reflection in situations where t1^2+t2*dot(N,N) becomes negative.

Roger Stafford

Roger,

I haven't read through all the discussion here so I might be jumping the gun a bit with this suggestion but hopefully it'll be useful...

I needed to rotate some data points for a task I was doing and checked to see if Matlab contained any functions for this.

Interestingly it has a function "rotate" which is geared towards rotating plotted data.

If you type "open rotate" at the command prompt you'll get the source code which can easily be rewritten to remove references to plot handles and add input syntax for vectors of data.

I have already done this myself and would put it on the File Exchange if it were not for the fact that the original code states "copyright Mathworks" so I don't want to go distributing it without permission.

But I will say this - my work, which I saved as a new function called "rotatedata" ended up with an input syntax that goes:

[newx,newy,newz]=rotatedata(x0,y0,z0,rotaxis,angle,origin)

where:

>> x0,y0,z0 are your original coordinates on the 3D line you want to rotate

>> rotaxis is a 3-element vector of zeros or ones that controls about which axis the points will rotate (e.g. for rotation about x-axis, rotaxis=[1,0,0] , for rotation about y and z axes put [0,1,1].

>> angle is the angle in degrees that the data will rotate by about all the axes specified by rotaxis

>> origin is a 3 element vector [x_origin, y_origin, z_origin] specifying the point about which rotation will occur

The code in "rotate" could also be easily rewritten to incorporate angle as a vector, specifying the rotation angle about the seperate axes, but I haven't done this - I just perfom the rotations separately if say I need to rotate about x by 45 degrees then about z by 20 degrees.


Hope this is useful.

Best,
Adam

"AJP" wrote in message <ire81k$sor$1@newscl01ah.mathworks.com>...
> Roger,
>
> I haven't read through all the discussion here so I might be jumping the gun a bit with this suggestion but hopefully it'll be useful...
>
> I needed to rotate some data points for a task I was doing and checked to see if Matlab contained any functions for this.
>
> Interestingly it has a function "rotate" which is geared towards rotating plotted data.
>
> If you type "open rotate" at the command prompt you'll get the source code which can easily be rewritten to remove references to plot handles and add input syntax for vectors of data.
>
> I have already done this myself and would put it on the File Exchange if it were not for the fact that the original code states "copyright Mathworks" so I don't want to go distributing it without permission.
>
> But I will say this - my work, which I saved as a new function called "rotatedata" ended up with an input syntax that goes:
>
> [newx,newy,newz]=rotatedata(x0,y0,z0,rotaxis,angle,origin)
>
> where:
>
> >> x0,y0,z0 are your original coordinates on the 3D line you want to rotate
>
> >> rotaxis is a 3-element vector of zeros or ones that controls about which axis the points will rotate (e.g. for rotation about x-axis, rotaxis=[1,0,0] , for rotation about y and z axes put [0,1,1].
>
> >> angle is the angle in degrees that the data will rotate by about all the axes specified by rotaxis
>
> >> origin is a 3 element vector [x_origin, y_origin, z_origin] specifying the point about which rotation will occur
>
> The code in "rotate" could also be easily rewritten to incorporate angle as a vector, specifying the rotation angle about the seperate axes, but I haven't done this - I just perfom the rotations separately if say I need to rotate about x by 45 degrees then about z by 20 degrees.
>
>
> Hope this is useful.
>
> Best,
> Adam

Thanks Adam,
If I define the axis of rotation as normal to the vector and define the angle, I think it would work. I will give that a try as well.
Thanks again.

"Roger Stafford" wrote in message <ire089$34u$1@newscl01ah.mathworks.com>...
> "Mahsa " <newsreader@mathworks.com> wrote in message <irdsal$jjq$1@newscl01ah.mathworks.com>...
> > Dear Roger,
> >
> > Thanks for the elaboration. I can only use this formula for when the ray is traveling inside the bubble and reaches the interface of bubble and liquid to enter the liquid(r<1).
> > I put the two criteria, but I get a different formula, I don't know maybe am doing something different:
> > Q=aP+bN
> > a=r
> > t1=r*dot(P,N) & t2=(1-r^2)*dot(P,P)
> > norm(Q)=norm(P) ==> dot(Q,Q)=dot(P,P)
> > dot(Q,Q)=r^2*dot(P,P)+b^2*dot(N,N)+2*b*r*dot(P,N) ==>
> > (1-r^2)*dot(P,P)=b^2*dot(N,N)+2*b*r*dot(P,N)
> > t2=b^2*dot(N,N)+2*b*t1
> > t2/dot(N,N)=(b^2+2*b*t1/dot(N,N))
> > b^2+2*b*t1/dot(N,N)+(t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
> > (b+t1/dot(N,N))^2-(t1/dot(N,N))^2=t2/dot(N,N)
> > (b+t1/dot(N,N))^2=t2/dot(N,N)+(t1/dot(N,N))^2=(t2*dot(N,N)+t1^2)/(dot(N,N))^2
> > ==>b+t1/dot(N,N)=1/dot(N,N)*sqrt(t1^2+t2*dot(N,N))
> > b=1/dot(N,N)*(-t1+sqrt(t1^2+t2*dot(N,N))
> > which is a bit different than your formula of: b=t2/(sqrt(t1^2+t2*dot(N,N))+t1)
> > I do not need to rationalize the denominator, as have no sqrt in there, Am I making a mistake somewhere?
> > Thanks,
> > Regards,
> > Mahsa.
> - - - - - - - - -
> You haven't made any mistakes, Mahsa, but yours and my expressions for b are identically equal. Multiply and divide your expression for b by
>
> sqrt(t1^2+t2*dot(N,N))+t1,
>
> which wouldn't change its value, and you will get my expression for b. To see that, first do the multiplication:
>
> 1/dot(N,N)*(-t1+sqrt(t1^2+t2*dot(N,N))*(sqrt(t1^2+t2*dot(N,N))+t1) =
> 1/dot(N,N)*(sqrt(t1^2+t2*dot(N,N))-t1)*(sqrt(t1^2+t2*dot(N,N))+t1) =
> 1/dot(N,N)*(t1^2+t2*dot(N,N)-t1^2) = 1/dot(N,N)*t2*dot(N,N) = t2
>
> Now divide by the same quantity and get
>
> t2/(sqrt(t1^2+t2*dot(N,N))+t1)
>
> which is my expression for b. They are exactly the same quantity.
>
> To settle any lingering doubts on this score try entering random vectors in the two expressions to see that the results are always equal (except for round off errors.)
>
> The above is a standard method of removing (rationalizing) the square root from the numerator (not denominator!) of expressions. The advantage of doing so in this case is that when r is very close to 1, your expression for b requires subtracting two sizable quantities which are nearly equal, which can created large round off errors, whereas in the other expression the round off error tends to be smaller.
>
> I stand by my statement that the same formula is valid for cases when r is greater than 1 as well those in which r is less than or equal to 1. Of course as I said before, when r is greater than 1 you need to cancel the refractive ray and use only reflection in situations where t1^2+t2*dot(N,N) becomes negative.
>
> Roger Stafford
Hi,
Thanks. I guess I misread the previous email (numerator vs. denominator). It is indeed correct. The rationale behind derivation is valid for refraction with any r, I agree with you. I put a condition of if dot(P,N)/[norm(P)*norm(N)]<theta_critical to avoid getting negative value under sqrt and avoid calculation for when total reflection occurs.

You have been very helpful. Thanks again.