aUCBLogo Demos and Tests / solve_cosines

to solve_cosines
action stop
   WindowMode
   setUpdateGraph "false
   a=0
   forever
    [   clearScreen
      ;pr repcount
      catch "error [ignore runResult [action]]
      err=Error
      if :err != []
      [   clearText
         print :err.2
         for [i 1 100]
         [   dispatchMessages
            waitMS 10
         ]
      ]
      updategraph
      updateVars
      for [i 1 100]
      [   dispatchMessages
         waitMS 10
         if Key? [stop]
      ]
   ]
end

to action
;   drawit
;   drawf
   solveit
;   solveitsymbolically
end

to drawit
   U=300
   delta=2
   for [phi 0 360 10]
   [   U1=U*(cos phi)
      U2=U*(cos phi+delta)*cos delta
      U3=U*(cos phi+2*delta)*sqr cos delta

      x=phi/360*800-400
      plot x U1 phi
      plot x U2 phi
      plot x U3 phi
   ]
end

to drawf
   clearScreen
   U=1
   phi=20
   delta=30
   U1=U*(cos phi)
   U2=U*(cos phi+delta)*cos delta
   U3=U*(cos phi+2*delta)*sqr cos delta
   for [phi 0.1 359.9 1]
   [   x=phi/360*800-400
      catch "error
      [   y=f x
         plot x (abs y)*1e16 phi
      ]
   ]
   updateGraph
end

to plot x y c
   setFloodColor HSB c 1 1
   pu   setXY x y pd fillCircle 2
end

to solveitsymbolically
   clearText
   WindowMode
   clearScreen
   U=1
   for [phi 0.1 45 20]
   [   for [delta 0.1 45 20]
      [   U1=U*(cos phi)
         U2=U*(cos phi+delta)*(cos delta)
         U3=U*(cos phi+2*delta)*(sqr cos delta)
         catch "error
         [   a=(sqr U2)*(2*U2-U1-U3)/(U3*(sqr U1-U2)+(sqr U2)*(2*U2-U1-U3))
            p = float real arccos (sqrt a)
            delta2 = float real arccos sqrt U3/(2*U2-U1)
            phi2=p-delta2
            V1=U*(cos phi2)
            V2=U*(cos phi2+delta2)*(cos delta2)
            V3=U*(cos phi2+2*delta2)*(sqr cos delta2)
            ifelse (abs U1-V1) > 1e-4
            or2 (abs U2-V2) > 1e-4
            or2 (abs U3-V3) > 1e-4
            [   pr "\ bug
            ][   (pr)
            ]
            (type "\    U1 "\    V1 "\    U2 "\    V2 "\    U3 "\    V3)
         ]
      ]
   ]
end

to solveit
   clearText
   WindowMode
   clearScreen
   U=1
   for [phi 0.1 45 20]
   [   for [delta 0.1 45 20]
      [   U1=U*(cos phi)
         U2=U*(cos phi+delta)*(cos delta)
         U3=U*(cos phi+2*delta)*(sqr cos delta)
         findZeros
         V1=U*(cos phi2)
         V2=U*(cos phi2+delta2)*(cos delta2)
         V3=U*(cos phi2+2*delta2)*(sqr cos delta2)
         ifelse (abs (U1-V1)/U1) > 1e-3
         or2 (abs (U2-V2)/U2) > 1e-3
         or2 (abs (U3-V3)/U3) > 1e-3
         [   pr "\ bug
         ][   (pr)
         ]
         (type "\    U1 "\    V1 "\    U2 "\    V2 "\    U3 "\    V3)
      ]
   ]
end

to findZeros
   epsilon=1e-4
   za=rnd*360
   zb=rnd*360
   y=f za zb
   norefresh
   i=0
   k=0
   fc=1/10
   while [(abs y) > epsilon]
   [   da=df_da za zb
      db=df_db za zb
;(pr da db)
      ifelse k>5000
      [   za=rnd*360
         zb=rnd*360
         y=f za zb
         k=0
pr "oh!
      ][   za1=za-min 10 max -10 y/da*fc
         zb1=zb-min 10 max -10 y/db*fc
         za1=(mod za1 360)
         zb1=(mod zb1 360)
         y1=f za1 zb1
         ifelse y1>y
         [   za1=za+min 10 max -10 y/da*fc
            zb1=zb-min 10 max -10 y/db*fc
            za1=(mod za1 360)
            zb1=(mod zb1 360)
            y1=f za1 zb1
            ifelse y1>y
            [   za1=za-min 10 max -10 y/da*fc
               zb1=zb+min 10 max -10 y/db*fc
               za1=(mod za1 360)
               zb1=(mod zb1 360)
               y1=f za1 zb1
               ifelse y1>y
               [   za1=za+min 10 max -10 y/da*fc
                  zb1=zb+min 10 max -10 y/db*fc
                  za=(mod za1 360)
                  zb=(mod zb1 360)
                  y=f za zb
               ][   za=za1
                  zb=zb1
                  y=y1
               ]
            ][   za=za1
               zb=zb1
               y=y1
            ]
         ][   za=za1
            zb=zb1
            y=y1
         ]
      ]
;pr y
      setPixelXY
         (float za)-180
         (float zb)-180
         HSB 100*abs y 1 1
      if i==1000
      [   i=0
         updateGraph
      ]
      i=i+1
      k=k+1
   ]
   phi2=za
   delta2=zb
end

to f a b
   z=sqrt (sqr (cos a+b)*(cos b)/(cos a)-U2/U1)
      +sqr ((cos a+2*b)*(sqr cos b)/(cos a)-U3/U1)
   output z
end

to df_da a b
   local [y c]
   y=f a b
   c=(y-(f a+1e-6 b))/1e-6
   output c
end

to df_db a b
   local [y d]
   y=f a b
   d=(y-(f a b+1e-6))/1e-6
   output d
end