%% nsection_Newton_spmd.m
% 
Nlabs=2;
%Nlabs=16;
% Example functions
%% 1)
 f=@sin;df=@cos;
 interval=[-5*pi-.1 5*pi];
 n=10;
%% 2) 
f=@(x) besselj(0,x); df=@(x) -besselj(1,x);
interval=[0 100];
n=20;
%% 3) 
 f=@(x) exp(x)-1.5-atan(x)
 interval=[30 3];
 syms x
 diff(f(x),x)
%{ 
ans = 
exp(x) - 1/(x^2 + 1)
%}
 df=@(x) exp(x) - 1./(x.^2 + 1);
%{ 
 close all
 subplot(2,1,1)
 ezplot(f,[ -6 2]),grid on    % x0=0.7
 subplot(2,1,2)
 ezplot(f,[ -16 -6]), grid on  % x0=-14
 %}
% ezplot(f,[-30 -16])% ezplot(f,[-100 -30 ]) % ezplot(f,[-1000 -100 ])
%   No more zeros, it seems.
%% 4)
 f=@(x) log(x+1)+tan(2*x)
%% 5)
interval=[90 200]
 f=@(x) 50*pi+sin(x)-100*atan(x)
%
%{ 
 subplot(2,2,1)
 ezplot(f,[50 100])
 grid on
 subplot(2,2,2)
 ezplot(f,[100 200])
 grid on
 subplot(2,2,3)
 ezplot(f,[200 300])
 grid on
 subplot(2,2,4)
 ezplot(f,[300 500])
 grid on
%}
%
% >> syms x
% >> diff(f(x),x) 
% ans = cos(x) - 100/(x^2 + 1)
df = @(x) cos(x) - 100 ./(x.^2 + 1)
Nlabs=2;
%%
% parpool('local',Nlabs)
% Devide "interval" into Nlabs subintervals.
% Find brackets for zeros in each.
%
%%
tic
spmd
    %n=10;
    n=20;   % Only 15 roots
    n=30;   %  25 roots, 6 roots missing
    n=40;   %  31 roots, all captured 
    n=50;   %  29, first 2 missing, why????
    n=51;   %  all 31
    n=41    % all 31
   
    n=100  % all
    n=1000 % all
    
   % n=49    % 29, misses first 2
   % n=48    % 29, misses 3rd and 4th
    
             % 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
         myroots=nan;  % NOTE: Don't use "roots" as a variable name !
         Niter=nan;
     else   
         [myroots, Niter]=mynewt(f,c,df);
     end       
end
%
T=toc
x=cat(1,myroots{:})   % Elegant way to avoid the following for-loop.
allroots=x(~isnan(x))
%% 
%{ 
clf
ezplot(f,[90 200])
hold on
plot(allroots,f(allroots),'*')
grid on
%}
% For Nlabs=2,n=50 
%%
%{
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)]
%%
%{
close all
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) 


