此函数产生的椭圆通过调整长径比可以用来模拟纤维混凝土中的纤维，当然也可以模拟类似的其它材料。- function myellipse
- % cheng wf 2007 6 11 19:51
- close all
- clear all
- clc
- %format compact
- depth=60;
- a=[3 2 1 0.5];b=[2 1 0.5 0.25];
- number=[10 10 100 20];
- x0=[];y0=[];xlat=[];ylon=[];
- [x0new1,y0new1,xlat1,ylon1]=myellipse_one(depth,a(1),b(1),number(1),x0,y0,xlat,ylon)
- x0=x0new1;y0=y0new1;xlat=xlat1;ylon=ylon1;
- for i=2:length(a)
- [x0new2,y0new2,xlat2,ylon2]=myellipse_rest(depth,a(i),b(i),number(i),x0,y0,xlat,ylon)
- x0=x0new2;y0=y0new2;xlat=xlat2;ylon=ylon2;
- end
- % 绘制椭圆图形
- [m,n]=size(xlat)
- plot(xlat,ylon)
- axis([0 60 0 60])
- function [x0new,y0new,xlat,ylon]=myellipse_one(depth,a,b,number,x0,y0,xlat,ylon)
- xxlen=length(x0);yylen=length(y0);
- x0(xxlen+1)=a+(depth-2*a)*rand;
- y0(yylen+1)=a+(depth-2*a)*rand;
- ecc=axes2ecc(a,b);
- [lat1,lon1]=ellipse1(x0(xxlen+1),y0(yylen+1),[a,ecc],rand*180);
- xlat=[xlat,lat1];ylon=[ylon,lon1];
- for i=2:number
- while 1
- x0(xxlen+i)=a+(depth-2*a)*rand;
- y0(yylen+i)=a+(depth-2*a)*rand;
- [lati,loni]=ellipse1(x0(xxlen+i),y0(yylen+i),[a,ecc],rand*180);
- xlat=[xlat,lati];ylon=[ylon,loni];
- for j=1:xxlen+i-1
- in1=[];in2=[];
- in1=inpolygon(xlat(:,xxlen+i),ylon(:,yylen+i),xlat(:,j),ylon(:,j))
- in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i))
- if all(in1==0)&&all(in2==0)
- n_j=j;
- if n_j==xxlen+i-1
- N=1
- break
- end
- continue
- else
- N=2;
- xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
- break
- end
- end
- if N==1
- break
- end
- end
- end
- x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;
- function [x0new,y0new,xlat,ylon]=myellipse_rest(depth,a,b,number,x0,y0,xlat,ylon)
- xxlen=length(x0);yylen=length(y0);
- ecc=axes2ecc(a,b);N=0;
- for i=1:number
- while 1
- x0(xxlen+i)=a+(depth-2*a)*rand;
- y0(xxlen+i)=a+(depth-2*a)*rand;
- [lat,lon]=ellipse1(x0(xxlen+i),y0(xxlen+i),[a,ecc],rand*180);
- xlat=[xlat,lat];ylon=[ylon,lon];
- for j=1:xxlen+i-1
- in1=[];in2=[];
- in1=inpolygon(xlat(:,xxlen+i),ylon(:,xxlen+i),xlat(:,j),ylon(:,j));
- in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i));
- if all(in1==0)&&all(in2==0)
- n_j=j;
- if n_j==xxlen+i-1
- N=1;
- break
- end
- continue
- else
- N=2;
- xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
- break
- end
- end
- if N==1
- break
- end
- end
- end
- x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;
|
