- 积分
- 8
- 注册时间
- 2006-9-14
- 仿真币
-
- 最后登录
- 1970-1-1
|
此函数产生的椭圆通过调整长径比可以用来模拟纤维混凝土中的纤维,当然也可以模拟类似的其它材料。
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;
备注:本程序里虽然采用椭圆个数作为判据,但是完全可以通过面积率来作为判据,稍加转化即可。当然三维亦可采用体积率作为判据。
参考文献:http://okok.org/cgi-bin/ut/topic_show.cgi?id=127992&h=1&bpg=1&age=30
[ 本帖最后由 chengweifeng 于 2007-7-10 16:49 编辑 ] |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
评分
-
1
查看全部评分
-
|