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.m 4898 2015-05-06 14:26:36Z 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.m 4898 2015-05-06 14:26:36Z 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 % TODO: 0039 % (aa May 2015) Works well for 16 electrodes, but not for 8,32 0040 0041 if ischar(I) && strcmp(I, 'UNIT_TEST'), do_unit_test; return; end 0042 0043 if ~all(abs(sum(I,1)) <1e-12); error('net I must be 0'); end 0044 [ll,nv] = size(I); 0045 0046 IF= fft(I); IF(1)= []; % Don't include zero term 0047 [SM,ll2] = eidors_cache(@solve_matrix,{ll, params},'analytic_2d_circle'); 0048 VF2= SM*IF(ll2); 0049 VF= zeros(size(I)); 0050 VF(1+[ll2, ll-ll2],:)= [VF2;conj(VF2)]; 0051 V = ifft(VF); 0052 0053 if norm(imag(V))>1e-9; error('Unexpected Imaginary output - probably a bug'); end 0054 V= real(V); 0055 0056 0057 function [SM,ll2] = solve_matrix(ll, params) 0058 ll1= ll-1; % IGNORE DC 0059 ll2 = 1:ceil(ll1/2); 0060 0061 s_h = params(1); 0062 s_i = params(2); 0063 b = params(3); 0064 a = params(4); 0065 alp = params(5); 0066 0067 mll2= max(ll2); 0068 K = ll2(:); 0069 D = spdiags(K,0,mll2,mll2); 0070 iD = spdiags(1./K,0,mll2,mll2); % inv D 0071 0072 if s_h == s_i 0073 SM= 1/s_h * iD; 0074 return 0075 end 0076 0077 mu = (s_h - s_i) / (s_h + s_i); 0078 0079 if b==0 0080 mua2k = mu*a.^(2*K); 0081 K = K.*(1 - mua2k)./(1 + mua2k); 0082 iD = spdiags(1./K,0,mll2,mll2); % inv(D) 0083 0084 SM= 1/s_h * iD; 0085 return 0086 end 0087 0088 0089 0090 T = sparse(mll2,mll2); 0091 for m=ll2; for n=ll2 0092 sm= 0; 0093 for p=1:min(m,n) 0094 % sm=sm+ nchoosek(m-1,p-1)*nchoosek(n,p)*ab2^p; 0095 % stop calculating nchoosek when too small, because 0096 dsm = nchoosek(m-1,p-1)*nchoosek(n,p)*a^(2*p)*b^(m+n-2*p); 0097 % disp([m,n,p,log10(dsm)]) 0098 if dsm < 1e-11; break; end 0099 sm=sm + dsm; 0100 end 0101 % T(i,j) = mu*(sm/m)* b^(m+n) * exp(1j*(m-n)*alp); 0102 T(m,n) = mu*(sm/m)* exp(1j*(m-n)*alp); 0103 end; end 0104 0105 I= speye(mll2); 0106 SM = (I-T*D)\(I+T*D)*iD; 0107 0108 function do_unit_test 0109 0110 N_elec = 16; 0111 imdl= mk_common_model('d2c0',N_elec); 0112 img = mk_image(imdl); 0113 img.fwd_solve.get_all_meas = 1; 0114 img.elem_data(1:256)= 0.1; 0115 vi= fwd_solve(img); 0116 0117 [jnk,maxl] = min(vi.volt(:,1)); 0118 vi= vi.volt(maxl:end,1); 0119 lv= length(vi); 0120 0121 vol = get_elem_volume(img.fwd_model); 0122 vt= sum(vol); va= sum(vol(img.elem_data ~=1)); 0123 rad = sqrt(va/vt); 0124 0125 I = zeros(lv,1); I(1+[0,lv/N_elec]) = [-25,25]*lv/N_elec; 0126 vsi= analytic_2d_circle(I, [1, 0.1, 0.0, rad, 0]); 0127 0128 plot([vi,vsi]); 0129 unit_test_cmp('16 elec compare', vi, vsi, 1);