function geom = prep_geom(rho,rhogauss, config)
% 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

%% Choose best option depending on number of grid points etc
if isempty(config.equi.gknots)
    if numel(rho)>=1
        nkts = ceil(numel(rho)/2); % default number of knots, 1/2 of input rho points
        assert(nkts>=2,'too few points for rhoin')
        knots = linspace(0,1,nkts)';
    else
        knots = [0 0.25 0.5 0.8 1]';
    end
else
    knots = config.equi.gknots; 
end

ord = config.equi.gsplineorder;% spline order
tens0 = config.equi.gspline_axis_tens;% tension for fit at 0
tensb = config.equi.gspline_edge_tens;% tension for fit at edge


%knots = [0 0.25 0.5 0.8 1]';
%knots = linspace(0,1,60)';
%knots = (linspace(0,1,20)').^(.45); % For refined grid

%knots = [0 0.2 0.4 0.6 0.78 0.87 0.95 1]'; % Good one for RT

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


% 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));
geom.LamG = S'*Tf; geom.LamGp= Sp'*Tf;


%%% Reduced set of matrices for kappa, epsilon etc geometry
knotsr = [0 0.5 1]'; % reduced set of knots
[Sr,Spr] = splineval(rho,knotsr,ord);
C0r = [Sr(:,1)';Spr(:,1)'];
Tfr = eye(size(Sr,1));
Tfr(:,1:size(C0r,1)) = pinv(full(C0r));

% save matrices to struct
geom.LamGr = Sr'*Tfr; geom.LamGpr= Spr'*Tfr;

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

[Sr,Spr] = splineval(rhogauss,knotsr,ord);
geom.LamGrgauss = Sr'*Tfr; geom.LamGprgauss = Spr'*Tfr;

%% matrix of second spline derivatives used for smoothing in geom2RAPTOR.m
% second derivative of splines at rho points: for smoothing
[~,~,Spp] = splineval(rho,knots,ord);
Lpp = Spp'*Tf;
tensgrid = (tensb-tens0)*rho+tens0; % varying tension from 0 to boundary
tLpp = bsxfun(@times,tensgrid,Lpp)/numel(rho); 
% divide by rho to get normalization independent of number of grid points used

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

geom.ind_H  = 1:geom.nsp;
geom.ind_gg0     = 1*geom.nsp + (1:geom.nsp);
geom.ind_gg1     = 2*geom.nsp + (1:geom.nsp);
geom.ind_gg2     = 3*geom.nsp + (1:geom.nsp);
geom.ind_gg3     = 4*geom.nsp + (1:geom.nsp);
geom.ind_F       = 5*geom.nsp + (1:geom.nsp);
geom.ind_Vp      = 6*geom.nsp + (1:geom.nsp);
geom.ind_epsilon = 7*geom.nsp + (1:geom.nspr);
geom.ind_Rgeom   = 7*geom.nsp + 1*geom.nspr + (1:geom.nspr);
geom.ind_kappa   = 7*geom.nsp + 2*geom.nspr + (1:geom.nspr);
geom.ind_delta   = 7*geom.nsp + 3*geom.nspr + (1:geom.nspr);
geom.ind_Phib    = 7*geom.nsp + 4*geom.nspr+ 1;
geom.ng = geom.ind_Phib; % size

% for debugging
geom.Tf = Tf;
geom.Tfr = Tfr;

geom.tLpp        = tLpp;
geom.knots = knots;
geom.knotsr = knotsr;
geom.ord = ord;
geom.rho = rho;
geom.rhogauss = rhogauss; 
return