声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1326|回复: 0

[编程技巧] 【分享】不相交的随机大小的椭圆

[复制链接]
发表于 2012-9-11 09:23 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
此函数产生的椭圆通过调整长径比可以用来模拟纤维混凝土中的纤维,当然也可以模拟类似的其它材料。
  1. function myellipse
  2. % cheng wf 2007 6 11 19:51
  3. close all
  4. clear all
  5. clc
  6. %format compact
  7. depth=60;
  8. a=[3 2 1 0.5];b=[2 1 0.5 0.25];
  9. number=[10 10 100 20];
  10. x0=[];y0=[];xlat=[];ylon=[];
  11. [x0new1,y0new1,xlat1,ylon1]=myellipse_one(depth,a(1),b(1),number(1),x0,y0,xlat,ylon)
  12. x0=x0new1;y0=y0new1;xlat=xlat1;ylon=ylon1;
  13. for i=2:length(a)
  14.     [x0new2,y0new2,xlat2,ylon2]=myellipse_rest(depth,a(i),b(i),number(i),x0,y0,xlat,ylon)
  15.     x0=x0new2;y0=y0new2;xlat=xlat2;ylon=ylon2;
  16. end
  17. %  绘制椭圆图形
  18. [m,n]=size(xlat)
  19.     plot(xlat,ylon)
  20.     axis([0 60 0 60])
复制代码
  1. function [x0new,y0new,xlat,ylon]=myellipse_one(depth,a,b,number,x0,y0,xlat,ylon)
  2. xxlen=length(x0);yylen=length(y0);
  3. x0(xxlen+1)=a+(depth-2*a)*rand;
  4. y0(yylen+1)=a+(depth-2*a)*rand;
  5. ecc=axes2ecc(a,b);
  6. [lat1,lon1]=ellipse1(x0(xxlen+1),y0(yylen+1),[a,ecc],rand*180);
  7. xlat=[xlat,lat1];ylon=[ylon,lon1];
  8. for i=2:number
  9.     while 1
  10.     x0(xxlen+i)=a+(depth-2*a)*rand;
  11.     y0(yylen+i)=a+(depth-2*a)*rand;
  12.     [lati,loni]=ellipse1(x0(xxlen+i),y0(yylen+i),[a,ecc],rand*180);
  13.     xlat=[xlat,lati];ylon=[ylon,loni];

  14.     for j=1:xxlen+i-1
  15.         in1=[];in2=[];
  16.         in1=inpolygon(xlat(:,xxlen+i),ylon(:,yylen+i),xlat(:,j),ylon(:,j))
  17.         in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i))
  18.         if all(in1==0)&&all(in2==0)
  19.             n_j=j;
  20.             if n_j==xxlen+i-1
  21.                N=1
  22.                 break
  23.             end
  24.             continue
  25.         else
  26.               N=2;
  27.               xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
  28.               break
  29.          end
  30.      end
  31.       if N==1
  32.         break
  33.       end
  34.    end
  35. end
  36. x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;
复制代码
  1. function [x0new,y0new,xlat,ylon]=myellipse_rest(depth,a,b,number,x0,y0,xlat,ylon)
  2. xxlen=length(x0);yylen=length(y0);
  3. ecc=axes2ecc(a,b);N=0;
  4. for i=1:number
  5.      while 1
  6.     x0(xxlen+i)=a+(depth-2*a)*rand;
  7.     y0(xxlen+i)=a+(depth-2*a)*rand;
  8.     [lat,lon]=ellipse1(x0(xxlen+i),y0(xxlen+i),[a,ecc],rand*180);
  9.    xlat=[xlat,lat];ylon=[ylon,lon];
  10.     for j=1:xxlen+i-1
  11.         in1=[];in2=[];
  12.         in1=inpolygon(xlat(:,xxlen+i),ylon(:,xxlen+i),xlat(:,j),ylon(:,j));
  13.         in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i));
  14.         if all(in1==0)&&all(in2==0)
  15.              n_j=j;
  16.              if n_j==xxlen+i-1
  17.                 N=1;
  18.                 break
  19.             end
  20.            continue
  21.        else
  22.              N=2;
  23.              xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
  24.               break
  25.        end
  26.     end
  27.        if N==1
  28.        break
  29.        end
  30.   end
  31. end
  32. x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;
复制代码

评分

1

查看全部评分

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 10:39 , Processed in 0.075778 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表