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 [1 100]
         
[   dispatchMessages
            
waitMS 10
         
]
      
]
      
updategraph
      
updateVars
      
for [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 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))
            
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 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-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