theta45=asin(sqrt(2)/sqrt(3))*180/pi;
sim_thetas=[0,45,theta45,45,90,90,90];
sim_phis=  [0, 0,     45,90, 0,45,90];

% Array of Names
% Array which says face=2, edge=4, corner=8 [Cases].
%fnl = strings(1,length(sim_phis));
fnl = "fort."+round(sim_thetas(:))+"_"+round(sim_phis(:))+".txt";

Rx=[-1,0,0;0,1,0;0,0,1]; %Have to do a coordinate transformation in asteroid system z-> -z (opposite spin direction)
% Swap Z coordinate and z component of the field (poloidal field). Toroidal field??
Ry=[-1,0,0;0,1,0;0,0,-1];
Rz=[1,0,0;0,1,0;0,0,-1];
I=eye(3);
sim_transformations = {{I,Rz},{I,Rx,Rz,Rx*Rz},{I,Rx,Ry,Rz,Rx*Ry,Rx*Rz,Ry*Rz,Rx*Ry*Rz},...
    {I,Ry,Rz,Ry*Rz},{I,Rx},{I,Rx,Ry,Rx*Ry},{I,Ry}};
cases = cellfun(@numel, sim_transformations);

%

%N = 25636;       % Number of elements, maybe take as input from file?
A = load('fort.90_45.txt');
pos = A(:,1:3);             % Postition of dipole's center
m = A(:,4:6);               % Magnetic field at dipole
% We are currently excluding the normalization to the dipole moment of each
% dipole element/cell
theta_=90; 
phi_=45; 


R = 125;
n = 101;
a_e = [125; 125*0.9];
b_e = [85; 85*0.9];


% [X,Y,Z] = sphere(n);
th = linspace(0,pi,n);
phi = linspace(0,2*pi,n);
rs = linspace(R-10,R,2);

[rr,tt,pp] = ndgrid(rs,th,phi);
for ii =1:2
    rr(ii,:,:)=1./sqrt((sin(tt(ii,:,:)).*cos(pp(ii,:,:))./a_e(ii)).^2+ ...
        (sin(tt(ii,:,:)).*sin(pp(ii,:,:))./a_e(ii)).^2+(cos(tt(ii,:,:))./b_e(ii)).^2);
end

% rr=rotz(phi_)*roty(theta_)*rr;
% tt=rotz(phi_)*roty(theta_)*tt;
% pp=rotz(phi_)*roty(theta_)*pp;

% xx=rr.*sin(tt).*cos(pp);
% yy=rr.*sin(tt).*sin(pp);
% zz=rr.*cos(tt);
%rsc=0.5*(rs(1:end-1)+rs(2:end));

%

rsc=0.25*(rs(2:end).^4-rs(1:end-1).^4)./((1.0/3.0)*(rs(2:end).^3-rs(1:end-1).^3));
phic=0.5*(phi(1:end-1)+phi(2:end));
thc=acos(0.5*(cos(th(1:end-1))+cos(th(2:end))));

[rrc,ttc,ppc] = ndgrid(rsc,thc,phic);
for ii =1:1
    rrc(ii,:,:)=1./sqrt((sin(ttc(ii,:,:)).*cos(ppc(ii,:,:))./a_e(ii)).^2+ ...
        (sin(ttc(ii,:,:)).*sin(ppc(ii,:,:))./a_e(ii)).^2+(cos(ttc(ii,:,:))./b_e(ii)).^2);
end
xxc=rrc.*sin(ttc).*cos(ppc);
yyc=rrc.*sin(ttc).*sin(ppc);
zzc=rrc.*cos(ttc);

% xxc=rotz(phi_)*roty(theta_)*(rrc.*sin(ttc).*cos(ppc));
% yyc=rotz(phi_)*roty(theta_)*(rrc.*sin(ttc).*sin(ppc));
% zzc=rotz(phi_)*roty(theta_)*(rrc.*cos(ttc));

%

%za = 45; ya = 90;
za = 0; ya = 0; xa = 0;
% pos_r = roty(-ya)*rotz(-za)*pos';
% m_r = roty(-ya)*rotz(-za)*m';
%pos_r = rotx(xa)*rotz(-za)*pos';
tn_=theta_+180;
pn_=phi_; %+180;
pos_r = roty(-tn_)*rotz(-pn_)*pos'; % for Rx switch sign of all z components.
m_r = roty(-tn_)*rotz(-pn_)*m';     % for Rx negate Bz AND (Bx AND By) --> might not matter

