%% INPUT SEGMENT
close all
k=2;         %Halbach coefficient (1=dipole, 2= quadrupole,...)
ri=40/1000;    %inner radius [m]
ro=60/1000;   %outer radius [m]
N=16;        %number of segments
type=4;      %type of Mandhala (1=cylindrical segment, 2 = circular, 3= trigonal,..., m= number of vertices...)
d=3/1000;    %distance between magnets [m] is only calculated if a=0 (next)
a=14/1000;    %magnet size [m] if specified d (above) is ignored
phi=pi/type; %rotation of each Mandhala element
L=3/1000;     %axial length [m] 
z0=(ro+ri)/2/sqrt(2*k+4);
Z0=[-z0,z0]; %position of multiple identical rings in z-dimension [m]
BR=1.4;      %remanence [T]
mr=1.05;     %relative permeability 

nx=100;      %number of points for field plots
ny=100;
scale=5;     %scaling for vector plot, pts = nx/scale 

%% CALCULATION 
R=(ro+ri)/2;

%calculation of ideal Halbach
if k==1 
    B0=BR*log(ro/ri);
    unit=' T';
else
    B0=BR*k/(k-1)*(ri^(1-k)-ro^(1-k));
    if k==2
        unit=' T/m';
    else
       unit=[' T/m^',num2str(k-1)];
    endif
endif
%calculation of permittivity effect
theta=(mr+1)/(mr-1);
fm=(1-theta)*theta*ro^2/(ri^2-theta^2*ro^2);
%calculation of segmentation
fs=sin((k+1)*pi/N)/((k+1)*pi/N);

%calculation of polygonal effect
if type ==1
    A=pi*(ro^2-ri^2)-d*N*(ro-ri);
    fM=A/pi/(ro^2-ri^2);
elseif type==2
    if a==0
       a=(ro-ri)/2-d/2;
     else
       d=ro-ri-2*a;
     endif
     A=a^2*pi;
     fM=N*A/pi/(ro^2-ri^2);
else
     if a==0
         a=(ro-ri-d/2)*sin(pi/type);
     elseif
         d=2*(ro-ri-a/sin(pi/type));
     endif
     A=type*a^2/4*cot(pi/type);
     fM=N*A/pi/(ro^2-ri^2);
endif

%calc truncation in 3rd dimension
denomi=(4*R^2+L^2);
switch k
    case 1
        fL=L*(6*R^2+L^2)*denomi^(-3/2);
    case 2
        fL=L*(L^4+10*L^2*R^2+30*R^4)*denomi^(-5/2);
    case 3
        fL=L*(L^6+14*L^4*R^2+70*L^2*R^4+140*R^6)*denomi^(-7/2);
    case 4
        fL=L*(L^8+18*L^6*R^2+126*L^4*R^4+420*L^2*R^6+630*R^8)*denomi^(-9/2);
    case 5
        fL=L*(L^10+22*L^8*R^2+198*L^6*R^4+924*L^4*R^6+2310*L^2*R^8+2772*R^10)*denomi^(-11/2);
    case 6
        fL=L*(L^12+26*L^10*R^2+286*L^8*R^4+1716*L^6*R^6+6006*L^4*R^8+12012*L^2*R^10+12012*R^12)*denomi^(-13/2)
    otherwise
        fL=0;
endswitch
Xx=linspace(-ri,ri,nx);
Yy=linspace(-ri,ri,ny);
Bx=zeros(nx,ny);
By=Bx;
for ix=1:nx
    for iy=1:ny
        if (Xx(ix)^2+Yy(iy)^2)<=ri^2
            rk=complex(Xx(ix),Yy(iy));
            if k==1
                Bx(ix,iy)=B0*fs*fm*fM*fL;
            else
                Bx(ix,iy)=B0*fs*fm*fM*fL*real(rk^(k-1));
                By(ix,iy)=-B0*fs*fm*fM*fL*imag(rk^(k-1));
            endif
        endif
    endfor
