function g = geom2RAPTOR(g1,g2,g3,g23or,F,Vp,epsilon,Rgeom,delta,Phib,geom)
% maps geometric profiles from grid of source code to RAPTOR g vector
% uses spline matrics in geom structure, which is calculated by prep_geom.m

%#codegen

% init
nsp = size(geom.LamGin,2);
nbt = size(g1,2);

coeff_g1 = zeros(nsp,1);
coeff_g2 = zeros(nsp,1);
coeff_g3 = zeros(nsp,1);
coeff_g23or = zeros(nsp,1);
coeff_Vp = zeros(nsp,1);
ng = nsp*9 + 1; % total g coefficients to store

g = zeros(ng,nbt);
for ii=1:nbt
    % impose boundary conditions directly via spline coefficients
    ifix = false(nsp,1); ifix([1,2]) = true; % value and first derivative 0 at 0;
    coeff_g1(~ifix,:) = geom.LamGin(:,~ifix)\g1(:,ii);
    coeff_g1(ifix) = 0; % g1(0),g1'(0) = 0;

    ifix = false(nsp,1); ifix([1,2]) = true; % value and first derivative 0 at 0;
    coeff_g2(~ifix,:) = geom.LamGin(:,~ifix)\g2(:,ii);
    coeff_g2(ifix) = 0; % g2(0),g2'(0) = 0;

    ifix = false(nsp,1);  % all free
    coeff_g3(~ifix,:) = geom.LamGin(:,~ifix)\g3(:,ii);
    coeff_g3(ifix) = 0; % no constraints

    ifix = false(nsp,1); ifix([1,3]) = true; % value and second derivative 0 at 0;
    coeff_g23or(~ifix,:) = geom.LamGin(:,~ifix)\(g23or(:,ii));
    coeff_g23or(ifix) = 0; %

    % F = RBphi: free shape, no constraints..
    coeff_F = geom.LamGin\F(:,ii);

    % dV/drhoN
    ifix = false(nsp,1); ifix([1,3]) = true; % value 0 at 0;
    coeff_Vp(~ifix,:) = geom.LamGin(:,~ifix)\(Vp(:,ii));
    coeff_Vp(ifix,:)  = 0; % Vp(0) = 0

    coeff_epsilon = geom.LamGin\epsilon(:,ii);
    coeff_Rgeom = geom.LamGin\Rgeom(:,ii);
    coeff_delta = geom.LamGin\delta(:,ii);

    %% assign to g
    g(geom.ind_g1,ii)    = coeff_g1;
    g(geom.ind_g2,ii)    = coeff_g2;
    g(geom.ind_g3,ii)    = coeff_g3;
    g(geom.ind_g23or,ii) = coeff_g23or;
    g(geom.ind_F,ii)     = coeff_F;
    g(geom.ind_Vp,ii)    = coeff_Vp;
    g(geom.ind_epsilon,ii)   = coeff_epsilon;
    g(geom.ind_Rgeom,ii) = coeff_Rgeom;
    g(geom.ind_delta,ii) = coeff_delta;
    g(geom.ind_Phib,ii)  = Phib(ii);
end

%
%% check
% check = false;
% if check
%     rhoin = geom.rhoin;
%     rhogauss = geom.rhogauss;
%     rho = geom.rho;
%
%     clf;
%     subplot(321);
%     plot(rhoin,g1,'bx',rho,geom.LamG*coeff_g1,'rs',rhogauss,geom.LamGgauss*coeff_g1,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_g1,'c')%,rho,geom.LamGpp*coeff_g1,'g');% derivatives
%     title('g1')
%
%     subplot(322);
%     plot(rhoin,g2,'bx',rho,geom.LamG*coeff_g2,'rs',rhogauss,geom.LamGgauss*coeff_g2,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_g2,'c')%,rho,geom.LamGpp*coeff_g2,'g');% derivatives
%     title('g2')
%
%     subplot(323);
%     plot(rhoin,g3,'bx',rho,geom.LamG*coeff_g3,'rs',rhogauss,geom.LamGgauss*coeff_g3,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_g3,'c')%,rho,geom.LamGpp*coeff_g3,'g'); % derivatives
%     title('g3')
%
%     subplot(324);
%     plot(rhoin,g23or,'bx',rho,geom.LamG*coeff_g23or,'rs',rhogauss,geom.LamGgauss*coeff_g23or,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_g23or,'c')%,rho,geom.LamGpp*coeff_g23or,'g'); % derivatives
%     title('g2g3/r')
%
%     subplot(325);
%     plot(rhoin,F,'bx',rho,geom.LamG*coeff_F,'rs',rhogauss,geom.LamGgauss*coeff_F,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_F,'c')%,rho,geom.LamGpp*coeff_F,'g'); % derivatives
%     title('F')
%
%     subplot(326);
%     plot(rhoin,Vp,'bx',rho,geom.LamG*coeff_Vp,'rs',rhogauss,geom.LamGgauss*coeff_Vp,'k^');
%     hold on; plot(rho,geom.LamGp*coeff_Vp,'c')%,rho,geom.LamGpp*coeff_Vp,'g'); % derivatives
%     title('Vp')
% end
