function g = update_g_from_profiledata(prof, B0, R0, config, model)
% Build the g vector corresponding to the new equilibrium
iR2_ch = prof.av_1_over_rpow2_av.data / R0^2;
iR_ch = prof.av_1_over_r_avprofile.data / R0;
dpsi_ch = prof.av_grad_psi_av.data * B0*R0;
dpsi2_ch = prof.av_grad_psi_pow2_av.data * B0^2*R0^2;
dpsi2oR2_ch = prof.av_bppow2_av.data * B0^2;
dpsidV_ch = 1./prof.dvdpsi.data * B0/R0;

F_ch = prof.t_eq_rbphi.data * B0*R0;

Psi_cheaseSi = prof.psichease_eq_psi_over_2pi.data * R0^2*B0;
s = prof.s_mesh.data;
C2_ch = prof.av_1_over_rpow2_av.data .* prof.int_rdlp_over_grad_psi_eq_int_jdchi.data ./ (R0*B0);
y = F_ch.*C2_ch;
psi_edge = Psi_cheaseSi(end);
psi_axis = Psi_cheaseSi(1);
Phi_ch = cumtrapz(s.*s.*(psi_edge-psi_axis),y);

Rgeom_ch = R0 * prof.rgeomprofile.data;
epsilon_ch = prof.a_over_rgeom.data;
delta_up = prof.delta_upper.data;
delta_b = prof.delta_bottom.data;
delta_ch = (delta_up + delta_b)*0.5;
kappa_ch = prof.elongation.data;

ij = 2:numel(iR_ch);

iR2 = iR2_ch(ij);
iR  = iR_ch(ij);
dpsi = dpsi_ch(ij);
dpsi2 = dpsi2_ch(ij);
dpsi2oR2 = dpsi2oR2_ch(ij);
dpsidV = dpsidV_ch(ij);
F = F_ch(ij);
Rgeom = Rgeom_ch(ij);
epsilon = epsilon_ch(ij);
kappa = kappa_ch(ij);
delta = delta_ch(ij);
Phib = Phi_ch(end);
rhotor = sqrt(Phi_ch(ij)/Phi_ch(end));

rho = model.rgrid.rho;
rhogauss = model.rgrid.rhogauss;
geom = prep_geom(rho, rhogauss, config);
geom.debug = config.debug.checkgeom;

g = geom2RAPTOR(rhotor, iR, iR2, dpsi, dpsi2, dpsi2oR2, dpsidV, F, epsilon, Rgeom, kappa, delta, Phib, geom);
end