%
% this file contains all procedures required to calculate longitudinal and transverse wake fields
%
% rfparams=[gradient,phase,wavelength,a,g,l,delta,delta_g]
% units=[eV/m,Degree,m,m,m,m,?,?]
%

%-----------------------------------------------------------------------------------

%
% transverse wakefield
% s is given in micro meters
% return value is in V/pCm^2
%
%function wt=w_transv(s,a,g,l)
%    s0=0.169*a^(1.79)*g^(0.38)*l^(-1.17);
%    wt=4.0*377.0*3.0e8*s0*1.0e-12/(pi*a^4)*(1.0-(1.0+sqrt(s*1.0e-6/s0))*exp(-sqrt(s*1.0e-6/s0)));
%endfunction;

%function wt=w_transv2(s,a,g,l)
%    tmp=g/l;
%    alpha=1.0-0.4648*sqrt(tmp)-(1.0-2.0*0.4648)*tmp;
%    s0=g/8.0*(a/(l*alpha))^2;
%    wt=4.0*377.0*3.0e8*s0*1.0e-12/(pi*a^4)*(1.0-(1.0+sqrt(s*1.0e-6/s0))*exp(-sqrt(s*1.0e-6/s0)));
%endfunction;

%
% longitudinal wakefield
% s is given in micro meters
% return value is in V/pCm
%
%function wl=w_long(s,a,g,l)
%    tmp=g/l;
%    alpha=1.0-0.4648*sqrt(tmp)-(1.0-2.0*0.4648)*tmp;
%    s0=g/8.0*(a/(l*alpha))^2;
%    wl=377.0*3.0e8*1.0e-12/(pi*a*a)*exp(-sqrt(s*1.0e-6/s0));
%endfunction;


%
% routines to calculate the wakefields
% a file is created with name $filename that can be read when a beam is created
% z_list contains longitudinal positions and weights of the slices
%
% unit to store in the file is MV/pCm^2
% multiply charge with factor for pC and MV
%
function wtlist=cavity_wake_trans(charge,z_matrix,a,g,l)
tic
    tmp=g/l;
    alpha=1.0-0.4648*sqrt(tmp)-(1.0-2.0*0.4648)*tmp;
    s0=g/8.0*(a/(l*alpha))^2;

    zstep=abs(z_matrix(1,1)-z_matrix(2,1))*1.0e-6/s0;
    ll=length(z_matrix);

    tmp=zeros((ll+1)*ll/2,1);
    tc=0;
    wl=zeros(ll,1);
    for i=1:ll
        xx=sqrt((i-1)*zstep);
        wl(i)=1.0-(1.0+xx)*exp(-xx);
	for j=1:i
            tc=tc+1;
            tmp(tc)=wl(i-j+1);
	end;
    end;
    wtlist=1.602177e-7*charge*1.0e-6*4.0*377.0*3.0e8*s0*1.0e-12/(pi*a^4)*tmp;
    size(tmp)
toc
endfunction;

function wllist=cavity_wake_long(charge,z_matrix,a,g,l)
tic
    gdl=g/l;
    alpha=1.0-0.4648*sqrt(gdl)-(1.0-2.0*0.4648)*gdl;
    s0=g/8.0*(a/(l*alpha))^2;
    aa=a*a;

    z_matrix(:,2)=charge*z_matrix(:,2);
    zstep=abs(z_matrix(1,1)-z_matrix(2,1))*1e-6/s0;
    n=length(z_matrix);

    tmp=zeros(n,2);
    tc=0;
    wl=zeros(n,1);
    for i=1:n
        wl(i)=exp(-sqrt((i-1)*zstep));
	tsum=0.0;
	for j=1:(i-1)
	    tsum=tsum+z_matrix(i-j,2)*wl(j+1);
	end;
        tc=tc+1;
	tmp(tc,1)=z_matrix(i,1);
        tmp(tc,2)=tsum+0.5*z_matrix(i,2)*wl(1);
    end;
    tmp(:,2)=-377.0*3.0e8*1.0e-12/(pi*aa)*tmp(:,2)*1.602177e-7*1e-6;
    wllist=tmp;
    size(tmp)
toc
endfunction;

%-----------------------------------------------------------------------------------

