%1) Go ahead and define constants we'll need and the dipole moment
const1 = ((4*pi)^(-1))*sqrt((4*pi*(10)^(-7))*(8.85*(10)^(-12))^(-1));
const2 = ((2*pi)^(-1));
const3 = 0.5*const2;
mu = 4*pi*((10)^(-7));
const4 = 180 * ((pi)^(-1));
epsilon = 8.85*((10)^(-12));
kwave = 0.001*5*pi;
const5 = (4*pi*epsilon)^(-1);

%2) Form spherical meshgrid, avoid redundancies and singularities

n = 25;
rmin = 0.2;
rmax = 2.0;
phi = linspace(-pi,pi - ((n)^(-1))*2*pi,n);
theta = linspace(-0.5*pi,0.5*pi - ((n)^(-1))*2*pi,n);
r = linspace(rmin,rmax,n);
p = 0.001*(r(n-1)-r(n-2));
p = [0;0;p];
%3) Create a meshgrid of sherical coordinates, extract vectors, and transpose


[Phi,Theta,R] = meshgrid(phi,theta,r);
SphericalPointsfromMesh = [Phi(:),Theta(:),R(:)];
SphericalPointsfromMeshT = SphericalPointsfromMesh';

%4) Create a repeated matrix of p

PmatC = repmat(p,[1 size(SphericalPointsfromMeshT,2)]);

%5) We want to convert PmatC to spherical points at each point given in SphericalPointsfromMeshT

%5-1) In order to do that we need to convert SphericalPointsfromMeshT to degrees

SphericalPointsfromMeshD = [const4*Phi(:),const4*Theta(:),R(:)];
SphericalPointsfromMeshDT = SphericalPointsfromMeshD';
azdegrees = SphericalPointsfromMeshDT(1,:);
eldegrees = SphericalPointsfromMeshDT(2,:);


%5-2) I think we're ready to do what we set out to do in (5)

for i1 = 1:size(SphericalPointsfromMeshDT,2)
    PmatS(:,i1) = cart2sphvec(p,azdegrees(1,i1),eldegrees(1,i1));
end 

%6) Now to create the B components in spherical basis

for i2 = 1:size(SphericalPointsfromMeshDT,2)
    Baz(1,i2) = - const1 * exp(1i*kwave*SphericalPointsfromMeshDT(3,i2))*(SphericalPointsfromMeshDT(3,i2)^(-3))*((kwave^2)*((SphericalPointsfromMeshDT(3,i2)^(2)) + 1i*kwave*SphericalPointsfromMeshDT(3,i2))*(PmatS(2,i2)));
end

for i3 = 1:size(SphericalPointsfromMeshDT,2)
    Bel(1,i3) =  const1 * exp(1i*kwave*SphericalPointsfromMeshDT(3,i3))*(SphericalPointsfromMeshDT(3,i3)^(-3))*((kwave^2)*((SphericalPointsfromMeshDT(3,i3)^(2)) + 1i*kwave*SphericalPointsfromMeshDT(3,i3))*(PmatS(1,i3)));
end

Bsr = zeros(size(Bel));

% Vertically concatenate

Bs = vertcat(Baz,Bel,Bsr);

% Convert Bs back to cartesian (at least try to)

for i4 = 1:size(Bs,2)
    BC(:,i4) = sph2cartvec(Bs(:,i4),azdegrees(1,i4),eldegrees(1,i4));
end


% Now to create E 

for i5 = 1:size(SphericalPointsfromMeshDT,2)
    Eaz(1,i5) = const5*exp(1i*kwave*SphericalPointsfromMeshDT(3,i5))*(SphericalPointsfromMeshDT(3,i5)^(-3))*((kwave^2)*((SphericalPointsfromMeshDT(3,i5)^(2)) + 1i*kwave*SphericalPointsfromMeshDT(3,i5) -1))*(PmatS(1,i5));
end

for i6 = size(SphericalPointsfromMeshDT,2)
    Eel(1,i6) = const5*exp(1i*kwave*SphericalPointsfromMeshDT(3,i6))*(SphericalPointsfromMeshDT(3,i6)^(-3))*((kwave^2)*((SphericalPointsfromMeshDT(3,i6)^(2)) + 1i*kwave*SphericalPointsfromMeshDT(3,i6) -1))*(PmatS(2,i6));
end

for i7 = 1:size(SphericalPointsfromMeshDT,2)
    Er(1,i7) = 0.5*const1*exp(1i*kwave*SphericalPointsfromMeshDT(3,i7))*(SphericalPointsfromMeshDT(3,i7)^(-3))*(-1i*kwave*SphericalPointsfromMeshDT(3,i7) + 1)*(PmatS(3,i7));
end

% Vertically concatenate

Es = vertcat(Eaz,Eel,Er);

% Convert Es back to Cartesian (at least try to)

for i8 = 1:size(Es,2)
    EC(:,i8) = sph2cartvec(Es(:,i8),azdegrees(1,i8),eldegrees(1,i8));
end

% Now calculate the Poynting Vector

H = ((mu)^(-1))*BC;
Hconj = conj(H);

EcrossHconj = cross(EC,Hconj);

Realpart = real(EcrossHconj);

FinalPoyntingVector = 0.5*Realpart;

% Now to plot the vector field and hopefully get something that is not ridiculous

%first have to convert the entire position grid to cartesian

[X,Y,Z] = sph2cart(Phi,Theta,R);

CartesianPointsfromMesh = [X(:),Y(:),Z(:)];

CartesianPointsfromMeshT = CartesianPointsfromMesh';

% Now to hopefully plot the vector field and not get something ridiculous

quiver3(CartesianPointsfromMeshT(1,:),CartesianPointsfromMeshT(2,:),CartesianPointsfromMeshT(3,:),FinalPoyntingVector(1,:),FinalPoyntingVector(2,:),FinalPoyntingVector(3,:));