AssignFieldErr

 function r=AssignFieldErr(r,refpos,N,rho,BNn,ANn)
 
 r : at lattice
 refpos : index of elements for wich the field error has to be applied
 N main component (2=quadrupole,3=sextupole)
 rho: reference radius

 the field errors are defined by the magnet designer with the multipole
 expansion
 B=B_N(rho0)*sum_(n=1)^(infty)((BNn+i*ANn)*(z/rho0)^(n-1))

 and by AT as
 B=Brho*sum_(n=1)^(infty)((b_n+i*a_n)*(z)^(n-1))
 
 the input BNn and ANn are the normal and skew field errors at rho0. 
    as defined by:
 B^(N)_(n) = radius^(n-N)*b_n/b_N

 b_n=radius^(N-n)*b_N*B^(N)_(n)

 optional output ,Pnew,Panew the PolynomB and PolynomA set in AT.

0001 function [r,Pnew,Panew]=AssignFieldErr(r,refpos,N,rho,BNn,ANn)
0002 % function r=AssignFieldErr(r,refpos,N,rho,BNn,ANn)
0003 %
0004 % r : at lattice
0005 % refpos : index of elements for wich the field error has to be applied
0006 % N main component (2=quadrupole,3=sextupole)
0007 % rho: reference radius
0008 %
0009 % the field errors are defined by the magnet designer with the multipole
0010 % expansion
0011 % B=B_N(rho0)*sum_(n=1)^(infty)((BNn+i*ANn)*(z/rho0)^(n-1))
0012 %
0013 % and by AT as
0014 % B=Brho*sum_(n=1)^(infty)((b_n+i*a_n)*(z)^(n-1))
0015 %
0016 % the input BNn and ANn are the normal and skew field errors at rho0.
0017 %    as defined by:
0018 % B^(N)_(n) = radius^(n-N)*b_n/b_N
0019 %
0020 % b_n=radius^(N-n)*b_N*B^(N)_(n)
0021 %
0022 % optional output ,Pnew,Panew the PolynomB and PolynomA set in AT.
0023 %
0024 
0025 
0026 if nargin<6
0027     ANn=BNn*0;
0028 end
0029 P=getcellstruct(r,'PolynomB',refpos,1,N);
0030 Pa=getcellstruct(r,'PolynomA',refpos,1,N);
0031 
0032 if N==1 && P(1)==0
0033     % is a dipole.
0034     P=getcellstruct(r,'BendingAngle',refpos);% ./getcellstruct(r,'Length',refpos);
0035 end
0036 
0037 if length(rho)==1
0038     Pnew=repmat(rho.^(N-[1:length(BNn)]),length(P),1).*repmat(P,1,length(BNn)).*repmat(BNn,length(P),1);
0039     Panew=repmat(rho.^(N-[1:length(ANn)]),length(Pa),1).*repmat(P,1,length(ANn)).*repmat(ANn,length(Pa),1); % refer to normal component also skew
0040 elseif length(rho)==size(P,1)
0041     RNB=[];
0042     RNA=[];
0043     
0044     for irrh=1:size(P,1)
0045         RNB(irrh,:)=rho(irrh).^(N-[1:length(BNn)]);
0046         RNA(irrh,:)=rho(irrh).^(N-[1:length(ANn)]);
0047     end
0048      
0049     Pnew=RNB.*repmat(P,1,length(BNn)).*repmat(BNn,length(P),1);
0050     Panew=RNA.*repmat(P,1,length(ANn)).*repmat(ANn,length(Pa),1);
0051 end
0052 
0053 for ir=1:length(refpos)
0054     pbold=r{refpos(ir)}.PolynomB;
0055     paold=r{refpos(ir)}.PolynomA;
0056    
0057     if ir==1 && length(pbold)>length(Pnew) % pad with zeros
0058         %disp('padding')
0059         Pnew=[Pnew zeros(length(refpos),length(pbold)-length(Pnew))];
0060         Panew=[Panew zeros(length(refpos),length(pbold)-length(Panew))];
0061     end
0062     
0063 %     size(pbold)
0064 %     size(paold)
0065 %     size(Pnew)
0066 %     size(Panew)
0067 
0068     Pnew(ir,1:length(pbold))=Pnew(ir,1:length(pbold))+pbold;
0069     Panew(ir,1:length(paold))=Panew(ir,1:length(paold))+paold;
0070     
0071     
0072     r{refpos(ir)}.PolynomB=Pnew(ir,:);
0073     r{refpos(ir)}.PolynomA=Panew(ir,:); % update both A and B fields!
0074     r{refpos(ir)}.MaxOrder=size(Pnew,2)-1; % -1 for c-integrators consitency
0075 end
0076 
0077 
0078 return