clc;
clear all
title('胞映射')
xl=-1.025;xu=1.525;yl=-2.535;yu=1.035;
nxc=101;nyc=101;
h1=(xu-xl)./nxc;h2=(yu-yl)./nyc;
nc=nxc.*nyc+1;n=1;
xl1=xl-h1./2.0;yl1=yl-h2./2.0;
for i=1:nyc
    for j=1:nxc
        ij=j+(i-1).*nxc;z1(ij)=j;z2(ij)=i;
    end
end
 z1(nc)=0;z2(nc)=0;
for i=1:nc
    gr(i)=0;
end
for ii=1:nc
    m=0;g=gr(ii);
    if g==0
        gr(ii)=-1;b1=z1(ii);b2=z2(ii);
        c1=b1;c2=b2;v=1;
        while v==1
            m=m+1;
            xc=c1.*h1-h1./2.0;yc=c2.*h2-h2./2.0;
            c1=fix((0.9*(yc+yl)+1.81*(xc+xl).^2-xl)./h1+0.5);
            c2=fix((-0.9*xc-0.9*xl-yl)./h2+0.5);
            if c1>nxc|c1<1|c2>nyc|c2<1;
                c1=0;c2=0;
          else
            end
          e1=c1;e2=c2;ij=c1+(c2-1).*nxc;
          if c1==0
              ij=nc;
          else
          end
          g=gr(ij);
          if g~=0;
              if g==-1;
                  i=0;
                  if b1~=e1|b2~=c2;
                      c1=b1;c2=b2;w=1;
                      while  w==1;
                          i=i+1;xc=c1.*h1-h1./2.0;yc=c2.*h2-h2./2.0;
      c1=fix((0.9*(yc+yl)+1.81*(xc+xl).^2-xl)./h1+0.5);
       c2=fix((-0.9*xc-0.9*xl-yl)./h2+0.5);
       if c1>nxc|c1<1|c2>nyc|c2<1;
           c1=0;c2=0;
       else
       end
       if c1==e1&c2==e2;
           break
       else
       end
             end
                  else
                  end
             if c1==0
                 gg=1;pe=1;
             else
                 n=n+1;pe=m-1;gg=n;
             end
             c1=b1;c2=b2;
             ij=c1+(c2-1).*nxc;
             if c1==0;
                 ij=nc;
             else
             end
             gr(ij)=gg;p(ij)=pe;s(ij)=i;
             if i~=0
                 for j=1:i-1
                     xc=c1.*h1-h1./2.0;    yc=c2.*h2-h2./2.0;
                     c1=fix((0.9*(yc+yl)+1.81*(xc+xl).^2-xl)./h1+0.5);
                    c2=fix((-0.9*xc-0.9*xl-yl)./h2+0.5);
                       if c1>nxc|c1<1|c2>nyc|c2<1;
                    c1=0;c2=0;
                     else
                       end
         ij=c1+(c2-1).*nxc;
            if c1==0|c2==0;
                ij=nc;
            else
            end
            gr(ij)=gg;p(ij)=pe;s(ij)=i-j;
                 end
                 i=1;
             end
             for j=1:m-1
                 xc=c1.*h1-h1./2.0;yc=c2.*h2-h2./2.0;
                 c1=fix((0.9*(yc+yl)+1.81*(xc+xl).^2-xl)./h1+0.5);
       c2=fix((-0.9*xc-0.9*xl-yl)./h2+0.5);
                    if c1>nxc|c1<1|c2>nyc|c2<1;
                    c1=0;c2=0;
                     else
                       end
          ij=c1+(c2-1).*nxc;
             if c1==0;
                 ij=nc;
             else
             end
             gr(ij)=gg;p(ij)=pe;s(ij)=0;   
             end
              else
                  g=gr(ij);pe=p(ij);st=s(ij);
                  c1=b1;c2=b2;ij=c1+(c2-1).*nxc;
                  if c1==0;
                      ij=nc;
                  else
                  end
                      gr(ij)=g;p(ij)=pe;s(ij)=st+m;
                      for j=1:m-1
                          xc=c1.*h1-h1./2.0;yc=c2.*h2-h2./2.0;
                   c1=fix((0.9*(yc+yl)+1.81*(xc+xl).^2-xl)./h1+0.5);
             c2=fix((-0.9*xc-0.9*xl-yl)./h2+0.5);
                    if c1>nxc|c1<1|c2>nyc|c2<1;
                    c1=0;c2=0;
                     else
                    end
                   ij=c1+(c2-1).*nxc;
             if c1==0;
                 ij=nc;
             else
             end  
               gr(ij)=g;p(ij)=pe;s(ij)=st+m-j;
                      end
              end
              break
          else
              gr(ij)=-1;
          end
        end
    else
    end
end
ax=112;ay=46;
for ii=2:n
 for i=1:nc-1
    if gr(i)==ii
        x=z1(i).*h1+xl1;y=z2(i).*h2+yl1;
        xi=x.*ax; yi=y.*ay;
        hold on
        plot(xi,yi,'.k');
        xc=xi./8;yc=yi./8;
    else
    end
end
end