endfor

%% OUTPUT VALUES
disp(['______________________________________________________________________________________________________________'])
pole={'DIPOLE','QUADRUPOLE','HEXUPOLE','OCTUPOLE','DECAPOLE','DODECAPOLE'};
if k<=6
   s{1}=['cylindrical inner Halbach  ',pole{k},' :    with k = ',num2str(k),' (',num2str(2*k),' poles)'];
else
   s{1}=['cylindrical inner Halbach  ',num2str(2*k),'-UPOLE :    with k = ',num2str(k),' (',num2str(2*k),' poles)'];
endif
s{2}=['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '];
s{3}=['inner radius = ',num2str(ri*1000),' mm,  outer radius = ',num2str(ro*1000),' mm, central radius = ',num2str(R*1000),' mm'];
s{4}=['distance between magnets = ',num2str(d*1000),' mm, magnet side length = ',num2str(a*1000),' mm;  remanence = ',num2str(BR),' T'];
s{5}=['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '];
flux='B';
for jj=1:k-1
    flux=[flux,char(39)];
endfor
s{6}=['ideal  ',flux,' = ',num2str(B0,'%.3f'),unit];
s{7}=['discretized in ',num2str(N),' segments  ',flux,' = ',num2str(B0*fs,'%.3f'),unit];
s{8}=['using a rel. permeability of ',num2str(mr),'  ',flux,' = ',num2str(B0*fs*fm,'%.3f'),unit];
shape={'cylindrical segments','circles','triangles','squares','pentagons','hexagons','heptagons','octagons','nonagons','decagons'};
if type<=10
    s{9}=['using ',num2str(N),' ',shape{type},' with side length ',num2str(a*1000),' mm  ',flux,' = ',num2str(B0*fs*fm*fM,'%.3f'),unit];
else
    s{9}=['using ',num2str(N),' ',num2str(type),'-ogons  with side length ',num2str(a*1000),' mm  ',flux,' = ',num2str(B0*fs*fm*fM,'%.3f'),unit];
endif
if fL==0
     s{10}=['truncated to a length of ',num2str(L*1000),' mm  ',flux,'  **** for this k no factor calculated '];
else
     s{10}=['truncated to a length of ',num2str(L*1000),' mm  ',flux,' = ',num2str(B0*fs*fm*fM*fL,'%.3f'),unit];
endif
s{11}=['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '];
s{12}='the vertices of the polygons are saved in the variables Px and Py';
s{13}=['2nd plot shows ',num2str(length(Z0)),' rings at positions in variable Z0'];

for jj=1:13
    sl=length(s{jj});
    str='';
    for kk=1:70-sl
      str=[str,' '];
    endfor
    disp([str,s{jj}]);
    
endfor
disp(['______________________________________________________________________________________________________________'])


%% PLOTS 
% make first plot of magnets in plane
subplot(2,3,1)
pa=linspace(0,2*pi,1000);
plot(ro*cos(pa),ro*sin(pa),'c');
hold on;
plot(ri*cos(pa),ri*sin(pa),'c');
plot(R*cos(pa),R*sin(pa),'g--');
Px=zeros(type,N);   %this is where the vertices are stored
Py=Px;