function wt=w_trans_old(s,a,g,l,delta,delta_g)
    tmp=g/l;
    alpha=1.0-0.4648*sqrt(tmp)-(1.0-2.0*0.4648)*tmp;
    s0=g/8.0*(a/(l*alpha))^2;
    wt=4.0*377.0*3.0e8*s0*1.0e-12/(pi*a^4)*(1.0-(1.0+sqrt(s*1.0e-6/s0))*exp(-sqrt(s*1.0e-6/s0)));
endfunction;

function wt=w_trans_new(s,a,g,l,delta,delta_g)
    s=s*1e-6;
    
    ax=a;
    gx=g;
    s0=0.169*ax^1.79*gx^0.38*l^-1.17;
    ss0=sqrt(s/s0);
    tmp=s0/(pi*ax^4)*(1.0-(1.0+ss0)*exp(-ss0));
    
    ax=a*(1.0-0.5*delta);
    gx=g-0.5*delta_g;
    s0=0.169*ax^1.79*gx^0.38*l^-1.17;
    ss0=sqrt(s/s0);
    tmp=tmp+s0/(pi*ax^4)*(1.0-(1.0+ss0)*exp(-ss0));

    ax=a*(1.0+0.5*delta);
    gx=g+0.5*delta_g;
    s0=0.169*ax^1.79*gx^0.38*l^-1.17;
    ss0=sqrt(s/s0);
    tmp=tmp+s0/(pi*ax^4)*(1.0-(1.0+ss0)*exp(-ss0));
    
    ax=a*(1.0-delta);
    gx=g-delta_g;
    s0=0.169*ax^1.79*gx^0.38*l^-1.17;
    ss0=sqrt(s/s0);
    tmp=tmp+0.5*s0/(pi*ax^4)*(1.0-(1.0+ss0)*exp(-ss0));
    
    ax=a*(1.0+delta);
    gx=g+delta_g;
    s0=0.169*ax^1.79*gx^0.38*l^-1.17;
    ss0=sqrt(s/s0);
    tmp=tmp+0.5*s0/(pi*ax^4)*(1.0-(1.0+ss0)*exp(-ss0));

    wt=0.25*4.0*377.0*3e8*1e-12*tmp;
endfunction;

function wtlist=cavity_wake_trans_new(charge,z_matrix,a,g,l,delta,delta_g)
tic
    ll=length(z_matrix);

    tmp=zeros((ll+1)*ll/2,1);
    tc=0;
    for i=1:ll
        z0=z_matrix(i,1);
	for j=1:i
            tc=tc+1;
            z=z_matrix(j,1);
            tmp(tc)=w_trans_new(z0-z,a,g,l,delta,delta_g);
	end;
    end;
    wtlist=1.602177e-7*charge*1.0e-6*tmp;
    size(tmp)
toc
endfunction;

function wtlist=cavity_wake_trans_new2(charge,z_matrix,a,g,l,delta,delta_g)
tic
    ax1=a;
    gx1=g;
    s01=0.169*ax1^1.79*gx1^0.38*l^-1.17;
    s01b=s01/(pi*ax1^4);

    ax2=a*(1.0-0.5*delta);
    gx2=g-0.5*delta_g;
    s02=0.169*ax2^1.79*gx2^0.38*l^-1.17;
    s02b=s02/(pi*ax2^4);

    ax3=a*(1.0+0.5*delta);
    gx3=g+0.5*delta_g;
    s03=0.169*ax3^1.79*gx3^0.38*l^-1.17;
    s03b=s03/(pi*ax3^4);
    
    ax4=a*(1.0-delta);
    gx4=g-delta_g;
    s04=0.169*ax4^1.79*gx4^0.38*l^-1.17;
    s04b=0.5*s04/(pi*ax4^4);
    
    ax5=a*(1.0+delta);
    gx5=g+delta_g;
    s05=0.169*ax5^1.79*gx5^0.38*l^-1.17;
    s05b=0.5*s05/(pi*ax5^4);

    z_matrix(:,1)=z_matrix(:,1)*1e-6;
    
    ll=length(z_matrix);

    tmp=zeros((ll+1)*ll/2,1);
    tc=1;
    for i=1:ll
        z0=z_matrix(i,1);

        s=z0-z_matrix(1:i,1);
        ss01=sqrt(s/s01);
        ss02=sqrt(s/s02);
        ss03=sqrt(s/s03);
        ss04=sqrt(s/s04);
        ss05=sqrt(s/s05);
        tmp(tc:tc+(i-1))=(s01b*(1.0-(1.0+ss01).*exp(-ss01))+...
                          s02b*(1.0-(1.0+ss02).*exp(-ss02))+...
                          s03b*(1.0-(1.0+ss03).*exp(-ss03))+...
                          s04b*(1.0-(1.0+ss04).*exp(-ss04))+...
                          s05b*(1.0-(1.0+ss05).*exp(-ss05)));
        tc=tc+i;
    end;
    wtlist=1.602177e-7*charge*1.0e-6*0.25*4.0*377.0*3e8*1e-12*tmp;
    size(tmp)
