%% INPUT SEGMENT k=2; %Halbach coefficient (1=dipole, 2= quadrupole,...) ri=5/100; %inner radius [m] ro=8/100; %outer radius [m] N=16; %number of segments type=4; %type of Mandhala (1=cylindrical segment, 2 = circular, 3= trigonal,..., m= number of vertices...) a=20/1000; %side length of magnet [m], if = 0, a is calculated from d (next) d=7/100; %distance between magnets [m], if a~=0 d is ignored 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.3; %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 phi=pi/type; %rotation of each Mandhala element R=(ro+ri)/2; %calculation of ideal Halbach if k==1 B0=BR*log(ro/ri)*1000; unit=' mT'; 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)]; end end %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 a~=0 if type ==1 a=0; elseif type==2 d=ro-ri-2*a; else d=2*(ro-ri-a/sin(pi/type)); end end if type ==1 A=pi*(ro^2-ri^2)-d*N*(ro-ri); fM=A/pi/(ro^2-ri^2); elseif type==2 a=(ro-ri)/2-d/2; A=a^2*pi; fM=N*A/pi/(ro^2-ri^2); else a=(ro-ri-d/2)*sin(pi/type); A=type*a^2/4*cot(pi/type); fM=N*A/pi/(ro^2-ri^2); end %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; end 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)); end end end end %% 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)']; end 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 size = ',num2str(a*1000),' mm; remanence = ',num2str(BR),' T']; s{5}=['- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ']; flux='B'; for jj=1:k-1 flux=[flux,char(39)]; end 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]; end 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]; end 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']; s=pad(s,60,'left'); for jj=1:13 disp(s{jj}); end 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),'k') 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)]; end line(px,py,'Color','k') end end title('Arrangement of magnets in xy plane') xlabel('x [m]') ylabel('y [m]') axis equal 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; end 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',':'); end hold off subplot(2,3,3) 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); end end 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 %%