function g = get_geom_for_new_kappa(gref,Ipref,RAPTOR_model,kappaopt,Ipopt)
% function g = get_geom_for_new_kappa(gref,Ipref,RAPTOR_model,kappaopt,Ipopt)
% Get geometrical parameters according to the specific kappa (boundary
% elongation) and plasma current Ip (total) values.
% gref: reference plasma geometry
% Ipref: reference total plasma current
% kappaopt: proposed boundary elongation
% Ipopt: proposed plasma current
% g: new plasma geometry

% Reference kappa
kapparef_rho = eval_kappa([],gref,[],RAPTOR_model,false);
kapparef = kapparef_rho(end,:);
g = zeros(size(gref));

% ggpsV ~ kappa^(0.0)*Ip^(1.0)
g(RAPTOR_model.geom.ind_H,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_H,:),(Ipopt./Ipref).^(1.0));
% ggpsV2 ~ kappa^(0.0)*Ip^(2.0)
g(RAPTOR_model.geom.ind_gg0,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_gg0,:),(Ipopt./Ipref).^(2.0));
% gg1 ~ kappa^(1.2)*Ip^(2.0)
g(RAPTOR_model.geom.ind_gg1,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_gg1,:),...
    (kappaopt./kapparef).^(1.2).*(Ipopt./Ipref).^(2.0));
% gg2 ~ kappa^(1.2)*Ip^(2.0)
g(RAPTOR_model.geom.ind_gg2,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_gg2,:),...
    (kappaopt./kapparef).^(1.2).*(Ipopt./Ipref).^(2.0));
% gg3 ~ kappa^(0.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_gg3,:) = gref(RAPTOR_model.geom.ind_gg3,:);
% F ~ kappa^(0.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_F,:) = gref(RAPTOR_model.geom.ind_F,:);
% Vp ~ kappa^(1.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_Vp,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_Vp,:),(kappaopt./kapparef).^(1.0));
% epsilon ~ kappa^(0.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_epsilon,:) = gref(RAPTOR_model.geom.ind_epsilon,:);
% Rgeom ~ kappa^(0.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_Rgeom,:) = gref(RAPTOR_model.geom.ind_Rgeom,:);
% delta ~ kappa^(0.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_delta,:) = gref(RAPTOR_model.geom.ind_delta,:);
% Phib ~ kappa^(1.0)*Ip^(0.0)
g(RAPTOR_model.geom.ind_Phib,:) = bsxfun(@times,gref(RAPTOR_model.geom.ind_Phib,:),(kappaopt./kapparef).^(1.0));
% kappa
% Construct profiles
% Rho dependence
rho_fnc = bsxfun(@rdivide,kapparef_rho,kapparef_rho(end,:));
kappaopt_rho = bsxfun(@times,rho_fnc,kappaopt);
nt = size(kappaopt_rho,2);
for ii=1:nt
    g(RAPTOR_model.geom.ind_kappa,ii) = RAPTOR_model.geom.LamGr\kappaopt_rho(:,ii);
end

return