%
% this file contains all procedures required to calculate longitudinal and transverse wake fields
%
% rfparams=[gradient,phase,wavelength,a,g,l]
% 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));


%
% 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 cavity_wake(filename,charge,sigmamin,sigmamax,sigmaz,nslice,a,g,l)
% create a list of longitudinal and transverse wake kicks
% required by the InjectorBeam command
    fid=fopen(filename,"rt");
    if fid==-1
        fid=fopen(filename,"wt");
        fprintf(fid,"%d\n",nslice);
        
        gausslong=gausslist(sigmamin,sigmamax,sigmaz,1.0,5*nslice);
        wakelong=cavity_wake_long(charge,gausslong,a,g,l);
        
        gausstrans=gausslist(sigmamin,sigmamax,sigmaz,1.0,nslice);
        waketrans=cavity_wake_trans(charge,gausstrans,a,g,l);
        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;
