V = analytic_2d_circle(J, [s_h, s_i, b, a, angl])
Voltage around a 2D circle with properties
conductivities: tank(s_h), inclusion(s_i)
tank radius = 1
inclusion: radius(a), distance(b) at angle(angl)
J = current density at a regular sample of nodes on the boundary
because analytic_2d_circle works on current density, it needs
to be scaled with the number of sample points, to get I
Based on eqn 21 in Seagar and Bates (1985).
(C) Andy Adler 2009. License: GPL V2 or V3
$Id: analytic_2d_circle.html 2819 2011-09-07 16:43:11Z aadler $
USAGE EXAMPLE: compare FEM to analytic model
imdl= mk_common_model('c2c0',16);
img = calc_jacobian_bkgnd(imdl);
img.fwd_solve.get_all_meas = 1;
img.elem_data(1:256)= 0.1;
vi= fwd_solve(img);
[jnk,maxl] = min(vi.volt(:,1));
vi= vi.volt(maxl:end,1);
lv= length(vi);
vol = get_elem_volume(img.fwd_model);
vt= sum(vol); va= sum(vol(img.elem_data ~=1));
rad = sqrt(va/vt);
I = zeros(lv,1); I(1+[0,lv/16]) = [-25,25]*lv/16;
vsi= analytic_2d_circle(I, [1, 0.1, 0.0, rad, 0]);
plot([vi,vsi]);0001 function V = analytic_2d_circle(I, params) 0002 % V = analytic_2d_circle(J, [s_h, s_i, b, a, angl]) 0003 % Voltage around a 2D circle with properties 0004 % conductivities: tank(s_h), inclusion(s_i) 0005 % tank radius = 1 0006 % inclusion: radius(a), distance(b) at angle(angl) 0007 % 0008 % J = current density at a regular sample of nodes on the boundary 0009 % because analytic_2d_circle works on current density, it needs 0010 % to be scaled with the number of sample points, to get I 0011 % 0012 % Based on eqn 21 in Seagar and Bates (1985). 0013 % 0014 % (C) Andy Adler 2009. License: GPL V2 or V3 0015 % $Id: analytic_2d_circle.html 2819 2011-09-07 16:43:11Z aadler $ 0016 % 0017 % USAGE EXAMPLE: compare FEM to analytic model 0018 % 0019 % imdl= mk_common_model('c2c0',16); 0020 % img = calc_jacobian_bkgnd(imdl); 0021 % img.fwd_solve.get_all_meas = 1; 0022 % img.elem_data(1:256)= 0.1; 0023 % vi= fwd_solve(img); 0024 % 0025 % [jnk,maxl] = min(vi.volt(:,1)); 0026 % vi= vi.volt(maxl:end,1); 0027 % lv= length(vi); 0028 % 0029 % vol = get_elem_volume(img.fwd_model); 0030 % vt= sum(vol); va= sum(vol(img.elem_data ~=1)); 0031 % rad = sqrt(va/vt); 0032 % 0033 % I = zeros(lv,1); I(1+[0,lv/16]) = [-25,25]*lv/16; 0034 % vsi= analytic_2d_circle(I, [1, 0.1, 0.0, rad, 0]); 0035 % 0036 % plot([vi,vsi]); 0037 0038 0039 if ~all(abs(sum(I,1)) <1e-12); error('net I must be 0'); end 0040 [ll,nv] = size(I); 0041 0042 IF= fft(I); IF(1)= []; % Don't include zero term 0043 [SM,ll2] = solve_matrix(ll, params); 0044 VF2= SM*IF(ll2); 0045 VF= zeros(size(I)); 0046 VF(1+[ll2, ll-ll2],:)= [VF2;conj(VF2)]; 0047 V = ifft(VF); 0048 0049 if norm(imag(V))>1e-9; error('Unexpected Imaginary output - probably a bug'); end 0050 V= real(V); 0051 0052 0053 function [SM,ll2] = solve_matrix(ll, params) 0054 ll1= ll-1; % IGNORE DC 0055 ll2 = 1:ceil(ll1/2); 0056 0057 SM = eidors_obj('get-cache', {ll,params}, 'analytic_2d_circle'); 0058 if ~isempty(SM) 0059 eidors_msg('analytic_2d_circle: using cached value', 3); 0060 return 0061 end 0062 0063 s_h = params(1); 0064 s_i = params(2); 0065 b = params(3); 0066 a = params(4); 0067 alp = params(5); 0068 0069 mll2= max(ll2); 0070 K = ll2(:); 0071 D = spdiags(K,0,mll2,mll2); 0072 iD = spdiags(1./K,0,mll2,mll2); % inv D 0073 0074 if s_h == s_i 0075 SM= 1/s_h * iD; 0076 return 0077 end 0078 0079 mu = (s_h - s_i) / (s_h + s_i); 0080 0081 if b==0 0082 mua2k = mu*a.^(2*K); 0083 K = K.*(1 - mua2k)./(1 + mua2k); 0084 iD = spdiags(1./K,0,mll2,mll2); % inv(D) 0085 0086 SM= 1/s_h * iD; 0087 return 0088 end 0089 0090 0091 0092 T = sparse(mll2,mll2); 0093 for m=ll2; for n=ll2 0094 sm= 0; 0095 for p=1:min(m,n) 0096 % sm=sm+ nchoosek(m-1,p-1)*nchoosek(n,p)*ab2^p; 0097 % stop calculating nchoosek when too small, because 0098 dsm = nchoosek(m-1,p-1)*nchoosek(n,p)*a^(2*p)*b^(m+n-2*p); 0099 % disp([m,n,p,log10(dsm)]) 0100 if dsm < 1e-11; break; end 0101 sm=sm + dsm; 0102 end 0103 % T(i,j) = mu*(sm/m)* b^(m+n) * exp(1j*(m-n)*alp); 0104 T(m,n) = mu*(sm/m)* exp(1j*(m-n)*alp); 0105 end; end 0106 0107 I= speye(mll2); 0108 SM = (I-T*D)\(I+T*D)*iD; 0109 0110 eidors_obj('set-cache', {ll,params}, 'analytic_2d_circle', SM); 0111 eidors_msg('analytic_2d_circle: setting cached value', 3);