toc
endfunction;




function wl=w_long_old(s,a,g,l,delta,delta_g)
    tmp=g/l;
    alpha=1.0-0.4648*sqrt(tmp)-(1.0-2.0*0.4648)*tmp;
    s0=g/8.0*(a/(l*alpha))^2;
    wl=377.0*3.0e8*1.0e-12/(pi*a*a)*exp(-sqrt(s*1.0e-6/s0));
endfunction;

function wl=w_long_new(s,a,g,l,delta,delta_g)
    s=s*1e-6;
    
    ax=a;
    gx=g;
    s0=0.41*ax^1.8*gx^1.6*l^-2.4;
    tmp=1/(pi*ax*ax)*exp(-sqrt(s/s0));

    ax=a*(1.0-0.5*delta);
    gx=g-0.5*delta_g;
    s0=0.41*ax^1.8*gx^1.6*l^-2.4;
    tmp=tmp+1/(pi*ax*ax)*exp(-sqrt(s/s0));

    ax=a*(1.0+0.5*delta);
    gx=g+0.5*delta_g;
    s0=0.41*ax^1.8*gx^1.6*l^-2.4;
    tmp=tmp+1/(pi*ax*ax)*exp(-sqrt(s/s0));

    ax=a*(1.0-delta);
    gx=g-delta_g;
    s0=0.41*ax^1.8*gx^1.6*l^-2.4;
    tmp=tmp+0.5/(pi*ax*ax)*exp(-sqrt(s/s0));

    ax=a*(1.0+delta);
    gx=g+delta_g;
    s0=0.41*ax^1.8*gx^1.6*l^-2.4;
    tmp=tmp+0.5/(pi*ax*ax)*exp(-sqrt(s/s0));

    wl=377.0*3e8*1e-12*tmp/4.0;
endfunction;

function wllist=cavity_wake_long_new(charge,z_matrix,a,g,l,delta,delta_g)
tic
    n=length(z_matrix);

    tmp=zeros(n,2);
    tc=0;
    for i=1:n
        z0=z_matrix(i,1);
        wgt0=charge*z_matrix(i,2);
	tsum=0.0;
	for j=1:(i-1)
            z=z_matrix(j,1);
            wgt=charge*z_matrix(j,2);
	    tsum=tsum+wgt*w_long_new(z0-z,a,g,l,delta,delta_g);
	end;
        tc=tc+1;
	tmp(tc,1)=z_matrix(i,1);
        tmp(tc,2)=tsum+0.5*wgt0*w_long_new(0,a,g,l,delta,delta_g);
    end;
    tmp(:,2)=-tmp(:,2)*1.602177e-7*1e-6;
    wllist=tmp;
    size(tmp)
toc
endfunction;

