Description of vb_solid_angle

0001 function [D, s0] = vb_solid_angle_grk(V,F,xx)
0002 % solid angle calculation by Gararkin method
0003 %   [D, s0] = vb_solid_angle_grk(V,F,xx)
0004 % 境界頂点係数行列の計算
0005 % (1次近似ポテンシャル)
0006 %
0007 % D     : 境界頂点係数行列 (重み付き立体角)
0008 % s0     : 三角面の面積法線
0009 %
0010 % V     : 境界三角面の座標
0011 % F     : 境界三角面頂点インデックス
0012 % xx    : 境界三角面の法線ベクトル
0013 %
0014 %
0015 % 2004-02-06 Taku Yoshioka
0016 % 2004-12-26 M. Sato modified
0017 % 2007-12-26 M. Sato modified
0018 %
0019 % Copyright (C) 2011, ATR All Rights Reserved.
0020 % License : New BSD License(see VBMEG_LICENSE.txt)
0021 
0022 Nskip = 100;
0023 NSKIP = 1000;
0024 
0025 %fprintf('--- Solid angle calculation\n');
0026 % 三角面の数　Number of triangle patch
0027 NF = size(F,1);
0028 NV = size(V,1);
0029 
0030 % 境界面係数行列
0031 D = zeros(NV,NV);
0032 %A = zeros(NV,NV);
0033 
0034 % 三角面の頂点　triangle vertex
0035 x1 = V(F(:,1),:);
0036 x2 = V(F(:,2),:);
0037 x3 = V(F(:,3),:);
0038 
0039 % 三角面の面積法線計算　(外積)
0040 z1 = x2-x1;
0041 z2 = x3-x1;
0042 s0 = vb_cross2(z1,z2);
0043 S  = sum(s0.^2,2);
0044 Sinv2 = 1./S;
0045 S = sqrt(S);
0046 
0047 % 法線を外向きにする
0048 xxS = (xx(F(:,1),:)+xx(F(:,2),:)+xx(F(:,3),:))./3;
0049 ix  = find( sum(s0.*xxS,2) < 0 );
0050 % 頂点入れ替え
0051 tmp      = x2(ix,:);
0052 x2(ix,:) = x3(ix,:);
0053 x3(ix,:) = tmp;
0054 % インデックス入れ替え
0055 tmp      = F(ix,2);
0056 F(ix,2)  = F(ix,3);
0057 F(ix,3)  = tmp;
0058 
0059 s0(ix,:) = -1*s0(ix,:);
0060 
0061 %
0062 % 各面の3つの内点について、それらの点から見た他の面の
0063 % (重み付き)立体角を計算し、その面の頂点でインデックスが与えられる
0064 % 3つの基底関数への寄与を計算する。
0065 %
0066 
0067 % 境界三角面(基準三角面)に関するループ
0068 for n = 1:NF
0069   % 内点
0070   ip = [(4*x1(n,:) +   x2(n,:) +   x3(n,:))/6; ...
0071         (  x1(n,:) + 4*x2(n,:) +   x3(n,:))/6; ...
0072         (  x1(n,:) +   x2(n,:) + 4*x3(n,:))/6];
0073   
0074   % 自分自身のインデックスを除く
0075   ix = [1:(n-1), (n+1):NF];
0076   
0077   Fn = F(ix,:);    % 対象三角面の頂点インデックス
0078   s  = s0(ix,:);% 三角面法線
0079   y1 = x1(ix,:);% 対象三角面の3つの頂点
0080   y2 = x2(ix,:);% 対象三角面の3つの頂点
0081   y3 = x3(ix,:);% 対象三角面の3つの頂点
0082 
0083   Sinv = Sinv2(ix);% 1/三角面の面積
0084   
0085   % 重み付き立体角の計算
0086   D0 = zeros(3,NV);
0087 
0088   % ３つの内点に関するループ
0089   for i = 1:3
0090     % x0 : 基準点(基準三角面の一つの内点)
0091     x0 = ip(i,:);
0092     y0 = x0(ones(NF-1,1),:);
0093     
0094     % 一つの内点から見た他の三角面(対象三角面)の重み付き立体角
0095     [omega1,omega2,omega3] = ...
0096             vb_solid_other(y0,y1,y2,y3,s,Sinv);
0097     
0098     % 対象三角面の３つの頂点に立体角寄与を分配
0099     for j = 1:NF-1
0100       D0(i,Fn(j,1)) = D0(i,Fn(j,1))+omega1(j);
0101       D0(i,Fn(j,2)) = D0(i,Fn(j,2))+omega2(j);
0102       D0(i,Fn(j,3)) = D0(i,Fn(j,3))+omega3(j);
0103     end
0104   end
0105   
0106   % 重み付き立体角の数値積分(3つの内点に関する重み付き和)
0107   % Sは面積の2倍なので、3*6*2=36で割っている
0108   D(F(n,1),:) = D(F(n,1),:)+S(n)/36*(4*D0(1,:)+D0(2,:)+D0(3,:));
0109   D(F(n,2),:) = D(F(n,2),:)+S(n)/36*(D0(1,:)+4*D0(2,:)+D0(3,:));
0110   D(F(n,3),:) = D(F(n,3),:)+S(n)/36*(D0(1,:)+D0(2,:)+4*D0(3,:));
0111   
0112   if mod(n,Nskip)==0, fprintf('-'); end;
0113   if mod(n,NSKIP)==0, fprintf(' %3.0f %% done\n', 100*n/NF); end;
0114 end
0115 
0116 % 立体角を(球の立体角/2)=２πで割る
0117 D   = (-1/(2*pi))*D;
0118 
0119 return
vb_solid_angle_grk

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
