function [pec_pergyro,jec_pergyro] = echcd_per_gyro(stap,geop,u,model,params)

assert(iscolumn(u),'RAPTOR:WronguSize','u must be a column vector');

nunits = numel(params.echcd.rdep);

pec_pergyro = zeros(model.rgrid.nrhogauss,nunits);
jec_pergyro = zeros(model.rgrid.nrhogauss,nunits); 

iupower = params.echcd.uindices(1:nunits);

for ii=1:nunits
    u_only_power_i = u;
    % set all power to zero, except one
    u_only_power_i(iupower) = 0;
    iup = iupower(ii);
    u_only_power_i(iup) = u(iup);

    [pec,jec,~,~,~,~] = echcd_gaussian(stap,geop,u_only_power_i,model,params.echcd); % this gets the jecperunit
    
    pec_pergyro(:,ii) = pec;
    jec_pergyro(:,ii) = jec;
end

end