for jj=1:N
    ux=R*cos(jj*2*pi/N);
    uy=R*sin(jj*2*pi/N);
    La=R/10;
    plot(ux,uy,'ro')
    alph=(k+1)*jj*2*pi/N;
    line([ux-La*cos(alph),ux+La*cos(alph)],[uy-La*sin(alph),uy+La*sin(alph)],'color','r')
    plot(ux+La*cos(alph),uy+La*sin(alph),'r.');
    if type==1
         pia=2*pi*jj/N-pi/N+atan(d/2/ri);
         pie=2*pi*jj/N+pi/N-atan(d/2/ri);
         pci=linspace(pia,pie,50);
         plot(ri*cos(pci),ri*sin(pci),'k')
         poa=2*pi*jj/N-pi/N+atan(d/2/ro);
         poe=2*pi*jj/N+pi/N-atan(d/2/ro);
         pco=linspace(poa,poe,50);
         plot(ro*cos(pco),ro*sin(pco),'k')
         line([ri*cos(pci(1)),ro*cos(pco(1))],[ri*sin(pci(1)),ro*sin(pco(1))],'color','k')
         line([ri*cos(pci(end)),ro*cos(pco(end))],[ri*sin(pci(end)),ro*sin(pco(end))],'color','k')
        %plot segments
    elseif type==2
        pc=linspace(0,2*pi,100);
        plot(ux+a*cos(pc),uy+a*sin(pc))
    else
      ru=a/2*csc(pi/type);
      px=[ux-ru*cos(phi-alph)];
      py=[uy+ru*sin(phi-alph)];
      for kk=1:type
         Px(kk,jj)=ux-ru*cos(2*kk*pi/type)*cos(phi-alph)+ru*sin(2*kk*pi/type)*sin(phi-alph);
         Py(kk,jj)=uy+ru*cos(2*kk*pi/type)*sin(phi-alph)+ru*sin(2*kk*pi/type)*cos(phi-alph);
         px=[px,Px(kk,jj)];
         py=[py,Py(kk,jj)];
      endfor
      line(px,py,'Color','k')
    endif
endfor

title('Arrangement of magnets in xy plane')
xlabel('x [m]')
ylabel('y [m]')
axis square
axis tight
hold off;

%second plot of rings in third dimension
z=linspace(min(Z0)-L*2, max(Z0)+L*2,2000);
Bz=zeros(1,length(z));
subplot(2,3,2); 
for kk=1:length(Z0)
    B=R^(2*k+3)*(R^2+(z-Z0(kk)).^2).^(-k-1.5)*B0*fs*fm*fM*fL;
    plot(z,B,'c');hold on
    Bz=Bz+B;
endfor
plot(z,Bz,'b','Linewidth',2)
title([flux,' along cylinder axis, z'])
ylabel([flux,' [',unit,']'])
xlabel('z [m]')
axis tight
A=B0*fs*fm*fM*fL/2;
for kk=1:length(Z0)
    rectangle('Position',[Z0(kk)-L/2,0,L,A],'Linestyle',':');
endfor
hold off

subplot(2,3,3)
colormap jet;
Ba=abs(sqrt(Bx.^2+By.^2))';
imagesc(Xx,Yy,Ba)
title('magnitude of flux')
xlabel('x [m]')
ylabel('y [m]')
set(gca,'YDir','normal')
axis equal
axis tight

subplot(2,3,4)
imagesc(Xx,Yy,Bx')
title('B_x component')
xlabel('x [m]')
ylabel('y [m]')
set(gca,'YDir','normal')
axis equal
axis tight

subplot(2,3,5)
imagesc(Xx,Yy,By')
title('B_y component')
xlabel('x [m]')
ylabel('y [m]')
set(gca,'YDir','normal')
axis equal
axis tight

subplot(2,3,6)
%hold on
imagesc(Xx,Yy,Ba);
hold on
dx=floor(nx/scale);
dy=floor(ny/scale);
Qx=zeros(dx,dy);
Qy=Qx;
Xq=zeros(dx,1);
Yq=zeros(dy,1);

jx=0;
for ix=1:scale:nx
    jx=jx+1;
    Xq(jx)=Xx(ix);
    jy=0;
    for iy=1:scale:ny
        jy=jy+1;
        Qx(jx,jy)=Bx(ix,iy);
        Qy(jx,jy)=By(ix,iy);
        Yq(jy)=Yy(iy);
    
    endfor
endfor
quiver(Xq,Yq,Qx',Qy','Color','k')
xlabel('x [m]')
ylabel('y [m]')
set(gca,'YDir','normal')
title('magnitude with vectors')
axis equal
axis tight
colorbar
hold off
%%