aa = 260/2;        % semi major axis (along z)
bb = 170/2;        % semi minor axis of 2 sides
R = aa;            % Asteroid Radius
n = 101;
nr=2;
% [X,Y,Z] = sphere(n);
th = linspace(0,pi,n);   % elevation node positions
phi = linspace(0,2*pi,2*n); % azimuthal node positions
dr=R*(phi(2)-phi(1));
rs = linspace(R-dr,R,nr); % radial node positions

[rr,tt,pp] = ndgrid(rs,th,phi); % Grid of node locations
xx=rr.*sin(tt).*cos(pp);        % X locations of nodes
yy=rr.*sin(tt).*sin(pp);        % Y locations of nodes
zz=rr.*cos(tt);                 % Z locations of nodes
%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)); % volume weighted zone 'radius'
phic=0.5*(phi(1:end-1)+phi(2:end)); % zone azimuth
thc=acos(0.5*(cos(th(1:end-1))+cos(th(2:end)))); %zone elevation

% zzc=0.125*((rs(2:end).^4).*(cos(th(2:end)).^2).*phi(2:end)-...
%     rs(1:end-1).^4.*cos(th(1:end-1)).^2.*phi(1:end-1))./((-1.0/3.0)*(...
%     rs(2:end).^3.*cos(th(2:end)).*phi(2:end)-rs(1: colorbar;end-1).^3).* ...
%     cos(th(1:end-1)).*phi(1:end-1));  colorbar;% volume weighted zone 'x'
% yyc=0.0625*(rs(2:end).^4.*(2.*th(2:end)-sin(2.*th(2:end))).*sin(phi(2:end))-...c
% 5

%     rs(1:end-1).^4.*(2.*th(1:end-1)-sin(2.*th(1:end-1))).*sin(phi(1:end-1)))./((1.0/3.0)*(...
%     rs(2:end).^3.*cos(th(2:end)).*phi(2:end)-rs(1:end-1).^3).* ...
%     cos(th(1:end-1)).*phi(1:end-1)); % zone y
% xxc=0.0625*(rs(2:end).^4.*(2.*th(2:end)-sin(2.*th(2:end))).*cos(phi(2:end))-...
%     rs(1:end-1).^4.*(2.*th(1:end-1)-sin(2.*th(1:end-1))).*cos(phi(1:end-1)))./((1.0/3.0)*(...
%     rs(2:end).^3.*cos(th(2:end)).*phi(2:end)-rs(1:end-1).^3).* ...
%     cos(th(1:end-1)).*phi(1:end-1)); % zone z

volc=-1/3*reshape((rs(2:end).^3-rs(1:end-1).^3),[numel(rs)-1,1,1]).*reshape(cos(th(2:end)) ...
    -cos(th(1:end-1)),[1,numel(th)-1,1]).*reshape(phi(2:end)-phi(1:end-1), [1,1,numel(phi)-1]) ;

zzc = (-0.125*reshape((rs(2:end).^4-rs(1:end-1).^4),[numel(rs)-1,1,1]).*reshape(cos(th(2:end)).^2 ...
    -cos(th(1:end-1)).^2,[1,numel(th)-1,1]).*reshape(phi(2:end)-phi(1:end-1), [1,1,numel(phi)-1]))./volc;
xxc = (0.0625*reshape((rs(2:end).^4-rs(1:end-1).^4),[numel(rs)-1,1,1]).*reshape(2*th(2:end)-sin(2*th(2:end)) ...
    -(2*th(1:end-1)-sin(2*th(1:end-1))),[1,numel(th)-1,1]).*reshape(sin(phi(2:end))-sin(phi(1:end-1)), [1,1,numel(phi)-1]))./volc;
yyc = (-0.0625*reshape((rs(2:end).^4-rs(1:end-1).^4),[numel(rs)-1,1,1]).*reshape(2*th(2:end)-sin(2*th(2:end)) ...
    -(2*th(1:end-1)-sin(2*th(1:end-1))),[1,numel(th)-1,1]).*reshape(cos(phi(2:end))-cos(phi(1:end-1)), [1,1,numel(phi)-1]))./volc;

