%
function cavitywake_new(wfilename,nslices,a,g,l,delta,delta_g)
    sstep=l/(nslices-1);
    
    wfid=fopen(wfilename,"wt");
    
    for s=0.0:sstep:l
        wl=w_long(s,a,g,l,delta,delta_g)*l;
        wt=w_trans(s,a,g,l,delta,delta_g)*l;
        wx=wt;
        wy=wt;
        fprintf(wfid,"%22.18e %22.18e %22.18e %22.18e\n",s/299792458,wl,wx,wy);
    end;
    fclose(wfid);
endfunction;


%
% longitudinal wakefield
% return value is in V/Cm
%
function wl=w_long(s,a,g,l,delta,delta_g)
    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*299792458*tmp/4.0;
endfunction;

%
% transverse wakefield
% return value is in V/Cm^2
%
function wt=w_trans(s,a,g,l,delta,delta_g)
    gx=g;
    s0=0.169*a^1.79*gx^0.38*l^-1.17;
    tmp=s0/(pi*a^4)*(1.0-(1.0+sqrt(s/s0))*exp(-sqrt(s/s0)));

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

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

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

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

    wt=0.25*4.0*377.0*299792458*tmp;
endfunction;

