Wednesday, October 13, 2010

Plot a heart curve in matlab

 
function ShowHeart(Ntype)
% a funny script to show different kind of heart curves
% Usage:
%        ShowHeart()
% Copyright: http://scriptdemo.blogspot.com
% Ref: http://www.mathematische-basteleien.de/heart.htm
%        http://mathworld.wolfram.com/HeartCurve.html

if nargin==0 Ntype=1; end

np=999;Is3D=0;IsFilled=1;
switch Ntype
    case {1}
       t=0:2*pi/np:2*pi;
       xx=(1-sin(t)).*cos(t);
       yy=(1-sin(t)).*sin(t);
    case {2}
       t=-1:2/(np*100):1;
       xx=sin(t).*cos(t).*log(abs(t));
       yy=sqrt(cos(t)).*abs(t).^0.3;
       yy(isnan(xx))=[];
       xx(isnan(xx))=[];
    case {3}
       t=-1:2/(np*100):1;
       xx=sin(t).*cos(t).*log(abs(t));
       yy=sqrt(abs(t)).*cos(t);
       yy(isnan(xx))=[];
       xx(isnan(xx))=[];
    case {4}
       t=0:2*pi/np:2*pi;
       r=sin(t).*sqrt(abs(cos(t)))./(sin(t)+7/5)-2*sin(t)+2;
       xx=r.*cos(t); yy=r.*sin(t); clear r
    case {5}
       t=0:60/np:60;
       r=0.01*(-t.*t+40*t+1200);
       xx=r.*sin(pi*t/180);
       yy=r.*cos(pi*t/180); clear r t
       xx=[xx fliplr(-xx)]; yy=[yy fliplr(yy)];
    case {6}
       t=-1:2/(np*100):1;
       r=(1-abs(t)).*(1+3*abs(t));
       xx=r.*sin(t); yy=r.*cos(t); clear r t
    otherwise
       disp(['Sorry, you just broke my heart!']);
       Is3D=1;
       t=0:2*pi/np:2*pi;
       r=sin(t).*sqrt(abs(cos(t)))./(sin(t)+7/5)-2*sin(t)+2;
       xx=r.*cos(t); yy=r.*sin(t); clear r
end
if ~Is3D
   if ~IsFilled
      plot(xx,yy,'color','r');
   else
      hf=fill(xx,yy,'r');set(hf,'linestyle','none');
   end
else
  hp=plot(xx,yy,'g');set(hp,'linestyle','none');
  figurenicer; hold on;
  myxlim=xlim;myylim=ylim;
  nrp=50;
  nxgrid=20;
  nygrid=20;
  xxn=linspace(myxlim(1),myxlim(2),nxgrid);
  yyn=linspace(myylim(1),myylim(2),nygrid);
  [xxn,yyn]=meshgrid(xxn,yyn);
  tmprand=randn(nrp,2);
  tmprand=tmprand/max(abs(tmprand(:)))/2;
  xxrand=zeros(nrp,nxgrid*nygrid);
  yyrand=xxrand;
  np=1;
  for Nx=1:nxgrid
      for Ny=1:nygrid
         xxrand(:,np)=xxn(Ny,Nx)+tmprand(:,1);
         yyrand(:,np)=yyn(Ny,Nx)+tmprand(:,2);
         np=np+1;
      end
  end
  Ind=inpolygon(xxn(:),yyn(:),xx(:),yy(:));
  tmpind=1:length(xxn(:));
  [jj,ii]=ind2sub(size(xxn),tmpind(Ind==1));
  Ind2=sub2ind(size(xxn),jj(:)+1,ii(:));
  Ind=sub2ind(size(xxn),jj(:),ii(:));
  xxrand(:,Ind2)=xxrand(:,Ind);
  yyrand(:,Ind2)=yyrand(:,Ind);
  
  [xxnew,yynew]=GetRandom(xx(:),yy(:),5);
  hp0=plot(xxnew,yynew,'g.');set(hp0,'linestyle','none','markersize',1);
  %plot(xxn(:),yyn(:),'.b');
  hp=plot(xxrand,yyrand,'.g');set(hp,'linestyle','none','markersize',1);
end
figurenicer;
axis equal;
axis tight;
function [xxnew,yynew]=GetRandom(xx,yy,nrp)
 np=length(xx);
 if (np>360)
    iii=linspace(1,np,360);
    iii=fix(iii);
    xx=xx(iii); yy=yy(iii);
    np=length(xx);
 end
 tmprand=randn(nrp,2);
 tmprand=tmprand/max(abs(tmprand(:)))*8;
 xxnew=zeros(nrp,length(xx));
 xx=reshape(xx,1,[]);
 yy=reshape(yy,1,[]);
 xxdif=diff(xx);xxdif=[xxdif(1) xxdif];
 yydif=diff(yy);yydif=[yydif(1) yydif];
 for N=2:length(xx)
     xxnew(:,N)=0.5*(xx(N-1)+xx(N))+tmprand(:,1)*xxdif(N);
     yynew(:,N)=0.5*(yy(N-1)+yy(N))+tmprand(:,2)*yydif(N);
 end
xxnew(:,1)=xxnew(:,2);
yynew(:,1)=yynew(:,2);

No comments:

ShowCalendar