function geom = prep_geom(rhoin,rho,rhogauss)
% function geom = prep_geom(rhoin,rho,rhogauss)
% prepare spline matrices for geometry-related profiles
%
% rhoin: grid of incoming profiles (rhotor)
% rho: RAPTOR rho grid (rhotor)
% rhogauss: RAPTOR gauss grid (rhotor)
%
% geom: structure containing splines and index vectors

ord = 3; % spline order

nkts = ceil(numel(rhoin)/2); % number of knots
assert(nkts>=2,'too few points for rhoin')

%%% denser packing near edge
%x=linspace(0,1,nkts);
%knots(x<0.5) = sqrt(2)*x(x<0.5);%
%knots(x>=0.5) = sqrt(1-2*(1-x(x>=0.5)).^2);

% no denser packing near the edge (in rhotor, this is not appropriate. Only for rhopsi)
knots = linspace(0,1,nkts);
%%

% matrices to map to RAPTOR rtorN grid
[S,Sp,Spp] = splineval(rho,knots,ord);
%
% define transformation matrices such that
% * first coefficient is value 0
% * second coefficient is 1st derivative at 0
% * third coefficient is 2nd derivative at 0
%  (similar (mirrored) for end conditions)
C0 = [S(:,1)';Sp(:,1)';Spp(:,1)'];
Tf = eye(size(S,1));
Tf(:,1:size(C0,1)) = pinv(full(C0));
%C1 = [Spp(:,end)';Sp(:,end)';S(:,end)'];
%Tf(:,end-size(C1,1)+1:end) = pinv(full(C1));

% matrices to
geom.LamG = S'*Tf; geom.LamGp= Sp'*Tf; geom.LamGpp = Spp'*Tf;

% same to RAPTOR rhogauss
[S,Sp] = splineval(rhogauss,knots,ord);
geom.LamGgauss = S'*Tf; geom.LamGpgauss = Sp'*Tf;

% same to rhoin: rhogrid of profiles (to map from equilibrium code)
[S] = splineval(rhoin,knots,ord);
geom.LamGin = S'*Tf;

% indices
geom.nsp = size(S,1);

geom.ind_g1    = 1:geom.nsp;
geom.ind_g2    = 1*geom.nsp + (1:geom.nsp);
geom.ind_g3    = 2*geom.nsp + (1:geom.nsp);
geom.ind_g23or = 3*geom.nsp + (1:geom.nsp);
geom.ind_F     = 4*geom.nsp + (1:geom.nsp);
geom.ind_Vp    = 5*geom.nsp + (1:geom.nsp);
geom.ind_epsilon   = 6*geom.nsp + (1:geom.nsp);
geom.ind_Rgeom = 7*geom.nsp + (1:geom.nsp);
geom.ind_delta = 8*geom.nsp + (1:geom.nsp);
geom.ind_Phib  = 9*geom.nsp + 1;

% for debugging
geom.rhoin = rhoin;
geom.rho = rho;
geom.rhogauss = rhogauss; 
return