function wllist=cavity_wake_long_new2(charge,z_matrix,a,g,l,delta,delta_g)
tic
    ax1=a;
    ax1b=1/(pi*ax1*ax1);
    gx1=g;
    s01=0.41*ax1^1.8*gx1^1.6*l^-2.4;

    ax2=a*(1.0-0.5*delta);
    ax2b=1/(pi*ax2*ax2);
    gx2=g-0.5*delta_g;
    s02=0.41*ax2^1.8*gx2^1.6*l^-2.4;

    ax3=a*(1.0+0.5*delta);
    ax3b=1/(pi*ax3*ax3);
    gx3=g+0.5*delta_g;
    s03=0.41*ax3^1.8*gx3^1.6*l^-2.4;

    ax4=a*(1.0-delta);
    ax4b=0.5/(pi*ax4*ax4);
    gx4=g-delta_g;
    s04=0.41*ax4^1.8*gx4^1.6*l^-2.4;

    ax5=a*(1.0+delta);
    ax5b=0.5/(pi*ax5*ax5);
    gx5=g+delta_g;
    s05=0.41*ax5^1.8*gx5^1.6*l^-2.4;

    wl0=0.5*(ax1b+ax2b+ax3b+ax4b+ax5b);
    
    z_matrix(:,1)=z_matrix(:,1)*1e-6;
    
    n=length(z_matrix);
    tmp=zeros(n,2);
    tc=0;
    for i=1:n
        s=z_matrix(i,1)-z_matrix(1:i-1,1);
        tsum=z_matrix(1:i-1,2)'*...
	     (ax1b*exp(-sqrt(s/s01))+...
	      ax2b*exp(-sqrt(s/s02))+...
	      ax3b*exp(-sqrt(s/s03))+...
	      ax4b*exp(-sqrt(s/s04))+...
	      ax5b*exp(-sqrt(s/s05)));
        tc=tc+1;
        tmp(tc,2)=tsum+z_matrix(i,2)*wl0;
    end;
    tmp(:,1)=z_matrix(:,1)*1e6;
    tmp(:,2)=-tmp(:,2)*charge*1.602177e-7*1e-6*377.0*3e8*1e-12/4.0;
    wllist=tmp;
    size(tmp)
toc
endfunction;

%-----------------------------------------------------------------------------------

function cavity_wake(filename,charge,sigmamin,sigmamax,sigmaz,nslice,a,g,l,delta,delta_g)
% create a list of longitudinal and transverse wake kicks
% required by the InjectorBeam command
    fid=fopen(filename,"rt");
    if fid~=-11234
        fid=fopen(filename,"wt");
        fprintf(fid,"%d\n",nslice);
        
        gausslong=gausslist(sigmamin,sigmamax,sigmaz,1.0,5*nslice);
        wakelong=cavity_wake_long_new2(charge,gausslong,a,g,l,delta,delta_g);
        
        gausstrans=gausslist(sigmamin,sigmamax,sigmaz,1.0,nslice);
        waketrans=cavity_wake_trans_new2(charge,gausstrans,a,g,l,delta,delta_g);
        wl=wakelong(floor(nslice/2)*5+3,2);
        fprintf(fid,"%22.18e\n",wl);
        for i=1:nslice
            z=wakelong(3+5*(i-1),1);
            wl=wakelong(3+5*(i-1),2);
            wgt=gausstrans(i,2);
	    fprintf(fid,"%20.16e %20.16e %20.16e\n",z,wl,wgt);
        end;
        
        nt=length(waketrans);
        for i=1:nt
    	    fprintf(fid,"%20.16e\n",waketrans(i));
        end;
        fclose(fid);
    end;
endfunction;


function cavity_wake_zero(filename,sigmamin,sigmamax,sigmaz,nslice)
% create a list of longitudinal and transverse wake kicks
% required by the InjectorBeam command
% all kicks are set to 0
    fid=fopen(filename,"wt");
    fprintf(fid,"%d\n",nslice);
    fprintf(fid,"1.0\n");

    gl=gausslist(sigmamin,sigmamax,sigmaz,1.0,nslice);
    n=length(gl);
    for i=1:n
        fprintf(fid,"%20.16e %20.16e %20.16e\n",gl(i,1),0.0,gl(i,2));
    end;
    maxi=(nslice+1)*nslice/2;
    for i=1:maxi
        fprintf(fid,"%20.16e\n",0.0);
    end;
    fclose(fid);
endfunction;


function gl=gausslist(sigmamin,sigmamax,sigmaz,charge,nslices)
    minz=sigmamin*sigmaz;
    maxz=sigmamax*sigmaz;
    zstep=(maxz-minz)/nslices;
    minz=minz+zstep/2.0;
    maxz=maxz-zstep/2.0;
    
    for i=1:nslices
        z=minz+(i-1)*zstep;
        gl(i,1)=z;
        gl(i,2)=charge*0.5*(erf((z+zstep/2)/(sigmaz*sqrt(2)))-erf((z-zstep/2)/(sigmaz*sqrt(2))));
    end;
endfunction;
