vb_solid_other

 solid angle calculation by Gararkin method
 3 solid angles from 3 points in one triangle are calculated
  [omega1,omega2,omega3] = vb_solid_other(x0,x1,x2,x3,s,Sinv)
 基準点から見た他の三角面の重み付き立体角計算

 x0         : 基準点(基準三角面の一つの内点)
 (x1,x2,x3) : 対象三角面の3つの頂点
 s          : 対象三角面法線
 Sinv       : 1/対象三角面の面積

 2004-02-06 Taku Yoshioka
 2004-12-26 M. Sato modified

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [omega1,omega2,omega3]=...
0002                vb_solid_other(x0,x1,x2,x3,s,Sinv)
0003 % solid angle calculation by Gararkin method
0004 % 3 solid angles from 3 points in one triangle are calculated
0005 %  [omega1,omega2,omega3] = vb_solid_other(x0,x1,x2,x3,s,Sinv)
0006 % 基準点から見た他の三角面の重み付き立体角計算
0007 %
0008 % x0         : 基準点(基準三角面の一つの内点)
0009 % (x1,x2,x3) : 対象三角面の3つの頂点
0010 % s          : 対象三角面法線
0011 % Sinv       : 1/対象三角面の面積
0012 %
0013 % 2004-02-06 Taku Yoshioka
0014 % 2004-12-26 M. Sato modified
0015 %
0016 % Copyright (C) 2011, ATR All Rights Reserved.
0017 % License : New BSD License(see VBMEG_LICENSE.txt)
0018 
0019 MIN_VAL = realmin*100;
0020 
0021 y1 = x1-x0;
0022 y2 = x2-x0;
0023 y3 = x3-x0;
0024 
0025 w1 = y2-y3;
0026 w2 = y3-y1;
0027 w3 = y1-y2;
0028 
0029 z1 = vb_cross2(y2,y3);
0030 z2 = vb_cross2(y3,y1);
0031 z3 = vb_cross2(y1,y2);
0032 
0033 d  = sum(z1.*y1,2);
0034 
0035 W1 = sqrt(sum(w1.^2,2));
0036 W2 = sqrt(sum(w2.^2,2));
0037 W3 = sqrt(sum(w3.^2,2));
0038 
0039 Y1 = sqrt(sum(y1.^2,2));
0040 Y2 = sqrt(sum(y2.^2,2));
0041 Y3 = sqrt(sum(y3.^2,2));
0042 
0043 YA = Y1.*Y2.*Y3;
0044 YB = Y1.*sum(y2.*y3,2)+Y2.*sum(y3.*y1,2)+Y3.*sum(y1.*y2,2);
0045 omega0 = 2*atan(d./(YA+YB));
0046 %
0047 % gamma1
0048 %
0049 tmp1 = W3.*Y1-sum(w3.*y1,2);
0050 tmp2 = W3.*Y2-sum(w3.*y2,2);
0051 
0052 %gamma1 = log(tmp2./tmp1)./W3;
0053 ix1 = find( tmp1 >  MIN_VAL & tmp2 >  MIN_VAL);
0054 ix2 = find( tmp1 <= MIN_VAL | tmp2 <= MIN_VAL);
0055 gamma1 = zeros(length(tmp1),1);
0056 gamma1(ix1) = log(tmp2(ix1)./tmp1(ix1))./W3(ix1);
0057 gamma1(ix2) = log(Y2(ix2)./Y1(ix2));
0058 
0059 %
0060 % gamma2
0061 %
0062 tmp1 = W1.*Y2-sum(w1.*y2,2);
0063 tmp2 = W1.*Y3-sum(w1.*y3,2);
0064 %gamma2 = log(tmp2./tmp1)./W1;
0065 ix1 = find( tmp1 >  MIN_VAL & tmp2 >  MIN_VAL);
0066 ix2 = find( tmp1 <= MIN_VAL | tmp2 <= MIN_VAL);
0067 gamma2 = zeros(length(tmp1),1);
0068 gamma2(ix1) = log(tmp2(ix1)./tmp1(ix1))./W1(ix1);
0069 gamma2(ix2) = log(Y3(ix2)./Y2(ix2));
0070 
0071 %
0072 % gamma3
0073 %
0074 tmp1 = W2.*Y3-sum(w2.*y3,2);
0075 tmp2 = W2.*Y1-sum(w2.*y1,2);
0076 %gamma3 = log(tmp2./tmp1)./W2;
0077 ix1 = find( tmp1 >  MIN_VAL & tmp2 >  MIN_VAL);
0078 ix2 = find( tmp1 <= MIN_VAL | tmp2 <= MIN_VAL);
0079 gamma3 = zeros(length(tmp1),1);
0080 gamma3(ix1) = log(tmp2(ix1)./tmp1(ix1))./W2(ix1);
0081 gamma3(ix2) = log(Y1(ix2)./Y3(ix2));
0082 
0083 omega4 = gamma2(:,ones(1,3)).*w1 ...
0084        + gamma3(:,ones(1,3)).*w2 ...
0085        + gamma1(:,ones(1,3)).*w3;
0086 
0087 omega1 = Sinv.*( sum(s.*z1,2).*omega0 ...
0088               + d.*(sum(w1.*omega4,2)) );
0089 omega2 = Sinv.*( sum(s.*z2,2).*omega0 ...
0090               + d.*(sum(w2.*omega4,2)) );
0091 omega3 = Sinv.*( sum(s.*z3,2) .*omega0 ...
0092               + d.*(sum(w3.*omega4,2)) );