%Mason Averill
%Spring 2019, PHYS 312

%Computational Project, flux through a cube with a point charge
%anywhere inside of it


%Note, this program is very similar to revA, the x,y,and z values 
%for R are just adjusted to be inclusive of one entire face of the cube
%then the program iterates the same operation for each face of the cube,
%then adds the flux through each surface for a final result.
%I wouldn't go past more than about 150 partitions, even at 100
%paritions the error is typically less than .01%

coulombConst=8.9875517873681764*10^9;%approximate value of k
epsilon0=8.85418781762039*10^-12;%approximate value of epsilon zero
qEnclosed=1*10^-6;%the value of the point charge enclosed, in coulombs
L=50;%length of one side of cube in meters

numberOfPartitions=100;%number of partitions to split each face of the cube into, in each 
                        %direction. In this case, setting this value at 10
                        %would result in 100 small areas per face of the
                        %cube, so 600 in total.
 
                        
chargeLocation=[ (4/10)*L (-3/10)*L (2.5/10)*L];%note, the furthest you can place the
                                                %charge in any direction is
                                                %L/2, since the cube is
                                                %centered at the origin.
                                                %Else the charge is outside

                                                
if(abs(chargeLocation(1,1))<L/2 && abs(chargeLocation(1,2))<L/2 && abs(chargeLocation(1,3))<L/2)                                           %the cube
            
%center of each partitioned area along y

yValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%center of each partitioned area along z
zValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%generate a similar array for x values
%xValues=-L/2*ones(numberOfPartitions,1);

i=1;
j=1;
k=1;
R=[-L/2-chargeLocation(1,1) yValues(1,1)-chargeLocation(1,2) zValues(1,1)-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;-L/2-chargeLocation(1,1) yValues(1,j)-chargeLocation(1,2) zValues(1,k+1)-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.





a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end


%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 


%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[-magArea 0 0];


%we now know E for each R, now lets sum up the dots
c=1;
sum1=0;
[d, not_needed]=size(E);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum1=sum1+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum1

%now on to face 2

%center of each partitioned area along y

yValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%center of each partitioned area along z
zValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%generate a similar array for x values
xValues=L/2*ones(numberOfPartitions,1);

i=1;
j=1;
k=1;
R=[L/2-chargeLocation(1,1) yValues(1,1)-chargeLocation(1,2) zValues(1,1)-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;L/2-chargeLocation(1,1) yValues(1,j)-chargeLocation(1,2) zValues(1,k+1)-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.

a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end

%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 

%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[magArea 0 0];


%we now know E for each R, now lets sum up the dots
c=1;
sum2=0;
d=size(E,1);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum2=sum2+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum2


%now on to face 3

%center of each partitioned area along y

yValues=L/2*ones(numberOfPartitions,1);

%center of each partitioned area along z
zValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%center of each partitioned area along x
xValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

i=1;
j=1;
k=1;
R=[xValues(1,1)-chargeLocation(1,1) L/2-chargeLocation(1,2) zValues(1,1)-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;xValues(1,j)-chargeLocation(1,1) L/2-chargeLocation(1,2) zValues(1,k+1)-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.


a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end


%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 

%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[0 magArea 0];


%we now know E for each R, now lets sum up the dots
c=1;
sum3=0;
d=size(E,1);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum3=sum3+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum3

%now on to face 4

%center of each partitioned area along y

yValues=-L/2*ones(numberOfPartitions,1);

%center of each partitioned area along z
zValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%generate a similar array for x values
xValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

i=1;
j=1;
k=1;
R=[xValues(1,1)-chargeLocation(1,1) -L/2-chargeLocation(1,2) zValues(1,1)-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;xValues(1,j)-chargeLocation(1,1) -L/2-chargeLocation(1,2) zValues(1,k+1)-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.

a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end

%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 

%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[0 -magArea 0];


%we now know E for each R, now lets sum up the dots
c=1;
sum4=0;
d=size(E,1);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum4=sum4+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum4

%now on to face 5

%center of each partitioned area along y

yValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%center of each partioned area along z
zValues=-L/2*ones(numberOfPartitions,1);

%generate a similar array for x values
xValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

i=1;
j=1;
k=1;
R=[xValues(1,1)-chargeLocation(1,1) yValues(1,1)-chargeLocation(1,2) -L/2-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;xValues(1,j)-chargeLocation(1,1) yValues(1,k+1)-chargeLocation(1,2) -L/2-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.

a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end

%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 

%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[0 0 -magArea];


%we now know E for each R, now lets sum up the dots
c=1;
sum5=0;
d=size(E,1);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum5=sum5+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum5



%now on to face 6

%center of each partitioned area along y

yValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

%center of each partitioned area along z
zValues=L/2*ones(numberOfPartitions,1);

%generate a similar array for x values
xValues=linspace(-L/2+(L/(2*numberOfPartitions)),L/2-(L/(2*numberOfPartitions)),numberOfPartitions);

i=1;
j=1;
k=1;
R=[xValues(1,1)-chargeLocation(1,1) yValues(1,1)-chargeLocation(1,2) L/2-chargeLocation(1,3)];
while(i<numberOfPartitions^2)
if(mod(i/numberOfPartitions,1)==0)
    j=j+1;
    k=0;
  
end 

    R=[R;xValues(1,j)-chargeLocation(1,1) yValues(1,k+1)-chargeLocation(1,2) L/2-chargeLocation(1,3)];
    k=k+1;
    i=i+1;
   


end 
%R is now a column vector, with each row being a vector pointing to the
%center of one of the partitions, or da.

a=2;
rMagStore=sqrt(R(1,1)^2+R(1,2)^2+R(1,3)^2);
while(a<=numberOfPartitions^2)
rMag=sqrt(R(a,1)^2+R(a,2)^2+R(a,3)^2);
rMagStore=[rMagStore;rMag];
a=a+1;
end

%rMagStore contains the magnitude of each R vector 

%now we need E for each R 

b=2;
E=(coulombConst*qEnclosed/rMagStore(1)^3)*R(1,:);

while(b<=numberOfPartitions^2)
 E=[E;(coulombConst*qEnclosed/rMagStore(b)^3)*R(b,:)];   
    
  b=b+1;  
    
    
end 

%now we need to generate the magnitude of each small area


magArea=(L/numberOfPartitions)^2;


areaVector=[0 0 magArea];


%we now know E for each R, now lets sum up the dots
c=1;
sum6=0;
d=size(E,1);

while(c<=d)
    
 dotter=dot(E(c,:),areaVector);   
 sum6=sum6+dotter; 
    
 c=c+1;
end 

%this has now summed up the dot products for one face of the cube,
%it is stored in the variable sum6


approximation=sum1+sum2+sum3+sum4+sum5+sum6;



actual=qEnclosed/epsilon0;

error=(abs(actual-approximation)/actual)*100;

 

disp("The approximated flux is: "+approximation+" V*m");

disp("The actual flux is: "+actual+" V*m");

disp("The percent error is: "+ error+"%");

fprintf("If you would like to reduce this error, increase the number of partitions,however,\nI've found 150 is the max it will handle in a reasonable amount of time.\nAlso, ensure the location of the point charge isn't \nin VERY close proximity to the edge of the cube.\n")


else 
    disp("nonono... the charge must be inside the cube(i.e. ensure you dont exceed L/2), if you don't the flux is simply 0");

end