%% nsection_Newton_spmd.m
% 
Nlabs=2;
%Nlabs=16;
interval=[-5*pi-.1 5*pi];
% Example functions
f=@sin;df=@cos;
%f=@(x) besselj(0,x); df=@(x) -besselj(1,x);

parpool('local',Nlabs)
% Devide "interval" into Nlabs subintervals.
% Find brackets for zeros in each.
%%
spmd
    n=10;   
             % Divide interval into n*Nlabs parts. 
             % Choose n 
             % 1) so as not to waste comp. power for extra accuracy as 
             %    starting points to Newton
             % 2) Big enough so as not to miss some of the brackets
             %    for zeros.
             
    Len=(interval(2)-interval(1))/Nlabs;
    a=interval(1)+(labindex-1)*Len;
    Labinterval=[a,a+Len];
    
    %diff(sign(f(Labinterval)))~=0
    [br,c]=nsection(f,Labinterval,n);   
     if prod(size(c))==0
         roots=nan;
         Niter=nan;
     else   
         [roots, Niter]=mynewt(f,c,df);
     end       
end
%%
x=cat(1,roots{:})   % Elegant way to avoid the following for-loop.
allroots=x(~isnan(x))
%% 
%{
x=zeros(1,Nlabs);
for k=1:Nlabs
         x(k)=roots{k};
end
%}

sortroots=sort(allroots);
format long
%%
uniqueroots=sortroots(logical([1; abs(diff(sortroots))>1e-10]))
%%
[uniqueroots f(uniqueroots)]
%%
ezplot(f, interval)
hold on
plot(uniqueroots,0,'xk')
shg

%%
%{
format long
uniqueroots{:}
R1=uniqueroots{1}'
R2=uniqueroots{2}'
%}
%{
T1=[R1 (-5*pi:pi:0)']
T2=[R2 (0:pi:5*pi)']


diff(T1,[2])
diff(T2,[2])

%}
%%
%(-5*pi:pi:0) - (roots{1})
%%
%[-5*pi:pi:0] - (roots{1})
%{
a=-5*pi:pi:0
b=[-5*pi:pi:0]
a-b

-5*pi:pi:0 - [-5*pi:pi:0]
%}
%%
delete(gcp)


