(*MATHEMATICA companion to Lab 201: Field from point charges*)

(*0. Commands and graphing objects used:
  Table
Graphics
Arrow - to plot vectors
Point
Disk
Line
Show - to combine plots and pictures
Export
Norm - length of a vector
VectorPlot (plots directional field)
StreamPlot - plots field lines
FullSimplify[...,Assumptions->
x ∈ Reals] - to simplify ... assuming x is real (not complex) 
%/.{a->..,b-> ..}- 
replacement command; mostly used when symbolic parameters a, b, ... in the previous expression (%) are to be replaced by speciific values
[[1]] - the 1st element of a list
/@ - advanced command, which indicates that the operation is applied to every element of the list. Useful for graphing multiple objects. E.g. if 'list' is a long list of {x,y} coordinates, the command Graphics[Point /@ list] will show an array of points with indicated coordinates. NOT CRUCIAL at this stage.(Also, see files IntroToMathematicaI and II on our web page https://web.njit.edu/~vitaly/121/ *)

(*Units: later we will learn to carry out transformations of equations and computations for variables with units. At the moment, we will use only SI base units*)

repSI:={ke->9. 10^9, nC->10.^-9}

(*I. Field from a single charge: (OverVector[E]) = (Subscript[k, e])q(OverVector[r])/(r^3)*)

efield[rvec_,q_]:=ke q rvec/Norm[rvec]^3

rvec:={x,y}

efield[rvec,1nC]/.repSI

ef=FullSimplify[%, Assumptions->{x ∈ Reals, y ∈ Reals}]

fig1P=VectorPlot[ef,{x,-1,1},{y,-1,1}]

(*and for negative charge of -1 nC:*)

fig1N=VectorPlot[-ef,{x,-1,1},{y,-1,1}]

(*Note: Mathematica gives all vectors the same length and indicates their magnitude by color. Why?? It would very hard to assign vectors proportional lenths due to strong singularity at r->0. Let us try by "hand" as described below (this is an advanced part which can be followed superficially)*)

(*first, I define a list of points from which to start vectors. To comply with symmetry, I arrange the points in 4 circles (12 points in each)with radii 1,2,3,4). Each circle is rotated by 1/4 of a notch to minimize overlap of future vectors. I use the real and imaginary parts of a complex number to locate a point, but basically this is the same as cos and sin in polar coordinates *)

list1=Table[{Re[Exp[2Pi
I n/12]],Im[Exp[2Pi I n/12]]},{n,0,11}];list2=Table[{2Re[Exp[2Pi
I n/12]],2Im[Exp[2Pi I n/12]]},{n,-3/4,11}];list3=Table[{3Re[Exp[2Pi
I n/12]],3Im[Exp[2Pi I n/12]]},{n,-1/2,11}];list4=Table[{4Re[Exp[2Pi
I n/12]],4Im[Exp[2Pi I n/12]]},{n,-1/4,11}];

(*now let's have a look:*)

Graphics[Point/@{list1,list2,list3,list4}]

(*now I assign vectors to each point; scaling of lengths is arbitrary and is given by parameter a. Note that Arrow ("vector") depends on 2 pairs of coordinates - beginning and end- hense the Transpose command*)

a=.5;arlist1=Transpose[{list1,(1+a*16)list1}];arlist2=Transpose[{list2,(1+a*1/2*4)list2}];arlist3=Transpose[{list3,(1+a*1/3*16/9)list3}];arlist4=Transpose[{list4,(1+1/4*a)list4}];

(*now I add colors to distinguish the origin of vectors:*)

Graphics[{Red,Point/@list1,Arrow/@arlist1,Orange,Point/@list2,Arrow/@arlist2,Green,Point/@list3,Arrow/@arlist3,Blue,Point/@list4,Arrow/@arlist4}]

(*Note the difficulty of plotting scaled vectors due to singularity at r=0. Things get even vorse for a negative charge or for a system of charges. 
Thus NEED field lines*)

StreamPlot[ef,{x,-2,2},{y,-2,2}]

fig1Pstream=Show[%,Graphics[{Red,Disk[{0,0},.1]}]]

(*much better; now the same for negative charge*)StreamPlot[-ef,{x,-2,2},{y,-2,2}];

fig1Nstream=Show[%,Graphics[{Blue,Disk[{0,0},.1]}]]

(*Two charges: field on x-axis, as in Fig.4 of the Manual *)

Clear[a,x,y,r1,r2];a=1;y=0; r1={-a,0}; r2={a,0}; q1=q2=1nC;

efield[rvec-r1,q1]+efield[rvec-r2,q2]/.y->0 
(* which is (Subscript[k, e])q*(((x-1)/(|x-1((|)^3)))+((x+1)/(|x+1((|)^3))))

ef2P=%[[1]]/.repSI

Plot[ef2P,{x,-2,2},PlotRange->{-1000,1000}]

Show[%, Graphics[{Dashed, Line[{{-1,-1000},{-1,1000}}],Line[{{1,-1000},{1,1000}}]}],GridLines->Automatic,PlotLabel->"Total electrric field vs x "]

figef2P=%;

(*Note enormous increase near charges!!*)

(*answer QUESTIONS 1- 4 on pp.23,24 and replot the graph for a dipole with q2=-1 nC*)

(*Part II. Field in 2D - see fig.5 in the Manual*)

Clear[a,x,y,r1,r2];a=1; r1={0,a}; r2={0,-a}; q1=q2=1nC;

efield[rvec-r1,q1]+efield[rvec-r2,q2]

ef2D=FullSimplify[%, Assumptions->{x ∈ Reals, y ∈ Reals}]/.repSI

(*Calculation of field at arb. point*)
ef2D/.{x->1,y->.45}(*as in fig.5*)

VectorPlot[ef2D,{x,-2,2},{y,-2,2}]

figVecPlot2D==%;

Export["figVecPlot2D.jpeg",%,"JPEG"]

(*TRY DIFFERENT FORMATS:   select the best. Print out  separate figure (or page from notebook). Sketch lines by hand.Compare with below*)

StreamPlot[ef2D,{x,-2,2},{y,-2,2}]

(*COMPLETE  THE TASKS 4-6 on p.26 for a dipole. Note that VectorPlot in Mathematica is the closest (not identical) to quiver  in MATLAB*)