% 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