analytic_2d_circle

PURPOSE ^

V = analytic_2d_circle(J, [s_h, s_i, b, a, angl])

SYNOPSIS ^

function V = analytic_2d_circle(I, params)

DESCRIPTION ^

 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]);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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);

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005