%% Iowa grid for AGU 2019
%{
dtheta=16;
dphi=29;
lp=100;
lq=200;
lphi=100;
altmin=80e3;
glat=41.5;   %38.9609;
glon=360-94.088;
gridflag=1;
flagsource=1;
iscurv=true;
%}

% %% Moore OK
% p.dtheta=20;
% p.dphi=27.5;
% p.lp=256;
% p.lq=256;
% p.lphi=210;
% p.altmin=80e3;
% p.glat=39;
% p.glon=262.51;
% p.gridflag=1;
% p.flagsource=1;
% p.iscurv=true;
% 
% %% HWM tests
% p.glat = 42.45;
% p.glon = 143.4;
% p.dtheta = 15;
% p.dphi=1;
% p.altmin = 80e3;
% p.lp = 264;
% p.lq = 512;
% p.lphi = 1;
% p.gridflag = 1;
% p.flagsource=1;
% p.iscurv=true;

% %% HWM tests
% p.glat = 42.45;
% p.glon = 143.4;
% p.dtheta = 15;
% p.dphi=1;
% p.altmin = 80e3;
% p.lp = 264;
% p.lq = 512;
% p.lphi = 1;
% p.gridflag = 1;
% p.flagsource=1;
% p.iscurv=true;

% %% midlat extension to wind tests
% p.glat = 52;
% p.glon = 143.4;
% p.dtheta = 27;
% p.dphi = 356.2887;
% p.altmin = 80e3;
% p.lp = 264;
% p.lq = 384;
% p.lphi = 24;
% p.gridflag = 1;

%% EQuatorial grid
%%p.dtheta=5.75;
%p.dtheta=7.75;
%p.dphi=35;
%p.lp=128;
%p.lq=256;
%p.lphi=48;
%p.altmin=80e3;
%p.glat=10;
%p.glon=360-76.9;     %Jicamarca
%p.gridflag=1;
%p.flagsource=0;

%Perkins
p.dtheta=25;
p.dphi=35;
p.lp=256;
p.lq=192;
p.lphi=32;
p.altmin=80e3;
p.glat=39;
p.glon=262.51;
p.gridflag=0;
p.flagsource=1;
p.iscurv=true;

%% Compute the grid
xg= gemini3d.grid.tilted_dipole(p);


%% DETERMINE GRID RESOLUTION
[dl1,dl2,dl3]= gemscr.setup.gridres(xg);


%% Plot the grid resolution on an appropriate set of physical axes (lat/lon)
ymd=[1981,05,22];
UTsec=4*3600+21*60;
saveplot_fmt={};
[theta,phi]= gemini3d.geog2geomag(38.96,360-94.088);
mlatsrc=90-theta*180/pi;
mlonsrc=phi*180/pi;
sourceloc=[mlatsrc,mlonsrc];

parm=log10(dl1/1e3);
parmlbl='log_{10} x_1 grid spacing (km)';
caxlims=[min(parm(:)),max(parm(:))];
gemini3d.plot.curv3d_long(datetime([ymd,0,0,UTsec]),xg,parm,parmlbl,caxlims,sourceloc);

parm=log10(dl2/1e3);
parmlbl='log_{10} x_2 grid spacing (km)';
caxlims=[min(parm(:)),max(parm(:))];
gemini3d.plot.curv3d_long(datetime([ymd,0,0,UTsec]),xg,parm,parmlbl,caxlims,sourceloc);



%% PLOTS
figure;
imagesc(log10(dl1(:,:,1)/1e3));
colorbar;
title('log_{10} Resolution along field line (min x3)')
datacursormode;

figure;
imagesc(xg.alt(:,:,1)/1e3);
colorbar;
title('altitude (min x3)')
datacursormode;

figure;
imagesc(dl2(:,:,1)/1e3);
colorbar;
title('Resolution perp. to field line (min x3)')
datacursormode;

figure;
imagesc(squeeze(dl3(1,:,:))/1e3);
colorbar;
title('ZONAL Resolution perp. to field line (min x1)')
datacursormode;

figure;
imagesc(squeeze(dl3(floor(end/2),:,:))/1e3);
colorbar;
title('ZONAL Resolution perp. to field line (mid x1)')
datacursormode;


%% FOR A GIVEN GRID ATTEMPT TO DEFINE A CONSTANT STRIDE IN THE X2 DIRECTION (VERY OFTEN DESIRABLE)
lx1=xg.lx(1);
lx2=xg.lx(2);
lx3=xg.lx(3);
h2=xg.h2(3:end-2,3:end-2,3:end-2);
x2=xg.x2(3:end-2);

%ix1=max(floor(lx1/2),1);
ix1=floor(lx1/4);
ix3=max(floor(lx3/2),1);
dl2trial=squeeze(dl2(ix1,:,ix3));                % step size for the original grid
%dl2target=6.75e3;                                % define a desired grid step size in x2
%dl2target=7.5e3;
dl2target=30e3;
dx2target=squeeze(dl2target./h2(ix1,:,ix3));     % what dx2 needs to be to hit target grid size

%A polynomial is being fitted to the target dL as a function of L, this
%can, in turn be used in the grid generation scripts which can handle
%nonuniform stepping.
coeffs=polyfit(x2(:),dx2target(:),2);          % most of the time a quadratic fit gets pretty close
dx2new=polyval(coeffs,x2);
dl2new=dx2new(:).*squeeze(h2(ix1,:,ix3))';

figure;
subplot(121);
plot(x2,dx2new,x2,dx2target);
xlabel('L shell');
ylabel('\Delta L');
legend('dL of nonuniform grid','target dL');
title(sprintf('poly coeffs:  %f %f %f',coeffs))

subplot(122);
plot(x2,dl2new,x2,dl2trial);
xlabel('L shell');
ylabel('grid step size (m)');
legend('dL of nonuniform grid','dL of source (original) grid');
