function fv=evaluateCenters(globalPar, par, cax)
    %par is a nfunc x npar matrix of coefficients
    %globalPar is a row vector of global coefficients
    %cax is a column vector of time instances
    
    oneVector=ones(1,size(cax,1));
    cx2=ones(1,size(par,1))'*cax';
    oneMatrix=ones(size(cx2));
    
    A=globalPar(1);
    tau=globalPar(2);
    alpha=globalPar(3);
    
    t0=globalPar(5);
    t1=globalPar(6);
    
    k1=par(:,1)*oneVector;
    B=par(:,2)*oneVector;
    k2=par(:,3)*oneVector;
    dt=par(:,4)*oneVector;
    
    kt=k1./(k2-oneMatrix/tau);
    ka=k1./(k2-oneMatrix*alpha);
    cx2p=cx2-dt;
    
    e0=exp(-cx2p/tau);
    e1=exp(-cx2p*alpha);
    ek=exp(-cx2p.*k2);
    
    r0=zeros(size(cx2));
    if alpha==1/tau
        fv=e0.*cx2p;
        rsel=oneMatrix/tau==k2;
        r0(rsel)=k1(rsel).*cx2p(rsel).*cx2p(rsel)*ek(rsel);
        r0(~rsel)=kt.*(fv-(e0-ek)./(k2-oneMatrix/tau));
        
    else
        q=tau/(alpha*tau-1);
        fv=q*(e0-e1);
        rsel=oneMatrix/tau==k2;
        r0(rsel)=r0(rsel)+q*k1(rsel).*cx2p(rsel);
        r0(~rsel)=r0(~rsel)+q*kt(~rsel).*(e0(~rsel)-ek(~rsel));
        rsel=alpha==k2;
        r0(rsel)=r0(rsel)-q*k1(rsel).*cx2p(rsel);
        r0(~rsel)=r0(~rsel)-q*ka(~rsel).*(e1(~rsel)-ek(~rsel));
        
    B1=oneMatrix-B;
    
    fv=B.*fv+B1.*r0;
   
    fv=fv*A;
    sel=cx2p<0;
    fv(sel)=0;
    

end

%par=ones(2,