function [shear,dshear_dx] = eval_shear(x,g,v,model,gaussgrid)
% function [iota,dshear_dx] = eval_shear(x,g,v,model,whichgrid)
% returns shear = rho/q dq/drho with rho=rhotorN
% #codegen

if ~gaussgrid;
    gauss = false;
    ngrid = model.rgrid.nrho;
else
    gauss = true;
    ngrid = model.rgrid.nrhogauss;
end

if nargout>1
        der = true; % compute derivative
        nx = size(x, 2);
        dshear_dx = zeros(ngrid,size(x, 1), nx); % init
else 
        der = false;
end

%shear = -2*V./(Vprime) .* diota_drho./iota;
% = -2*V./(Vprime) .* d/drho(1/rho dpsi/drho) ./ (1/rho dpsi/drho)
% note that the Phib factors cancel

%shear = -rhotor .* diota_drho./iota;
% = -rhotor .* d/drho(1/rho dpsi/drho) ./ (1/rho dpsi/drho)
% note that the Phib factors cancel
if ~gauss
    % avoid singularity at rho=0 by using l'hopital's rule for
    % first rho point.
    diota_drho = model.psi.Lamppr * x(model.psi.xind,:); % actually 2*Phib
    iota       = model.psi.Lampor * x(model.psi.xind,:); % actually 2*Phib
    shear =  - bsxfun(@times, model.rgrid.rho, diota_drho./iota);
    
    if der
       % dshear/dPsihat = -s/(diota/rho) *Ciotaprime + s/iota* Ciota
       for xi = 1:nx
           dshear_dpsihat = bsxfun(@times,shear(:,xi)./diota_drho(:,xi),model.psi.Lamppr)...
                - bsxfun(@times,shear(:,xi)./iota(:,xi),model.psi.Lampor);
           dshear_dx(:,model.psi.xind,xi) = dshear_dpsihat;   
       end
    end
else
    % all points at rho>0, no need for tricks
    %%
    diota_drho = model.psi.Lampprgauss * x(model.psi.xind,:); % actually 2*Phib
    iota       = model.psi.Lamporgauss * x(model.psi.xind,:); % actually 2*Phib
    shear = - model.rgrid.rhogauss.*diota_drho./iota;
    
    if der
        % dshear/dPsihat = -s/(diota/rho) *Ciotaprime + s/iota* Ciota
        dshear_dpsihat = bsxfun(@times,shear./diota_drho,model.psi.Lampprgauss)...
            - bsxfun(@times,shear./iota,model.psi.Lamporgauss);
        dshear_dx(:,model.psi.xind,:) = dshear_dpsihat;
    end
end



return
