vb_make_sphere_point

 Make points on a hemisphere
 Pos で指定された球面座標の範囲に座標点を作成する

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

0001 function [X ,NDlist ,Jlist ] = vb_make_sphere_point(Pos)
0002 % Make points on a hemisphere
0003 % Pos で指定された球面座標の範囲に座標点を作成する
0004 %
0005 % Copyright (C) 2011, ATR All Rights Reserved.
0006 % License : New BSD License(see VBMEG_LICENSE.txt)
0007 
0008 % Pos の構造体
0009 %   theta  : 仰角    : [theta_min theta_max] (0:水平 pi/2:真上)
0010 %   phi    : 方位角  : [  phi_min   phi_max] (0:pi)
0011 %   r      : Radius
0012 %   Ntheta : Number of divisions in theta
0013 %   Nphi   : Number of divisions in phi
0014 
0015 % X            = [rx, vphi, vtheta] : (N x 9) matrix
0016 % rx           = (theta,phi,r)で指定される座標ベクトル
0017 % vphi, vtheta : rx に垂直な２つの接線ベクトル
0018 % Jlist{n} : theta(n) に対応する点のインデックス
0019 % NDlist(n) : Jlist{n} の点数
0020 
0021 theta_step = (Pos.theta_max - Pos.theta_min);
0022 phi_range  = (Pos.phi_max   - Pos.phi_min)/2;
0023 
0024 if theta_step > 0,
0025     theta_step = theta_step/Pos.Ntheta;
0026 end;
0027 
0028 Nphi   = Pos.Nphi;
0029 
0030 theta  = theta_step*[0:Pos.Ntheta]+Pos.theta_min;
0031 Ntheta = length(theta);
0032 NDlist = zeros(Ntheta,1);
0033 
0034 Phi    = [];
0035 Theta  = [];
0036 Nstart = 1;
0037 Jlist  = cell(Ntheta,1);
0038 
0039 for n=1:Ntheta,
0040     DNphi=fix(Nphi*cos(theta(n))/2);
0041     
0042     if DNphi> 0,
0043         if DNphi <2,
0044             DNphi=2;
0045         end;
0046         dphi=phi_range/DNphi;
0047     else
0048         dphi=0;
0049     end;
0050     
0051     if (phi_range < pi) | ( DNphi == 0 ),
0052         jj          = 0:(2*DNphi);
0053         phi       = jj*dphi + Pos.phi_min;
0054         NDlist(n) = 2*DNphi+1;
0055         Jlist{n}  = jj + Nstart;
0056         Nstart      = Nstart + 2*DNphi+1;
0057     else
0058         jj          = 0:(2*DNphi-1);
0059         phi       = jj*dphi + Pos.phi_min;
0060         NDlist(n) = 2*DNphi;
0061         Jlist{n}  = jj + Nstart;
0062         Nstart      = Nstart + 2*DNphi;
0063     end;
0064     
0065     Phi    = [Phi, phi];
0066     Theta  = [Theta, theta(n)*ones(1,length(phi)) ];
0067 end;
0068 
0069 N=length(Theta);
0070 
0071 r      = Pos.r;
0072 
0073 [x ,y ,z ] = sph2cart(Phi,Theta,r); 
0074 [x1,y1,z1] = sph2cart(Phi+0.5*pi, zeros(1,length(Theta)), 1); 
0075 [x2,y2,z2] = sph2cart(Phi, Theta+0.5*pi, 1); 
0076 
0077 X          = [x(:),y(:),z(:), x1(:),y1(:),z1(:), x2(:),y2(:),z2(:)];
0078