0001 if nargin<4, sOption='All'; end
0002
0003 bxaxis=true; byaxis=true; bzaxis=true;
0004 if strcmpi(sOption,'x'), byaxis=false; bzaxis=false; end
0005 if strcmpi(sOption,'y'), bxaxis=false; bzaxis=false; end
0006 if strcmpi(sOption,'z'), bxaxis=false; byaxis=false; end
0007
0008 nNbG=length(a);
0009
0010 stEllGeom.a=a;
0011 stEllGeom.b=b;
0012 stEllGeom.c=c;
0013 stEllGeom.type='General';
0014
0015
0016 zeroa=0*a;
0017 stEllGeom.S=zeroa;
0018 if (bxaxis), stEllGeom.L1=zeroa;
0019 stEllGeom.exix2ave=zeroa; stEllGeom.exix4ave=zeroa; end
0020 if (byaxis), stEllGeom.L2=zeroa;
0021 stEllGeom.exiy2ave=zeroa; stEllGeom.exiy4ave=zeroa; end
0022 if (bzaxis), stEllGeom.L3=zeroa;
0023 stEllGeom.exiz2ave=zeroa; stEllGeom.exiz4ave=zeroa; end
0024 if strcmpi(sOption,'All')
0025 stEllGeom.exix2y2ave=zeroa;
0026 stEllGeom.exix2z2ave=zeroa;
0027 stEllGeom.exiy2z2ave=zeroa;
0028 end
0029
0030
0031 for indg=1:nNbG
0032
0033 abc=a(indg)*b(indg)*c(indg);
0034
0035 fsqrts=@(t,f) (b(indg)^2*c(indg)^2*sin(t).^2.*cos(f).^2 ...
0036 +a(indg)^2*c(indg)^2*sin(t).^2.*sin(f).^2+a(indg)^2*b(indg)^2*cos(t).^2).^(1/2);
0037
0038
0039 dfint = @(t,f) sin(t).*fsqrts(t,f);
0040 S=dblquad(dfint,0,pi,0,2*pi);
0041 stEllGeom.S(indg)=S;
0042
0043 fstovr3=@(t,f) sin(t)./(a(indg)^2*sin(t).^2.*cos(f).^2 + b(indg)^2*sin(t).^2.*sin(f).^2 ...
0044 + c(indg)^2*cos(t).^2).^(3/2);
0045
0046
0047 if (bxaxis)
0048
0049 fxova= @(t,f) sin(t).*cos(f);
0050
0051 dfint = @(t,f) abc/(4*pi)*fxova(t,f).^2 .* fstovr3(t,f);
0052 stEllGeom.L1(indg)=dblquad(dfint,0,pi,0,2*pi);
0053 end
0054 if (byaxis)
0055
0056 fyovb= @(t,f) sin(t).*sin(f);
0057 dfint = @(t,f) abc/(4*pi)*fyovb(t,f).^2 .* fstovr3(t,f);
0058 stEllGeom.L2(indg)=dblquad(dfint,0,pi,0,2*pi);
0059 end
0060 if (bzaxis)
0061
0062 fzovc= @(t,f) cos(t);
0063 dfint = @(t,f) abc/(4*pi)*fzovc(t,f).^2 .* fstovr3(t,f);
0064 stEllGeom.L3(indg)=dblquad(dfint,0,pi,0,2*pi);
0065 end
0066
0067 if ~strcmpi(sOption,'L')
0068
0069 if (bxaxis)
0070 dfint= @(t,f) (b(indg)*c(indg))^2 *fxova(t,f).^2 ./ fsqrts(t,f) .* sin(t);
0071 stEllGeom.exix2ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0072 dfint= @(t,f) (b(indg)*c(indg))^4 * fxova(t,f).^4 ./ fsqrts(t,f).^3 .* sin(t);
0073 stEllGeom.exix4ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0074 end
0075 if (byaxis)
0076 dfint= @(t,f) (a(indg)*c(indg))^2 *fyovb(t,f).^2 ./ fsqrts(t,f) .* sin(t);
0077 stEllGeom.exiy2ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0078 dfint= @(t,f) (a(indg)*c(indg))^4 * fyovb(t,f).^4 ./ fsqrts(t,f).^3 .* sin(t);
0079 stEllGeom.exiy4ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0080 end
0081 if (bzaxis)
0082 if (bxaxis && byaxis)
0083 stEllGeom.exiz2ave(indg)=1-stEllGeom.exix2ave(indg)-stEllGeom.exiy2ave(indg);
0084 else
0085
0086 dfint= @(t,f) (a(indg)*b(indg))^2 *fzovc(t,f).^2 ./ fsqrts(t,f) .* sin(t);
0087 stEllGeom.exiz2ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0088 end
0089 dfint= @(t,f) (a(indg)*b(indg))^4 * fzovc(t,f).^4 ./ fsqrts(t,f).^3 .* sin(t);
0090 stEllGeom.exiz4ave(indg)=dblquad(dfint,0,pi,0,2*pi)/S;
0091 end
0092
0093 if strcmpi(sOption,'All')
0094 stEllGeom.exix2y2ave(indg)=(1-2*stEllGeom.exiz2ave(indg)+stEllGeom.exiz4ave(indg) ...
0095 -stEllGeom.exix4ave(indg)-stEllGeom.exiy4ave(indg))/2;
0096 stEllGeom.exix2z2ave(indg)=(1-2*stEllGeom.exiy2ave(indg)+stEllGeom.exiy4ave(indg) ...
0097 -stEllGeom.exix4ave(indg)-stEllGeom.exiz4ave(indg))/2;
0098 stEllGeom.exiy2z2ave(indg)=(1-2*stEllGeom.exix2ave(indg)+stEllGeom.exix4ave(indg) ...
0099 -stEllGeom.exiy4ave(indg)-stEllGeom.exiz4ave(indg))/2;
0100 end
0101
0102
0103
0104
0105
0106
0107
0108
0109 end
0110 end
0111
0112
0113
0114