%% 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 %%