% volc=reshape((rs(2:end)-rs(1:end-1)),[numel(rs)-1,1,1]).*reshape((0.5*(rs(2:end)+rs(1:end-1))).^2,[numel(rs)-1,1,1]).*...
%     reshape(sin(0.5*(th(2:end)+th(1:end-1))),[1,numel(th)-1,1]).*reshape((th(2:end)-th(1:end-1)),[1,numel(th)-1,1]) ...
%     .*reshape(phi(2:end)-phi(1:end-1), [1,1,numel(phi)-1]) ;

[rrc,ttc,ppc] = ndgrid(rsc,thc,phic); % Mesh of zone centers

%% Uncomment this for random mesh
%rrc=rrc*0+R;
%ttc=acos(rand(size(ttc))*2-1);
%ppc=2*pi*rand(size(ppc));
%volc=volc.*0+1;
%%
%xxc=rrc.*sin(ttc).*cos(ppc);          % X locations of zonal mesh
%yyc=rrc.*sin(ttc).*sin(ppc);          % Y
%zzc=rrc.*cos(ttc);   

pos_rd = cat(2,xxc(:),yyc(:),zzc(:))';  % Creates a 3xN array of zone center positions
% m_rd = pos_rd*0;               % Magnetic field at dipole

%% Calculate dipole field at zone centers
% AM = [0;0;1e14];            % Dipole moment of the core
AM = [1e14;0;0];            % Dipole moment of the core
r=pos_rd-[0;0;0];
rm = sum(r.^2,1).^(1/2);
m_rd=3*r.*sum(AM.*r,1)./rm.^5-AM./rm.^3;

%%
%sum(r(1,:))
%sum(r(2,:))
%sum(r(3,:))
%sum(r.*r(1,:)./rm.^5.*volc(1:end),2)
%sum(r.*r(2,:)./rm.^5.*volc(1:end),2)
%sum(r.*r(3,:)./rm.^5.*volc(1:end),2)

%% If we want to do a pure dipole
%pos_rd=pos_rd*0;
%m_rd=m_rd*0+AM;

%% Plot the dipole locations
figure,
% quiver3(pos_rd(1,1:1:end),pos_rd(2,1:1:end),pos_rd(3,1:1:end),m_rd(1,1:1:end), ...
%     m_rd(2,1:1:end),m_rd(3,1:1:end));
scatter3(pos_rd(1,1:1:end),pos_rd(2,1:1:end),pos_rd(3,1:1:end),3,volc(1:1:end));
axis equal;
xlabel('X'); ylabel('Y'); zlabel('Z');
%% Define the grid
npoints=101;
xx1=linspace(-5*R,5*R,npoints);
yy1=xx1;
zz1=xx1;
sz=length(xx1);
[X,Y,Z]=ndgrid(xx1,yy1,zz1);
b=zeros(sz,sz,sz,3);  
mask=(X.^2+Y.^2+Z.^2)<(R+4*dr)^2;
r2=0;
pos_mesh = cat(2,X(:),Y(:),Z(:))'; 
%% So we have two 3D meshes - and we want to calculate the position vector between each pair
%dx=xx1(:)-pos_rd(1,:);
%dy=yy1(:)-pos_rd(2,:);
%dz=zz1(:)-pos_rd(3,:);

%% Calculate the multiple dipole field
bd1d=zeros(size(pos_mesh));
for i=1:size(pos_rd,2) % Loop over dipoles
%     pos_rd(:,i)=[250;0;0];
%     m_rd(:,i)=[0;1e14;0];
    r=pos_mesh-pos_rd(:,i);   
    rm=sum(r.^2,1).^(1/2);
    bd1d=bd1d+(3*r.*sum(m_rd(:,i).*r,1)./rm.^5-m_rd(:,i)./rm.^3).*volc(i);   
    %bd1d=bd1d+(3*r.*sum(m_rd(:,i).*r,1)./rm.^3-m_rd(:,i)./rm).*volc(i);
    %bd1d=r;    
end
% bd1d(:,mask)=0;
%% Debugging Center calculation
% figure,
% subplot(1,2,1),
% plot3(xxc(1,5,4),yyc(1,5,4),zzc(1,5,4),'b+');
% scatter3(reshape(xxc(1:end,1:end,1:end),1,[]),reshape(yyc(1:end,1:end,1:end),1,[]),reshape(zzc(1:end,1:end,1:end),1,[]));
% hold on;
% scatter3(reshape(xx(1:2,5:6,4:5),1,[]),reshape(yy(1:2,5:6,4:5),1,[]),reshape(zz(1:2,5:6,4:5),1,[]));