F1 = scatteredInterpolant(pos_r',m_r(1,:)','linear','nearest');
F2 = scatteredInterpolant(pos_r',m_r(2,:)','linear','nearest');
F3 = scatteredInterpolant(pos_r',m_r(3,:)','linear','nearest');
mx = reshape(F1(xxc(:),yyc(:),zzc(:)),size(xxc));
my = reshape(F2(xxc(:),yyc(:),zzc(:)),size(xxc));
mz = reshape(F3(xxc(:),yyc(:),zzc(:)),size(xxc));

%

% vc=(-1/3)*(rr(2:end,1:end-1,1:end-1).^3-rr(1:end-1,1:end-1,1:end-1).^3).* ...
%     (cos(tt(1:end-1,2:end,1:end-1))-cos(tt(1:end-1,1:end-1,1:end-1))).* ...
%     (pp(1:end-1,1:end-1,2:end)-pp(1:end-1,1:end-1,1:end-1));    % Volume of cells

figure,
%scatter3(xx(:),yy(:),zz(:),'r.');
scatter3(pos_r(1,1:2:end),pos_r(2,1:2:end),pos_r(3,1:2:end),'r.');
hold on;
scatter3(xxc(:),yyc(:),zzc(:),'b.'); %,vc(:));
hold off;
axis equal;

%

figure,
% quiver3(xxc(:),yyc(:),zzc(:),mx(:),my(:),mz(:));
% hold on;
quiver3(pos_r(1,1:4:end),pos_r(2,1:4:end),pos_r(3,1:4:end),m_r(1,1:4:end),m_r(2,1:4:end),m_r(3,1:4:end));
% quiver3(pos(1:1:end,1),pos(1:1:end,2),pos(1:1:end,3),m(1:1:end,1),m(1:1:end,2),m(1:1:end,3));
% axis equal;
xlabel('X'); ylabel('Y'); zlabel('Z');

%

figure,
% qcp = quiver3(pos(1:4:end,1),pos(1:4:end,2),pos(1:4:end,3),m(1:4:end,1),m(1:4:end,2),m(1:4:end,3));
% qcp = quiver3(pos(1:3:end,1),pos(1:3:end,2),pos(1:3:end,3),m(1:3:end,1),m(1:3:end,2),m(1:3:end,3));
qcp = quiver3(xxc(:),yyc(:),zzc(:),mx(:),my(:),mz(:));

%// Compute the magnitude of the vectors
mags = sqrt(sum(cat(2, qcp.UData(:), qcp.VData(:), ...
            reshape(qcp.WData, numel(qcp.UData), [])).^2, 2));

%// Get the current colormap
currentColormap = colormap(gca);

%// Now determine the color to make each arrow using a colormap
[~, ~, ind] = histcounts(mags, size(currentColormap, 1));

%// Now map this to a colormap to get RGB
cmap = uint8(ind2rgb(ind(:), currentColormap) * 255);
cmap(:,:,4) = 255;
cmap = permute(repmat(cmap, [1 3 1]), [2 1 3]);

%// We repeat each color 3[666.7754 2.0313e+03 2.4064e+03] times (using 1:3 below) 
% because each arrow has 3 vertices
set(qcp.Head, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:3,:,:), [], 4).');   %'

%// We repeat each color 2 times (using 1:2 below) because each tail has 2 vertices
set(qcp.Tail, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:2,:,:), [], 4).');
% Colorbar of the inside sphere is messing with the plot colorbar
% hold on
% [Xsp,Ysp,Zsp] = sphere(n);
% hsp = surf(R*0.95*Xsp,R*0.95*Ysp,R*0.95*Zsp);
% set(hsp,'edgecolor','none','facecolor',[.3 .3 .3]);
% colormap bone; shading interp;
axis equal; colorbar;
xlabel('X'); ylabel('Y'); zlabel('Z');

%

m_sph=zeros(numel(mx),3);
for ii=1:numel(mx)
    m_sph(ii,:) = cart2sphvec([mx(ii);my(ii);mz(ii)],atan2d(yyc(ii),xxc(ii)),...
        asind(zzc(ii)/rrc(ii)));
end

%

%% Averaging (How to handle arbitrary spin? Weights?)
% Azimuthal averaging - average over phi
nt=n-1; np=n-1;
tmp = 0;
m_sphr = reshape(m_sph(:,3),size(xxc));
m_sphth = reshape(m_sph(:,2),size(xxc));
m_sphp = reshape(m_sph(:,1),size(xxc));

