function [X_grid, Y_grid, Z_grid, alpha] = synth_grid_ellipsoid(gh, theta_grid_model, phi_grid_model, varargin)
% [X_grid, Y_grid, Z_grid] = synth_grid_ellipsoid(gh, theta_grid_model, phi_grid_model)
% [X_grid, Y_grid, Z_grid] = synth_grid_ellipsoid(gh, theta_grid_model, phi_grid_model, height)

a = 6371.2; 
rad= pi/180;

theta_grid_model = flipud(theta_grid_model(:));
phi_grid_model = phi_grid_model(:);

if nargin > 3
    h = varargin{1};
else
    h = 0; % altitude above ellipsoid
end

% synthesize model vector values on 3D equidistant grid in spherical coordinates
% create extended fine grid values
% theta_grid_model = [0:0.5:180]'; phi_grid_model = [-180:0.5:180]';
% theta_grid_model(1) = 1e-8; theta_grid_model(end) = 180-1e-8;  % to avoid singularity at poles
r_equator = 6378.14;
r_pole = 6356.75;
r_grid_model = linspace(h+r_pole-.2, h+r_equator+.2, 5)';

% model values on spherical grid, for radii between r_pole and r_equator
[B_r_grid, B_theta_grid, B_phi_grid] = synth_model(gh, r_grid_model/a, theta_grid_model, phi_grid_model);

% augment B_grid

% Step = .25;
Step = abs(mean(diff(theta_grid_model)));
% % equal geodetic latitude grid
alpha = [90:-Step:-90]';
% phi   = [-180+Step/2:Step:180-Step/2]; % longitude
% [alpha_grid, phi_grid] = ndgrid(alpha, phi);

[alpha_grid, phi_grid] = ndgrid(alpha, phi_grid_model);

% convert these geodetic coordinates (geodetic grid) to geocentric coordinates (grid)
[r_grid, theta_grid] = geod2geoc(alpha_grid*rad, h);
theta_grid = theta_grid/rad;

% 3D linear interpolation to spherical components on geodetic grid
B_r     = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_r_grid, phi_grid, theta_grid, r_grid, '*linear');
B_theta = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_theta_grid, phi_grid, theta_grid, r_grid, '*linear');
Y_grid  = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_phi_grid, phi_grid, theta_grid, r_grid, '*linear');
% B_r     = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_r_grid, phi_grid, 90-alpha_grid, r_grid, '*linear');
% B_theta = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_theta_grid, phi_grid, 90-alpha_grid, r_grid, '*linear');
% Y_grid  = interp3(phi_grid_model, theta_grid_model, r_grid_model, B_phi_grid, phi_grid, 90-alpha_grid, r_grid, '*linear');

% transform from spherical to geodetic components
[~, ~, X_grid, Z_grid] = geoc2geod(r_grid, theta_grid*rad, B_r, B_theta);
X_grid = flipud(X_grid);
Y_grid = flipud(Y_grid);
Z_grid = flipud(Z_grid);

