Find the angle between two tangents by iterating en2

Find the angle between two tangents by iterating the distant between the
focal point and the vertex of the parabola, d=en2. I have the following
parabolic equation x^2+d*y*4+d^2*(-4) and want to find d=en2 when the angle
between two tangents is 0.1674317785. I know manually that.

 

(%i1) x^2+d*y*4+d^2*(-4)$

(%i2) e1:  rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));

(%i3) en1: 10$

(%i4) en2: 5$

(%i5) %phi[1] : first(e2)$

(%i6) %phi[1] : first(e2)$

(%i7) %phi[2] : second(e2)$

(%i8) [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))]$

          [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))]$

 
((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*
(-1)));

(%i9) %,numer;

 

This gives me the angle (%i8) between the two tangents as 1.5707963268. Now
if I keep incrementing d=en2 by 10 it will get me very close, i.e.,

en2: 15 gives me ang: 0.6435011088

en2: 25 gives me ang: 0.3947911197

etc. en2.

 en2: 55 gives me ang: 0.1813197744

 

Now if I switch from a for loop method to say the Newton method then this
should get closer to a value with say 10 decimal places. Problem is I do not
know how to set this up. I think it is something like the following, but I
am not sure. 

 

x^2+d*y*4+d^2*(-4)$

e1:  rhs(first(solve(x^2+d*y*4+d^2*(-4),y)));

en1: 10$                        /* set x = 10   */

f1(ang) := block([n, en2:5],

              e2: subst([x=en1, d=en2],[first(diff(e1,x)), diff(e1,x)]),

              %phi[1] : first(e2),

              %phi[2] : second(e2),

              ang : [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))],

                       [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))],

 
((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*
(-1))),

              for i:1 thru n while ang # 0.178092938 do en2 : en2+i,

              return (en2));         /*Should return en2:55. I do not
understand what's going on here. */

newton(exp,var,x0,eps):=

   block([en2:5,e2,xn,s,numer:true],

     s:diff(exp,var),

     xn:x0,

   loop,

     if abs(subst(xn,var,exp))<eps then return(xn),

     xn:xn-subst(xn,var,exp)/subst(xn,var,s),

     e2: subst([x=en1, d=en2],[first(diff(e1,x)),diff(e1,x)]),

     %phi[1] : first(e2),

     %phi[2] : second(e2),

     ang : [(%phi[1]*d*(-2)),(d+((%phi[1])^(2)*d*(-1)))],

           [(%phi[2]*d*(-2)),(d+((%phi[2])^(2)*d*(-1)))],

 
((%pi*1/2)+(atan((((%phi[1]+(%phi[2]*(-1))))^((-1))*(1+(%phi[2]*%phi[1]))))*
(-1))),

     go(loop))$

h:expand(e1);   /*I am not sure what goes in here? In either case it does
not work.  */

newton(h,x,3.5,1.0e-10);         /*Last value should be the angle
0.1674317785 */

 

Can someone help me with the algorithms?

 

 

William Porter

(425) 343-7711

wporter at omegapar.com