function [Phib,rhotor,psi] = calcPhi(psiN,iota,psib,psia)

%% iota-based Phi (numerically best?)
Phi = zeros(numel(psiN),size(iota,2));
psi = bsxfun(@plus,psiN*(psib-psia),psia);
for ir=1:(size(Phi,1)-1)
    Phi(ir+1,:) = Phi(ir,:) + ...
        2*(psi(ir+1,:)-psi(ir,:))./(iota(ir+1,:)+iota(ir,:));
    % recursive equation obtained from trapeze integration rule
end
rhotor = sqrt(bsxfun(@times,Phi,1./Phi(end)));
Phib = Phi(end, :);


%Interpos integration. This does not make much difference
% [~,~,~,Phi] = interpos( psi, 1./iota);
% rhotor = sqrt(bsxfun(@times,Phi,1./Phi(end)));
% Phib = Phi(end,:);

return