source_monopole

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 monopole source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 % monopole source
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 
0006 
0007 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009 
0010 % a conductive medium
0011 sigma = 0.2;
0012 
0013 % with a point source current magnitude
0014 Io = 10;
0015 
0016 % located at the origin
0017 IoR = [0,0,0];
0018 
0019 % consider a spherical surface
0020 sphereR = 1;
0021 FV = sphere_tri('ico',4,sphereR);
0022 
0023 % By symmetry, the current flow is radially outward.
0024 % At a distance r > 0 from the source, we have
0025 
0026 % radial location, vector coordinates
0027 R = FV.vertices;
0028 R(:,1) = R(:,1) - IoR(1);
0029 R(:,2) = R(:,2) - IoR(2);
0030 R(:,3) = R(:,3) - IoR(3);
0031 r = vector_magnitude(R,IoR);
0032 r = repmat(r,1,3);
0033 
0034 % extracellular currents at surface
0035 Je = [ Io / (4 * pi) ] * [ R ./ (r.^3) ];
0036 
0037 % electric potential at surface
0038 V  = [ Io / (4 * pi * sigma) ] * [ 1 ./ r(:,1) ];
0039 
0040 FV.Cdata = V;
0041 
0042 figure;
0043 H = patch('faces',FV.faces,'vertices',FV.vertices,...
0044   'FaceVertexCData',FV.Cdata,'facecolor','interp',...
0045   'edgecolor','none');
0046 
0047 
0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 
0051 clear FV
0052 
0053 % a conductive medium
0054 sigma = 0.2;
0055 
0056 % with a point source current magnitude
0057 Io = 10;
0058 
0059 % located at the origin
0060 IoR = [0,0,0];
0061 
0062 FV.vertices = ones(10*10,3);
0063 x = linspace(-10,10,10);
0064 y = x;
0065 n = 0;
0066 for i = 1:length(x),
0067   for j = 1:length(y),
0068     n = n + 1;
0069     FV.vertices(n,:) = [ x(i), y(j), 10 ];
0070   end
0071 end
0072 
0073 FV.faces = delaunay(FV.vertices(:,1),FV.vertices(:,2));
0074 
0075 % radial location, vector coordinates
0076 R = FV.vertices;
0077 R(:,1) = R(:,1) - IoR(1);
0078 R(:,2) = R(:,2) - IoR(2);
0079 R(:,3) = R(:,3) - IoR(3);
0080 r = vector_magnitude(R,IoR);
0081 r = repmat(r,1,3);
0082 
0083 % extracellular currents at surface
0084 Je = [ Io / (4 * pi) ] * [ R ./ (r.^3) ];
0085 
0086 % electric potential at surface
0087 V  = [ Io / (4 * pi * sigma) ] * [ 1 ./ r(:,1) ];
0088 
0089 FV.Cdata = V;
0090 
0091 figure;
0092 H = patch('faces',FV.faces,'vertices',FV.vertices,...
0093   'FaceVertexCData',FV.Cdata,'facecolor','interp',...
0094   'edgecolor','none');
0095 colormap(hot); colorbar
0096 
0097 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0099 
0100 clear FV
0101 
0102 % a conductive medium
0103 sigma = 0.2;
0104 
0105 % with a point source current magnitude
0106 Io = 10;
0107 
0108 % located at the origin
0109 IoR = [0,0,0];
0110 
0111 dim = 20;
0112 
0113 FV.vertices = ones(dim*dim*dim,3);
0114 x = linspace(-10,10,dim);
0115 y = x;
0116 z = x;
0117 n = 0;
0118 for i = 1:length(x),
0119   for j = 1:length(y),
0120     for k = 1:length(z),
0121       n = n + 1;
0122       FV.vertices(n,:) = [ x(i), y(j), z(k) ];
0123     end
0124   end
0125 end
0126 
0127 FV.faces = delaunay3(FV.vertices(:,1),FV.vertices(:,2),FV.vertices(:,3));
0128 
0129 % radial location, vector coordinates
0130 R = FV.vertices;
0131 R(:,1) = R(:,1) - IoR(1);
0132 R(:,2) = R(:,2) - IoR(2);
0133 R(:,3) = R(:,3) - IoR(3);
0134 r = vector_magnitude(R,IoR);
0135 r = repmat(r,1,3);
0136 
0137 % extracellular currents at surface
0138 Je = [ Io / (4 * pi) ] * [ R ./ (r.^3) ];
0139 
0140 % electric potential at surface
0141 V  = [ Io / (4 * pi * sigma) ] * [ 1 ./ r(:,1) ];
0142 
0143 
0144 % create Analyze volume of potential
0145 avw = avw_hdr_make;
0146 avw.hdr.dime.dim(2:4) = [dim,dim,dim];
0147 
0148 avw.img = zeros(dim,dim,dim);
0149 n = 0;
0150 for i = 1:length(x),
0151   for j = 1:length(y),
0152     for k = 1:length(z),
0153       n = n + 1;
0154       avw.img(i,j,k) = V(n);
0155     end
0156   end
0157 end
0158 avw_view(avw);