% 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