%%
writebov3dvector('BxField_51_y125',min(pos_mesh,[],2), max(pos_mesh,[],2), 0, reshape(bd1d.',[sz,sz,sz,3]), 'B');

function writebov3dvector(name, lower, upper, tnow, data, varname)
    f=fopen([name,'.bov'],'w');
    datafile=[name,'.dat'];
    fprintf(f, 'TIME: %f\n', tnow);
    fprintf(f, 'DATA_FILE: %s\n', datafile);
    fprintf(f, 'DATA_SIZE: %d %d %d\n', size(data,1), size(data, 2), size(data,3));
    fprintf(f, 'DATA_FORMAT: DOUBLE\n');
    fprintf(f, 'VARIABLE: %s\n',varname);
    fprintf(f, 'DATA_ENDIAN: LITTLE\n');
    fprintf(f, 'CENTERING: zonal\n');
    fprintf(f, 'BRICK_ORIGIN: %e %e %e\n',lower);
    fprintf(f, 'BRICK_SIZE: %e %E %e\n', upper-lower);
    fprintf(f, 'BYTE_OFFSET: 0\n');
    fprintf(f, 'DATA_COMPONENTS: %d', size(data, 4));
    fclose(f);
    f=fopen(datafile, 'wb');
    fwrite(f, permute(data,[4,1,2,3]),'double');
    fclose(f);
end

% %% Old stuff...  follows
% bd=reshape(bd1d,[3,sz,sz,sz]);
% %%
% bdd=b*0;
% for ii=1:sz
%     for jj=1:sz
%         for kk=1:sz
%             
% %            if sqrt(sum(X(ii,jj,kk)^2+Y(ii,jj,kk)^2+Z(ii,jj,kk)^2)) < R*1.01
% %                continue
% %            end
%             for ll=1:length(m_rd)
%                 r = [X(ii,jj,kk),Y(ii,jj,kk),Z(ii,jj,kk)]-pos_rd(:,ll)';
%                 r2 = sqrt(sum(r.^2));
%                 rm = r2;
%                 %if (r2 < 1e-6)
%                 %    continue
%                 %end
%                 tmp=((3*r*dot(m_rd(:,ll)',r)/rm^5)-m_rd(:,ll)'/rm^3);
%                 %tmp=((3*r*dot([0,1,0],r)/rm^5)-[0,1,0]/rm^3);    % Test
%                 % with uniform field
%                 bdd(ii,jj,kk,:) = bdd(ii,jj,kk,:)+reshape(tmp,size(bdd(ii,jj,kk,:)));
%             end
%         end
%     end
% end
% %bd = k*bd;
% figure,
% %pcolor(squeeze(X),squeeze(Z),squeeze(log10(sqrt(sum(bd(:,1,:,:).^2,4))))); colorbar;
% 
% %% Direct Dipole Field Around The Asteroid
% k = 1e-7;
% bdd = b*0;
% for ii=1:sz
%     for jj=1:sz
%         for kk=1:sz
%             if sqrt(sum(X(ii,jj,kk)^2+Y(ii,jj,kk)^2+Z(ii,jj,kk)^2)) < R*1.01
%                 continue
%             end
%             r = [X(ii,jj,kk),Y(ii,jj,kk),Z(ii,jj,kk)]-[0,0,0];
%             r2 = sqrt(sum(r.^2));
% %             if (r2 < 1e-6)
% %                 continue
% %             end
%             tmp=(3*r*dot(AM',r)/r2^5)-AM'/r2^3;
%             bdd(ii,jj,kk,:) = bdd(ii,jj,kk,:)+reshape(tmp,size(bdd(ii,jj,kk,:)));
%         end
%     end
% end
% %%
% figure
% clf
% icomp=3
% iz=1
% subplot(1,2,1)
% imagesc(squeeze(bd(icomp,:,:,iz)))
% colorbar
% subplot(1,2,2)
% imagesc(squeeze(bdd(:,:,iz,icomp)))
% colorbar
% %bd = k*bd;
% %%
% figure,
% pcolor(squeeze(X),squeeze(Z),squeeze(log10(sqrt(sum(bdd(:,1,:,:).^2,4))))); colorbar;