Superferric ILC-Optimized CESR-c Wiggler:  12 Poles, Period=32cm, Gap=86mm, Width=238mm, Bpeak=1.95T

Introduction

All sections must be evaluated in the order of presentation. A section may be evaluated several   times with the same or different parameters if the previous sections have been evaluated before. The only limitation is memory. The Radia.exe memory can be re-initialized by re-executing the section entitled "Load and   Initialize Radia". The Kernel memory can be re-initialized by calling Exit[] and then re-executing all sections of this example. If the Front-End is   running out of memory, close some windows or close the Front-End and start it again. Note that if by mistake any section is executed before the previous ones are executed, it may be necessary to exit the Kernel and re-start everything from the very beginning.

Load and Initialize Radia

The following instruction loads the Radia package and ret the Radia  version number.

<<Radia`;

Build the Geometry

These instructions define 1 wiggler, with current density j1 in the coil;

radUtiDelAll[];

RadUtiMem[];

(* Coil Presentation Parameters *)
n1=100;
c1={0,1,1};
thcn=0.001;

(* Coil geometry parameters-in mm *)
Gap=43.2; 
Coilheight=19.0; 
Coilctr=Gap+Coilheight/2; 
Coilradin=29.0; 
Coilthick=25.4; 
Coilradout=Coilthick+Coilradin; 
Coillength=180.;  
Coilwidth=51.2; 
Coilarea=Coilthick*Coilheight; 
Turns=660; 
Endturns=182;
Endcoilwidth=0.1;
Nearendcoilwidth=11.2;
Endcoilthick=25.4*Endturns/Turns;
Endcoilradin=29.0;
Endcoilradout=Endcoilthick+Endcoilradin;
Endcoilheight=19.0;
Endcoilarea=Endcoilthick*Endcoilheight;
Nearendcoilthick=25.4;
Nearendcoilradin=29.0;
Nearendcoilradout=Nearendcoilthick+Nearendcoilradin;
Chamfer=0.0;
Chamferdepth=0.0;
Coillength=Coillength*1.0;

(* Current Densities in A/mm^2 *)
Current=141; 
j1=Current*Turns/Coilarea;
j1end=Current*Endturns/Endcoilarea;
cN=4; 
Clear[nsub]; nsub={3,3,2};
Clear[nsubb]; nsubb={3,3,2};

createapole[length_,width_,ctr_,height_,radin_,radout_,cham_,chamdepth_,gap2_,cN2_,nsub2_,nsubb2_,n12_,c12_,thcn2_]:=

((*Create curved interior pole piece,negative x,negative y*)
flist={{1,2,3},{5,2,3,6},{4,1,3,6},{5,2,1,4},{5,4,6}};
cc1={-length/2,-width/2,ctr-height/2};
cc2={-length/2,-width/2,ctr+height/2};
Clear[polfcnxny];
polfcnxny[n_]:=(plist1={cc1+{-radin*Sin[Pi*n/cN2/2],-radin*Cos[Pi*n/cN2/2],
              0},cc1+{-radin*Sin[Pi*(n+1)/cN2/2],-radin*
                Cos[Pi*(n+1)/cN/2],0},cc1};
      plist2={cc2+{-radin*Sin[Pi*n/cN2/2],-radin*Cos[Pi*n/cN2/2],0},
          cc2+{-radin*Sin[Pi*(n+1)/cN2/2],-radin*Cos[Pi*(n+1)/cN2/2],0},cc2};
      plist=Flatten[{plist1,plist2},1];
      radObjPolyhdr[plist,flist,{0,0,0}]);
Clear[Polenxny];
Polenxny=radObjCnt[Table[polfcnxny[n],{n,0,cN2-1}]];
radObjDrwAtr[Polenxny,{1,0,0},thcn2];

(* Create curved interior pole piece, negative x, positive y *)
flist={{1,2,3},{5,2,3,6},{4,1,3,6},{5,2,1,4},{5,4,6}};
cc1={-length/2,width/2,ctr-height/2};
cc2={-length/2,width/2,ctr+height/2};
Clear[polfcnxpy];
polfcnxpy[n_]:=(plist1={cc1+{-radin*Sin[Pi*n/cN2/2],radin*Cos[Pi*n/cN2/2],0},
            cc1+{-radin*Sin[Pi*(n+1)/cN2/2],radin*Cos[Pi*(n+1)/cN2/2],0},cc1};
          plist2={cc2+{-radin*Sin[Pi*n/cN2/2],radin*Cos[Pi*n/cN2/2],0},
             cc2+{-radin*Sin[Pi*(n+1)/cN/2],radin*Cos[Pi*(n+1)/cN2/2],0},cc2};
          plist=Flatten[{plist1,plist2},1];
          radObjPolyhdr[plist,flist,{0,0,0}]);
Clear[Polenxpy]; 
Polenxpy=radObjCnt[Table[polfcnxpy[n],{n,0,cN2-1}]];
radObjDrwAtr[Polenxpy,{1,0,0},thcn2];

(*Create curved interior pole piece,positive x,negative y*)
flist={{1,2,3},{5,2,3,6},{4,1,3,6},{5,2,1,4},{5,4,6}};
cc1={length/2,-width/2,ctr-height/2};
cc2={length/2,-width/2,ctr+height/2};
Clear[polfcpxny];
polfcpxny[n_]:=(plist1={cc1+{radin*Sin[Pi*n/cN2/2],-radin*Cos[Pi*n/cN2/2],
              0},cc1+{radin*Sin[Pi*(n+1)/cN2/2],-radin*
                Cos[Pi*(n+1)/cN2/2],0},cc1};
      plist2={cc2+{radin*Sin[Pi*n/cN2/2],-radin*Cos[Pi*n/cN2/2],0},
          cc2+{radin*Sin[Pi*(n+1)/cN2/2],-radin*Cos[Pi*(n+1)/cN2/2],0},cc2};
      plist=Flatten[{plist1,plist2},1];
      radObjPolyhdr[plist,flist,{0,0,0}]);
Clear[Polepxny];
Polepxny=radObjCnt[Table[polfcpxny[n],{n,0,cN2-1}]];
radObjDrwAtr[Polepxny,{1,0,0},thcn2];

(*Create curved interior pole piece,positive x,positive y*)
flist={{1,2,3},{5,2,3,6},{4,1,3,6},{5,2,1,4},{5,4,6}};
cc1={length/2,width/2,ctr-height/2};
cc2={length/2,width/2,ctr+height/2};
Clear[polfcpxpy];
polfcpxpy[n_]:=(plist1={cc1+{radin*Sin[Pi*n/cN2/2],radin*Cos[Pi*n/cN2/2],0},
          cc1+{radin*Sin[Pi*(n+1)/cN2/2],radin*Cos[Pi*(n+1)/cN2/2],0},cc1};
      plist2={cc2+{radin*Sin[Pi*n/cN2/2],radin*Cos[Pi*n/cN2/2],0},
          cc2+{radin*Sin[Pi*(n+1)/cN2/2],radin*Cos[Pi*(n+1)/cN2/2],0},cc2};
      plist=Flatten[{plist1,plist2},1];
      radObjPolyhdr[plist,flist,{0,0,0}]);
Clear[Polepxpy];
Polepxpy=radObjCnt[Table[polfcpxpy[n],{n,0,cN2-1}]];
radObjDrwAtr[Polepxpy,{1,0,0},thcn2];

(*Create interior rectangular pole piece,central, positive x*)
Polemid0px=radObjRecMag[{(length)/4+radin/2,0.,ctr},{(length)/2+radin,width,height},{0,0,0}];
radObjDivMag[Polemid0px,nsub2];
radObjDrwAtr[Polemid0px,{1,0,0},thcn2];

(*Create interior rectangular pole piece,central, negative x*)
Polemid0nx=radObjRecMag[{-(length)/4-radin/2,0.,ctr},{(length)/2+radin,width,height},{0,0,0}];
radObjDivMag[Polemid0nx,nsub2];
radObjDrwAtr[Polemid0nx,{1,0,0},thcn2];

(*Create interior rectangular pole piece, edge,positive x, negative y*)
Polemidpxny=radObjRecMag[{(length)/4,-width/2.0-radin/2,ctr},{(length)/2,radin,height},{0,0,0}];
radObjDivMag[Polemidpxny,nsub2];
radObjDrwAtr[Polemidpxny,{1,0,0},thcn2];

(*Create interior rectangular pole piece, positive x, positive y *)
Polemidpxpy=radObjRecMag[{(length)/4,width/2.0+radin/2,ctr},{(length)/2,radin,height},{0,0,0}];
radObjDivMag[Polemidpxpy,nsub2];
radObjDrwAtr[Polemidpxpy,{1,0,0},thcn2];

(*Create interior rectangular pole piece, negative x,negative y*)
Polemidnxny=radObjRecMag[{-(length)/4,-width/2.0-radin/2,ctr},{(length)/2,radin,height},{0,0,0}];
radObjDivMag[Polemidnxny,nsub2];
radObjDrwAtr[Polemidnxny,{1,0,0},thcn2];

(*Create interior rectangular pole piece, positive y *)
Polemidnxpy=radObjRecMag[{-(length)/4,width/2.0+radin/2,ctr},{(length)/2,radin,height},{0,0,0}];
radObjDivMag[Polemidnxpy,nsub2];
radObjDrwAtr[Polemidnxpy,{1,0,0},thcn2];

(*Create base plate*)
Baseheight=145-gap2-height; 
Polebase=radObjRecMag[{0.,0.,ctr+height/2+Baseheight/2},{length+2*radin,2*radout+width,Baseheight},{0,0,0}];
radObjDivMag[Polebase,nsubb2];
radObjDrwAtr[Polebase,{1,0,0},thcn2];
(*Create the full pole configuration*)

Poles0=radObjCnt[{Polenxny,Polenxpy,Polepxny,Polepxpy,Polemid0px,Polemid0nx,Polemidpxpy,Polemidpxny,Polemidnxpy,Polemidnxny,Polebase}];
        
(* Iron material properties *)
(* ironmat=radMatSatIso[{20000,2},{0.1,2},{0.1,2}]; *)
ironmat=RadMatXc06[];
(*ironmatopera=radMatSatIso[{1.24649,0.34022,0.572764},{2699.29,118.794,14.817}];*)  (*gnuplot fit*)
ironmatopera=radMatSatIso[{1.43196,0.580746,1.25843},{2659.05,29.8208,1.19571}];  (*mathematica fit*)
(* Apply the iron material properties *)
radMatApl[Poles0,ironmat];

Return[Poles0])

CenterPole = createapole[Coillength,Coilwidth,Coilctr,Coilheight,Coilradin,Coilradout,Chamfer,Chamferdepth,Gap,cN,nsub,nsubb,n1,c1,thcn];
EndPole = createapole[Coillength,Endcoilwidth,Coilctr,Coilheight,Endcoilradin,Endcoilradout,Chamfer,Chamferdepth,Gap,cN,nsub,nsubb,n1,c1,thcn];
NearEndPole = createapole[Coillength,Nearendcoilwidth,Coilctr,Coilheight,Nearendcoilradin,Nearendcoilradout,Chamfer,Chamferdepth,Gap,cN,nsub,nsubb,n1,c1,thcn];

shift={-11,-9,-7,-5,-3,-1,1,3,5,7,9,11};
Pole={0,0,0,0,0,0,0,0,0,0,0,0};
npoles=12;

Do[Pole[[i]]=radObjDpl[CenterPole];
poshift=radTrfTrsl[{0,shift[[i]]*(Coilradout+Coilwidth/2),0}];
radTrfOrnt[Pole[[i]],poshift];,{i,3,npoles-2}]

Pole[[1]]=radObjDpl[EndPole];
Pole[[2]]=radObjDpl[NearEndPole];
Pole[[npoles-1]]=radObjDpl[NearEndPole];
Pole[[npoles]]=radObjDpl[EndPole];
nposhiftend=radTrfTrsl[{0,-((npoles-4)*(Coilradout+Coilwidth/2)+2*(Nearendcoilradout+Nearendcoilwidth/2)+1*(Endcoilradout+Endcoilwidth/2)),0}];
nposhiftnearend=radTrfTrsl[{0,-((npoles-4)*(Coilradout+Coilwidth/2)+1*(Nearendcoilradout+Nearendcoilwidth/2)),0}];
pposhiftend=radTrfTrsl[{0,+((npoles-4)*(Coilradout+Coilwidth/2)+2*(Nearendcoilradout+Nearendcoilwidth/2)+1*(Endcoilradout+Endcoilwidth/2)),0}];
pposhiftnearend=radTrfTrsl[{0,+((npoles-4)*(Coilradout+Coilwidth/2)+1*(Nearendcoilradout+Nearendcoilwidth/2)),0}];
radTrfOrnt[Pole[[1]],nposhiftend];
radTrfOrnt[Pole[[2]],nposhiftnearend];
radTrfOrnt[Pole[[npoles-1]],pposhiftnearend];
radTrfOrnt[Pole[[npoles]],pposhiftend];

Poles=radObjCnt[Pole];

Coil={0,0,0,0,0,0,0,0,0,0,0,0};
sign={1,-1,1,-1,1,-1,1,-1,1,-1,1,-1};
factor={1,1,1,1,1,1,1,1,1,1,1,1};
j={j1end,j1,j1,j1,j1,j1,j1,j1,j1,j1,j1,j1end};

Do[Coil[[i]]=radObjRaceTrk[{0.,shift[[i]]*(Coilradout+Coilwidth/2),Coilctr},{Coilradin,Coilradout},{Coillength,Coilwidth},Coilheight,n1,sign[[i]]*j[[i]]*factor[[i]]];
radObjDrwAtr[Coil[[i]],c1,thcn];,{i,3,npoles-2}]

Coil[[1]]=radObjRaceTrk[{0.,-((npoles-4)*(Coilradout+Coilwidth/2)+2*(Nearendcoilradout+Nearendcoilwidth/2)+1*(Endcoilradout+Endcoilwidth/2)),Coilctr},{Endcoilradin,Endcoilradout},{Coillength,Endcoilwidth},Coilheight,n1,sign[[1]]*j[[1]]*factor[[1]]];
Coil[[2]]=radObjRaceTrk[{0.,-((npoles-4)*(Coilradout+Coilwidth/2)+1*(Nearendcoilradout+Nearendcoilwidth/2)),Coilctr},{Nearendcoilradin,Nearendcoilradout},{Coillength,Nearendcoilwidth},Coilheight,n1,sign[[2]]*j[[2]]*factor[[2]]];
Coil[[npoles-1]]=radObjRaceTrk[{0.,+((npoles-4)*(Coilradout+Coilwidth/2)+1*(Nearendcoilradout+Nearendcoilwidth/2)),Coilctr},{Nearendcoilradin,Nearendcoilradout},{Coillength,Nearendcoilwidth},Coilheight,n1,sign[[npoles-1]]*j[[npoles-1]]*factor[[npoles-1]]];
Coil[[npoles]]=radObjRaceTrk[{0.,+((npoles-4)*(Coilradout+Coilwidth/2)+2*(Nearendcoilradout+Nearendcoilwidth/2)+1*(Endcoilradout+Endcoilwidth/2)),Coilctr},{Endcoilradin,Endcoilradout},{Coillength,Endcoilwidth},Coilheight,n1,sign[[npoles]]*j[[npoles]]*factor[[npoles]]];
radObjDrwAtr[Coil[[1]],c1,thcn];
radObjDrwAtr[Coil[[2]],c1,thcn];
radObjDrwAtr[Coil[[npoles-1]],c1,thcn];
radObjDrwAtr[Coil[[npoles]],c1,thcn];

Coils=radObjCnt[Coil];

(* Complete magnet *)    
Wig=radObjCnt[{Poles, Coils}];
SinglePole=radObjCnt[{Pole[[4]],Coil[[4]]}];
SinglePoleC=radObjCnt[{Coil[[4]]}];
SinglePoleP=radObjCnt[{Pole[[4]]}];

(* Define Mirror Coils *)
RadTrfZerPara[Wig,{0,0,0},{0,0,1}];
RadTrfZerPara[SinglePole,{0,0,0},{0,0,1}];
RadTrfZerPara[SinglePoleC,{0,0,0},{0,0,1}];
RadTrfZerPara[SinglePoleP,{0,0,0},{0,0,1}];

Plot the Geometry

These instructions plot the geometry of coils previously created.

(* Special Options for the Plot[]*)
RadPlot3DOptions[];

(* Save the Geometry in "dr" *)
dr=radObjDrw[SinglePole];

(* Draw the Geometry *)
Show[Graphics3D[dr
,PlotLabel->FontForm["Single Pole of Wiggler",{"Times",14}]
,DefaultFont -> {"Times",10}]];

Show[Graphics3D[dr
,PlotLabel->FontForm["Single Pole of Wiggler",{"Times",14}]
,DefaultFont -> {"Times",10},ViewPoint->{0,-2000,0}]];

(* Special Options for the Plot[]*)
RadPlot3DOptions[];

(* Save the Geometry in "dr" *)
dr2=radObjDrw[Wig];

(* Draw the Geometry *)
Show[Graphics3D[dr2
,PlotLabel->FontForm["12 Pole Wiggler",{"Times",14}]
,DefaultFont -> {"Times",10}]];

Solve the problem

t1 = 
        AbsoluteTime[];
res = RadSolve[Wig, 0.00001, 150000];
Print["Maximum Magnetic Field anywhere in Object is ",  
      res[[2]],  " T"]
t2 = AbsoluteTime[];
Print["Solved in ", Round[t2 - t1], " seconds"]

lpos = 
      80 ;    
peak = 
      radFld[Wig, "bz", {0, lpos, 0}];
db = 
        radFld[Wig, "bz", {10, lpos, 0}];
dbb = Abs[peak - db]/peak;
Print["Peak field = ", peak, " T"]
Print["Vertical field horizontal rolloff = ", 
      dbb*100, " %, at x=10mm"]

list = 
        radFldPtcTrj[Wig, 5.00, {0. , 0. , 0. , 0. }, {-1040. , 1040. }, 
          500];
ListPlot[list[[All, 2]]];
xout = list[[500, 2]];
pxout = list[[500, 3]];
Print["Trajectory offset at exit = ",  
      xout,  " mm"]

RadPlotOptions[];
p1 = 
        ListPlot[list[[All, 2]], 
          PlotStyle -> {Thickness[0.02], RGBColor[1, 0, 0]}, 
          FrameLabel -> {"Z,steps", "X,mm", "12 Poles, Simple \
Optimized", ""}, 
          RotateLabel -> True];
SetDirectory["/home/jtu2/raid/redesign/radia/addpoles/"];
Export["simple12_poles2_per32cm_gap_trajectory.eps", 
        p1];

RadPlotOptions[];
p1 = 
        Plot[radFld[Wig, "bz", {x, 100, 0}], {x, -100, 100}, 
          PlotRange -> {1, 2}, 
          PlotStyle -> {Thickness[0.02], RGBColor[1, 0, 0]}, 
          FrameLabel -> {"X,mm", "Bz,Tesla", "12 Poles, Simple \
Optimized", ""}, 
          RotateLabel -> True];
Export["simple12_poles2_per32cm_gap_profile.eps", 
        p1];

RadPlotOptions[];
p1 = 
        Plot[100*(Abs[peak - radFld[Wig, "bz", {x, 100, 0}]])/
              peak, {x, -10, 10}, PlotRange -> {- .01, 0.1}, 
          PlotStyle -> {Thickness[0.02], RGBColor[1, 0, 0]}, 
          FrameLabel -> {"X,mm", "dBB,%", "12 Poles, Simple \
Optimized", ""}, 
          RotateLabel -> True];
Export["simple12_poles2_per32cm_gap_dbb.eps", 
        p1];

tmp3 = 
        OpenWrite["simple12_poles2_per32cm_gap_field.dat", 
          FormatType -> OutputForm];
Do[Write[tmp3, NumberForm[i/5. ], "          ", 
        NumberForm[radFld[Wig, "bz", {i/5. , 100, 0}]], "         ", 
        FortranForm[
          Abs[peak - radFld[Wig, "bz", {i/5. , -100, 0}]]/
            peak]], {i, -500, 500}]
Close[tmp3]

Plot the Magnetic Field

The magnetic field magnitude inside the wiggler  is plotted along the z axis.

RadPlotOptions[];
Plot[radFld[Wig,"bz",{0,y,0}],{y,-1200,0}
,FrameLabel->{"Z [mm]","By[T]","Y=X=0",""}];

RadPlotOptions[];
Plot[radFld[Wig,"bz",{0,y,0}],{y,0,1200}
,FrameLabel->{"Z [mm]","By[T]","Y=X=0",""}];

RadPlotOptions[];
Plot[radFld[Wig,"bz",{0,y,0}],{y,30,170}
,FrameLabel->{"Z [mm]","By[T]"," Y=X=0",""}];

RadPlotOptions[];
Plot[radFld[Wig,"bz",{x,100,0}],{x,-150,150}
,FrameLabel->{"X [mm]","By[T]"," Z=Y=0",""}];

RadPlotOptions[];
Plot[radFld[Wig,"bz",{x,100,0}],{x,0,50}
,FrameLabel->{"X [mm]","By[T]"," Z=Y=0",""}];

Calculate  and  Write out  of  Magnetic Field Map

\( (*\ \ \ Format\ \ \  : \ \ x\ \ y\ \ z\ \ Bx\ \ By\ \ Bz\ \ \ *) 

    Clear[mtFld]\)
mtFld[x_, y_, z_] := radFld[Wig, "b", {x, y, z}];

t0 = AbsoluteTime[];
bzmap0 = 
        Table[{x, y, z, mtFld[x, y, z]}, {x, 0, 48, 4}, {z, 0, 26, 2}, {y, 0, 
            1040, 2}];
huhu1 = Flatten[bzmap0];
bzmap1 = Partition[huhu1, 6];
t1 = AbsoluteTime[];
ddd = Date[];
\(\(Print[" ", ddd[\([3]\)], ".", ddd[\([2]\)], ".", 
        ddd[\([1]\)], "   ", 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ddd[\([4]\)], ":", ddd[\([5]\)], ":", 
        ddd[\([6]\)], 
"
 Solved  in : ", 
        Round[t1 - t0], " sec"];\)\)

SetDirectory["/home/jtu2/raid/redesign/radia/addpoles/"];
tmp2 = 
        OpenWrite["simple12polemapfromradia_poles2_per32cm_gap.dat", 
          FormatType -> OutputForm];
tab = bzmap1;
SetOptions[tmp2, PageWidth -> 200];
Write[tmp2, "# Full Gap = ", 
        2*Gap,  " mm"];
Write[tmp2, "# Pole Width = ",  
        Coillength + 2*Coilradin,  " mm"];
Write[tmp2, "# Current = ",  
        Current,  " Amps"];
Write[tmp2, "# End Turns = ",  
        Endturns];
Write[tmp2, "# Peak Field = ",  
        peak,  " T"];
Write[tmp2, "# db/b @ x=10mm = ",  
        100*dbb,  " %"];
Write[tmp2, "# xout = ",  xout*1000,  " um"];
\(\(Do[
        Write[tmp2, 
          NumberForm[ .1*
              tab[[i, 1]]], "     ", 
          NumberForm[ .1*
              tab[[i, 3]]], "     ", 
          NumberForm[ .1*
              tab[[i, 
                2]]], "     ", 
		FortranForm[
            10000*tab[[i, 
                4]]], "     ", 
		FortranForm[
            10000*tab[[i, 
                6]]], "     ", 
		FortranForm[
            10000*tab[[i, 5]]]], {i, 
          Length[tab]}];\)\)
Close[tmp2]