mra = repmat(mean(m_sphr,3),[1,1,size(m_sphr,3)]);
% mra = m_sphr, try mean along 2nd dim, [1, size(m_sphr,2), 1]
mtha = repmat(mean(m_sphth,3),[1,1,size(m_sphth,3)]);
mpa = repmat(mean(m_sphp,3),[1,1,size(m_sphp,3)]);

% mra = repmat(mean(m_sphr,2),[1, size(m_sphr,2), 1]);
% mtha = repmat(mean(m_sphth,2),[1, size(m_sphr,2), 1]);
% mpa = repmat(mean(m_sphp,2),[1, size(m_sphr,2), 1]);

%

m_c=zeros(numel(mx),3);
for ii=1:numel(mx)
%     m_c(ii,:) = sph2cartvec([m_sph(ii,1);m_sph(ii,2);m_sph(ii,3)],atan2d(yyc(ii),xxc(ii)),...
%         asind(zzc(ii)/rrc(ii)));
    m_c(ii,:) = sph2cartvec([mpa(ii);mtha(ii);mra(ii)],atan2d(yyc(ii),xxc(ii)),...
        asind(zzc(ii)/rrc(ii)));
end

%

figure,
surf(reshape(xxc,[100 100]),reshape(yyc,[100 100]),...
    reshape(zzc,[100 100]),reshape(mpa,[100 100]));%reshape(m_c(:,1),size(xxc)));
shading interp;
axis equal; colorbar;
xlabel('X'); ylabel('Y'); zlabel('Z');

%

figure,
surf(reshape(xxc,[100 100]),reshape(yyc,[100 100]),...
    reshape(zzc,[100 100]),reshape(m_c(:,1),[100 100]));%reshape(m_c(:,1),size(xxc)));
shading interp;
axis equal; colorbar;

figure,
surf(reshape(xxc,[100 100]),reshape(yyc,[100 100]),...
    reshape(zzc,[100 100]),reshape(m_c(:,2),[100 100]));%reshape(m_c(:,1),size(xxc)));
shading interp;
axis equal; colorbar;

figure,
surf(reshape(xxc,[100 100]),reshape(yyc,[100 100]),...
    reshape(zzc,[100 100]),reshape(m_c(:,3),[100 100]));%reshape(m_c(:,1),size(xxc)));
shading interp;
axis equal; colorbar;

%

figure,
% qcp = quiver3(pos(1:4:end,1),pos(1:4:end,2),pos(1:4:end,3),m(1:4:end,1),m(1:4:end,2),m(1:4:end,3));
% qcp = quiver3(pos(1:3:end,1),pos(1:3:end,2),pos(1:3:end,3),m(1:3:end,1),m(1:3:end,2),m(1:3:end,3));
qcp = quiver3(xxc(:),yyc(:),zzc(:),m_c(:,1),m_c(:,2),m_c(:,3));

%// Compute the magnitude of the vectors
mags = sqrt(sum(cat(2, qcp.UData(:), qcp.VData(:), ...
            reshape(qcp.WData, numel(qcp.UData), [])).^2, 2));

%// Get the current colormap
currentColormap = colormap(gca);

%// Now determine the color to make each arrow using a colormap
[~, ~, ind] = histcounts(mags, size(currentColormap, 1));

%// Now map this to a colormap to get RGB
cmap = uint8(ind2rgb(ind(:), currentColormap) * 255);
cmap(:,:,4) = 255;
cmap = permute(repmat(cmap, [1 3 1]), [2 1 3]);
cb=colorbar;
% cb.Limits=([min(mags) max(mags)]);  
%// We repeat each color 3[666.7754 2.0313e+03 2.4064e+03] times (using 1:3 below) because each arrow has 3 vertices
set(qcp.Head, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:3,:,:), [], 4).');   %'

%// We repeat each color 2 times (using 1:2 below) because each tail has 2 vertices
set(qcp.Tail, 'ColorBinding', 'interpolated', 'ColorData', reshape(cmap(1:2,:,:), [], 4).');
% hold on
% [Xsp,Ysp,Zsp] = sphere(n);
% hsp = surf(R*0.95*Xsp,R*0.95*Ysp,R*0.95*Zsp);
% set(hsp,'edgecolor','none','facecolor',[.3 .3 .3]);
% colormap bone; shading interp;
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z');