% 6.581 PS3 Qsn2
% computeReactionPotential.m
% calculate reaction pontential at charge locations
% of a spherical molecule centered at origin
%infile = 'sphere48.qif';
infile = 'sphere192.qif';
%infile = 'sphere768.qif';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% YOU ONLY NEED TO MODIFY THIS SECTION
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% radius of sphere
radius = 1;
% charge coordinates and magnitudes
qMag = [1 1 1 1]'; % qMag = [q1 q2 q3 ... qN]'
qCoord = [ % qCoord = [q1x q1y q1z; q2x q2y q2z; ...; qNx qNy qNz]
0.4 -0.2 0;
0.4 0.2 0;
0.8 0.2 0;
0.8 -0.2 0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% END OF SECTION YOU NEED TO MODIFY
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = length(qMag);
rPot = zeros(N,1); % reaction potential at qCoord will be stored here
% Permitivity of dielectric regions
E_0 = 8.854187818E-12;
E_in = 4*E_0;
E_out = 80*E_0;
E_rel = E_out/E_in;
% Read in the panels for a unit sphere
[panels] = readpanels(infile);
% Adjust panels for appropriate size of sphere
numpanels = size(panels,3);
for i=1:numpanels
numverts = panels(1,1,i);
panel = panels(2:numverts+1,:,i);
[TH,PHI,R] = cart2sph(panel(:,1),panel(:,2),panel(:,3));
[X,Y,Z] = sph2cart(TH,PHI,radius);
panels(2:numverts+1,1,i) = X;
panels(2:numverts+1,2,i) = Y;
panels(2:numverts+1,3,i) = Z;
end
done = 'read input file'
% Compute the panel centroids and areas
[centroids, normals, areas] = gencolloc(panels);
for i = 1:numpanels
% make sure normal is pointing outward, assuming center of sphere is at origin
if dot(centroids(i,:),normals(i,:)) < 0
normals(i,:) = -normals(i,:);
end
end
done = 'generated collocation points'
% Generate the collocation matrix
A = zeros(numpanels);
for i=1:numpanels
numverts = panels(1,1,i);
panel = panels(2:numverts+1,:,i);
[area, centroid, normal, p, dip, dpdn] = calcp(panel, centroids, normals);
A(:,i) = (1 - E_rel) / (4*pi) * dpdn';
A(i,i) = (1 + E_rel) / 2;
end
done = 'generated matrix'
% Create the rhs
b = zeros(numpanels,1);
for k = 1:N
for i = 1:numpanels
r = centroids(i,:) - qCoord(k,:);
r2 = sum(r.^2); % r2 = |r|^2
b(i) = b(i) + (1 - E_rel) / (4*pi) * qMag(k) * dot(r/norm(r),normals(i,:)) / r2;
end
end
% Solve for the charge density vector
sigma = A \ b;
done = 'solved for charge'
% visualize surface charge density
for i = 1:numpanels
numverts = panels(1,1,i);
panel = panels(2:numverts+1,:,i);
fill3(panel(:,1),panel(:,2),panel(:,3),sigma(i))
hold on
end
plot3(qCoord(:,1),qCoord(:,2),qCoord(:,3),'b.','MarkerSize',20)
hold off
colorbar
shading flat
alpha(0.5)
axis off
% Integrate the charge over the surface to compute the reaction potential
for k = 1:N
for i = 1:numpanels
numverts = panels(1,1,i);
panel = panels(2:numverts+1,:,i);
[area, centroid, normal, p] = calcp(panel, qCoord(k,:));
rPot(k) = rPot(k) + sigma(i)*p;
end
end
done = 'reaction potential computed'
rPot = rPot/(4*pi) % normalized by E_in