% 本程序為用胞化積分軌跡法進行全局性態分析

% 初始化，劃分胞化空間

clear
xcellnumber=10;
ycellnumber=10;
xleft=-5;
xright=5;
ytop=10;
ybottom=-10;
di=1;
width=(xright-xleft)/xcellnumber;
height=(ytop-ybottom)/ycellnumber;
center(1:xcellnumber/di,1)=[(xleft+width/2):di*width:(xright-width/2)]';
center(1:ycellnumber/di,2)=[(ybottom+height/2):di*height:(ytop-height/2)]';
Id=zeros(xcellnumber*ycellnumber,1);
Ib=zeros(xcellnumber*ycellnumber,1);
Y=zeros(xcellnumber*ycellnumber,2);
dpx=0.001;
dpy=0.001;
nm=10001;
ipsilong=0.0001;
row=0;
periodcells=0;
for i=1:10
    for k=1:10
        initvalue((i-1)*10+k,1)=center(i,1);
        initvalue((i-1)*10+k,2)=center(k,2);
    end
end
plot(initvalue(:,1),initvalue(:,2),'.');
% omiga=4.7547;
% B=58.4;
% G=311.5;
% eita=0.1;
% K=1;
% U=1e-6;
G=2;
omiga=3.9311;
B=2;
K=1;
U=1;
eita=0.1;

% 從每個胞的中心點開始進行映射。

for i=1:ycellnumber/di
   i;
   cput=cputime;
    for j=1:xcellnumber/di
         j;
        trajectory=zeros(1,2);             
        sequence=zeros;
        initial(1,1)=center(j,1);
        initial(1,2)=center(i,2);
        ci=(i-1)*di*xcellnumber+(j-1)*di+1;
        Ib(ci)=Ib(ci)+1;

       
        if Id(ci)>0
            if abs(Y(ci,1)-initial(1,1))<=dpx&abs(Y(ci,2)-initial(1,2))<=dpy
                continue;
            else
                Id(ci)=-Id(ci);
            end
        else               
            Id(ci)=-nm;
            Y(ci,:)=[initial(1,1),initial(1,2)];
        end
        trajectory(1,:)=initial(1,:);
        sequence(1)=ci;
        count1=0;
        count2=1;
        for count=2:2000
            if count==2000
                Id(sequence)=-Id(sequence);
                break
            end
                
%             [t,x]=ode45('duffing',[0 2*pi/omiga],[initial(1,1),initial(1,2)]);
%             
            opts=simset('Initialstate',[initial(1,1),initial(1,2)]);
%             [t,x,y]=sim('vanderpol',[2],opts);
            [t,x,y]=sim('duffingscm',[2*pi/omiga],opts);
            [m,e]=size(t);
            
               % 當映射的像點不在選定區域
 
            if x(m,1)>=xright|x(m,1)<xleft|x(m,2)<ybottom|x(m,2)>=ytop      
                if count1~=count-1
                    k=1;
                    count1=0;
                else
                    k=k+1;
                end                              
                if k>12
                    for index=1:count2
                        cj=sequence(index);
                        if Id(cj)==-nm
                            Id(cj)=nm+1;
                        else if Id(cj)<0
                                Id(cj)=-Id(cj);
                            end
                        end
                    end
                   break
                end
                initial=x(m,:);
                count1=count;
                continue;
            end
            
            count2=count2+1;
            indexx=fix(abs(x(m,1)-xleft)/width)+1;
            indexy=fix(abs(x(m,2)-ybottom)/height)+1;
            
            ci=(indexy-1)*ycellnumber+indexx;
            sequence(count2)=ci;
            trajectory(count2,:)=[x(m,1),x(m,2)];
            Ib(ci)=Ib(ci)+1;
            
            
         % 定義四種狀態
         
            if Id(ci)==0
                state=1;
            else
                if Id(ci)==-nm
                    state=2;
                else
                    if Id(ci)>0
                        state=3;
                    else
                        state=4;
                    end
                end
            end
                    
         % 對四個狀態分別進行處理
            
            switch state
            case 1
                Id(ci)=-nm;
                Y(ci,:)=[x(m,1),x(m,2)];
                initial=x(m,:);
            case 2
                if sqrt(((Y(ci,1)-x(m,1))^2+(Y(ci,2)-x(m,2))^2))<=ipsilong
                    period=0;                 
                   for index=count2-1:-1:1
                        if sequence(index)==ci
                            period(1,1:(count2-index))=sequence(index+1:count2);
                            break
                        end
                    end
                    [r2,c2]=size(period);
                    isold=0;
                    C=find(Id<xcellnumber*ycellnumber);
                    Idd=Id(C);
                    for index=max(Idd(1,:)):-1:1
                        oldperiod=isoldperiod(period,periodcells);
                        if oldperiod
                            Id(sequence)=index;  
                            isold=1;
                            break
                        end
                    end
                    if isold~=1
                        periodcells(max(Idd)+1,1:c2)=period;
                        Id(sequence)=max(Idd)+1;
                        
                        break
                    end
                else 
                    Y(ci,:)=[x(m,1),x(m,2)];
                    initial=x(m,:);
                end
            case 3
                if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
                    for index=1:count2
                        cj=sequence(index);
                        if Id(cj)==-nm|Id(cj)==-Id(ci)
                            Id(cj)=Id(ci);
                        else 
                            if Id(cj)<0&Id(cj)~=-nm&Id(cj)~=-Id(ci)
                                Id(cj)=-Id(cj);
                            end
                        end
                    end
                    break
                else
                    Id(ci)=-Id(ci);
                    initial=x(m,:);
                end
            otherwise
                if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
                    Id(ci)=-Id(ci);
                    for index=1:count2
                        cj=sequence(index);
                        if Id(cj)==-nm|Id(cj)==-Id(ci)
                            Id(cj)=Id(ci);
                        else
                            if Id(cj)<0&Id(cj)>-nm&Id(cj)~=-Id(ci)
                                Id(cj)=-Id(cj);
                            end
                        end
                    end
                    break
                else
                    period=0;
                    for index=count2-1:-1:1
                        if sequence(index)==ci
                            period(1,1:(count2-index))=sequence(index+1:count2);
                            Y1=trajectory(index,:);
                            break
                        end
                    end
                    if sqrt(((Y1(1,1)-x(m,1))^2+(Y1(1,2)-x(m,2))^2))<=ipsilong
                        [r2,c2]=size(period);                  
                        isold=0;
                        C=find(Id<xcellnumber*ycellnumber);
                        Idd=Id(C);
                        for index=max(Idd):-1:1
                            oldperiod=isoldperiod(period,periodcells(index,:));
                            if oldperiod
                                Id(sequence)=index;  
                                isold=1;
                                break
                            end
                        end
                        if isold~=1
                            periodcells(max(Idd)+1,1:c2)=period;
                            Id(sequence)=max(Idd)+1;                                                   
                        end
                        break
                    end
                end
                initial=x(m,:);
            end
        end
    end
    et=(cputime-cput)/60;
end