equilvol_mc_jd

0001 function pV = equilvol_mc_jd(Rp,theta,ptx,pty,N,mfac,method)
0002 %
0003 if nargin < 7,
0004     method = 'linear';
0005 end
0006 %
0007 if nargin < 6,
0008     mfac = 50;
0009 end
0010 %
0011 if nargin < 5,
0012     N = 10;
0013 end
0014 %
0015 [np,nt] = size(ptx);
0016 %
0017 if any(theta ~= linspace(0,2*pi,nt)),
0018     error('case of nonuniform theta grid not yet implemented')
0019 end
0020 %
0021 pmx = interpft(ptx(:,1:nt-1),(nt-1)*mfac,2); %FFT interpolation
0022 pmy = interpft(pty(:,1:nt-1),(nt-1)*mfac,2); %FFT interpolation
0023 %
0024 pmx = [pmx,pmx(:,1)]; %2*pi point retrieved
0025 pmy = [pmy,pmy(:,1)]; %2*pi point retrieved
0026 pmr = sqrt(pmx.^2 + pmy.^2);
0027 mtheta = linspace(0,2*pi,mfac*(nt-1)+1);
0028 %
0029 xmax = max(pmx(end,:),[],2);
0030 ymin = min(pmy(end,:),[],2);
0031 ymax = max(pmy(end,:),[],2);
0032 %
0033 Rmax = Rp + xmax;
0034 %
0035 Xmin = -Rmax;
0036 Xmax = Rmax;
0037 Ymin = -Rmax;
0038 Ymax = Rmax;
0039 Zmin = ymin;
0040 Zmax = ymax;
0041 %
0042 Vtot = (Xmax - Xmin)*(Ymax - Ymin)*(Zmax - Zmin);
0043 %
0044 F = rand([N,3]);
0045 %
0046 X = Xmin + F(:,1)*(Xmax - Xmin);
0047 Y = Ymin + F(:,2)*(Ymax - Ymin);
0048 Z = Zmin + F(:,3)*(Zmax - Zmin);
0049 R = sqrt(X.^2 + Y.^2);
0050 %
0051 x = R - Rp;
0052 y = Z;
0053 r = sqrt(x.^2 + y.^2);
0054 theta = atan(y./x) + pi*(x < 0) + 2*pi*(x > 0 & y < 0);
0055 theta(x == 0) = pi/2 + pi*(y(x == 0) < 0);
0056 %
0057 pr = interp1(mtheta,pmr.',theta,method).';
0058 %
0059 pV = Vtot*sum(repmat(r.',[np,1]) < pr,2)/N;
0060 
0061 
0062 
0063 
0064 
0065 
0066 
0067 
0068 
0069 
0070 
0071