Skip to content

geomag2geog, geog2geomag update #22

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
paulinchin wants to merge 16 commits into
base: main
Choose a base branch
from
Open

geomag2geog, geog2geomag update #22

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
dffcf9c
Update geog2geomag.m
paulinchin Nov 25, 2020
5140de0
Update geomag2geog.m
paulinchin Nov 25, 2020
99b253c
Update geog2geomag.m
paulinchin Nov 25, 2020
ed7cda4
Update geog2UEN.m
paulinchin Nov 25, 2020
c1c6380
Update TestAUnit.m
paulinchin Nov 25, 2020
a73242d
Update UEN2geog.m
paulinchin Nov 25, 2020
50d1cc1
Update geog2geomag.m
paulinchin Nov 25, 2020
35d458e
Update geomag2geog.m
paulinchin Nov 25, 2020
51e4050
Update makegrid_tilteddipole_3D.m
paulinchin Nov 25, 2020
ec5e4df
Update makegrid_cart_3D.m
paulinchin Nov 25, 2020
539d740
Update plot_mapgrid.m
paulinchin Nov 25, 2020
f6ff29b
Update TECplot_map.m
paulinchin Nov 25, 2020
b7e2721
Update geog2geomag.m
paulinchin Nov 25, 2020
11cdbaa
Update geomag2geog.m
paulinchin Nov 25, 2020
497bedc
Update grid_examples.m
paulinchin Nov 25, 2020
fdf5e57
Update magplot_fort_map.m
paulinchin Nov 25, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 10 additions & 3 deletions +gemini3d/+setup/+gridgen/makegrid_cart_3D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,15 @@
function xgf = makegrid_cart_3D(p)
function xgf = makegrid_cart_3D(p,year)
%% make 3D grid and save to disk

if ~exist('year','var')
% If year is not specified, set it 1985 by default
year = 1985;
end

arguments
p (1,1) struct
end

%% create altitude grid
%p.alt_min = 80e3;
%p.alt_max = 1000e3;
Expand Down Expand Up @@ -63,7 +70,7 @@
%DISTANCE EW AND NS (FROM ENU (or UEN in our case - cyclic permuted) COORD. SYSTEM) NEED TO BE CONVERTED TO DIPOLE SPHERICAL AND THEN
%GLAT/GLONG - BASICALLY HERE WE ARE MAPPING THE CARTESIAN GRID ONTO THE
%SURFACE OF A SPHERE THEN CONVERTING TO GEOGRAPHIC.
[thetactr,phictr] = gemini3d.geog2geomag(p.glat, p.glon); %get the magnetic coordinates of the grid center, based on user input
[thetactr,phictr] = gemini3d.geog2geomag(p.glat, p.glon, year); %get the magnetic coordinates of the grid center, based on user input

%% Center of earth distance
r=Re+z;
Expand All @@ -84,7 +91,7 @@
phi=repmat(phi,[lx1,1,lx3]);

%% COMPUTE THE GEOGRAPHIC COORDINATES OF EACH GRID POINT
[glatgrid,glongrid] = gemini3d.geomag2geog(theta,phi);
[glatgrid,glongrid] = gemini3d.geomag2geog(theta,phi,year);

%% COMPUTE ECEF CARTESIAN IN CASE THEY ARE NEEDED
xECEF=r.*sin(theta).*cos(phi);
Expand Down
14 changes: 11 additions & 3 deletions +gemini3d/+setup/+gridgen/makegrid_tilteddipole_3D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function xgf = makegrid_tilteddipole_3D(dtheta,dphi,lpp,lqp,lphip,altmin,glat,glon,gridflag)
function xgf = makegrid_tilteddipole_3D(dtheta,dphi,lpp,lqp,lphip,altmin,glat,glon,gridflag,year)

%NOTE THAT INPUTS DTHETA AND DPHI ARE INTENDED TO REPRESENT THE FULL THETA
%AND PHI EXTENTS OF
Expand All @@ -13,11 +13,18 @@
% properly.
%(5) Terrestrial mag. moment and mass are hard-coded in.
%(6) We assume target latitude is northern hemisphere
%(7)
%
%OTHER NOTES:
% - Fortran code will interpret this grid as including ghost cells, so
% if you want a dimension to be size "n" adjust requested grid size so that
% it is "n+4"

if ~exist('year','var')
% If year is not specified, set it 1985 by default
year = 1985;
end

arguments
dtheta {mustBeNumeric}
dphi {mustBeNumeric}
Expand All @@ -28,6 +35,7 @@
glat {mustBeNumeric}
glon {mustBeNumeric}
gridflag {mustBeInteger}
year {mustBeInteger}
end

%% PAD GRID WITH GHOST CELLS
Expand All @@ -42,7 +50,7 @@


%TD SPHERICAL LOCATION OF REQUESTED CENTER POINT
[thetatd,phid]=gemini3d.geog2geomag(glat,glon);
[thetatd,phid]=gemini3d.geog2geomag(glat,glon,year);


%SETS THE EDGES OF THE GRID
Expand Down Expand Up @@ -368,7 +376,7 @@

%xg.glat=(pi/2-theta)*180/pi; xg.glon=phi*180/pi*ones(lx(1),lx(2));
for iphi=1:lphi
[glats,glons]=gemini3d.geomag2geog(xg.theta(:,:,iphi),xg.phi(1,1,iphi)*ones(lq,lp));
[glats,glons]=gemini3d.geomag2geog(xg.theta(:,:,iphi),xg.phi(1,1,iphi)*ones(lq,lp),year);
% only meant to work for one latitude at at time
xg.glat(:,:,iphi)=glats;
xg.glon(:,:,iphi)=glons;
Expand Down
12 changes: 8 additions & 4 deletions +gemini3d/+setup/plot_mapgrid.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,9 @@
function ha = plot_mapgrid(xg,flagsource,neuinfo)
function ha = plot_mapgrid(xg,flagsource,neuinfo,year)

if ~exist('year','var')
% If year is not specified, set it 1985 by default
year = 1985;
end

arguments
xg (1,1) struct
Expand All @@ -12,7 +17,6 @@
dmlon=max(mlon(:))-min(mlon(:));
dmlat=max(mlat(:))-min(mlat(:));


%% ORGANIZE INPUT STRUCTURE
if flagsource ~= 0
neugridtype=neuinfo.neugridtype;
Expand All @@ -29,7 +33,7 @@
rhomax=neuinfo.rhomax;
end %if

[sourcetheta,sourcephi]=geog2geomag(sourcelat,sourcelong);
[sourcetheta,sourcephi]=geog2geomag(sourcelat,sourcelong,year);
sourcemlat=90-sourcetheta*180/pi;
sourcemlon=sourcephi*180/pi;

Expand Down Expand Up @@ -421,7 +425,7 @@
hold on;
%ax=axis;
%load coastlines;
[thetacoast,phicoast]=geog2geomag(coastlat,coastlon);
[thetacoast,phicoast]=geog2geomag(coastlat,coastlon,year);
mlatcoast=90-thetacoast*180/pi;
mloncoast=phicoast*180/pi;

Expand Down
8 changes: 4 additions & 4 deletions +gemini3d/+tests/TestAUnit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,16 +119,16 @@ function test_max_mpi(tc)

function test_coord(tc)

[lat, lon] = gemini3d.geomag2geog(pi/2, pi/2);
[lat, lon] = gemini3d.geomag2geog(pi/2, pi/2, 2020);
tc.verifyEqual([lat, lon], [0, 19], 'AbsTol', 1e-6, 'RelTol', 0.001)

[theta, phi] = gemini3d.geog2geomag(0,0);
[theta, phi] = gemini3d.geog2geomag(0,0, 2020);
tc.verifyEqual([theta, phi], [1.50863496978059, 1.24485046147953], 'AbsTol', 1e-6, 'RelTol', 0.001)

[alt, lon, lat] = gemini3d.UEN2geog(0,0,0, pi/2, pi/2);
[alt, lon, lat] = gemini3d.UEN2geog(0,0,0, pi/2, pi/2, 2020);
tc.verifyEqual([alt, lat, lon], [0, 0, 19], 'AbsTol', 1e-6, 'RelTol', 0.001)

[z,x,y] = gemini3d.geog2UEN(0, 0, 0, pi/2, pi/2);
[z,x,y] = gemini3d.geog2UEN(0, 0, 0, pi/2, pi/2, 2020);
tc.verifyEqual([z,x,y], [0, -2076275.16205889, 395967.844181141], 'AbsTol', 1e-6, 'RelTol', 0.001)

end
Expand Down
5 changes: 3 additions & 2 deletions +gemini3d/UEN2geog.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
function [alt,glon,glat] = UEN2geog(z,x,y,thetactr,phictr)
function [alt,glon,glat] = UEN2geog(z,x,y,thetactr,phictr,year)
arguments
z {mustBeNumeric}
x {mustBeNumeric}
y {mustBeNumeric}
thetactr {mustBeNumeric}
phictr {mustBeNumeric}
year {mustBeNumeric}
end

% converts magnetic up, north, east coordinates into geographic
Expand Down Expand Up @@ -43,6 +44,6 @@


%Now convert the magnetic to geographic using our simple transformation
[glat,glon] = gemini3d.geomag2geog(theta,phi);
[glat,glon] = gemini3d.geomag2geog(theta,phi,year);

end
4 changes: 2 additions & 2 deletions +gemini3d/geog2UEN.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [z,x,y] = geog2UEN(alt,glon,glat,thetactr,phictr)
function [z,x,y] = geog2UEN(alt,glon,glat,thetactr,phictr,year)
arguments
alt {mustBeNumeric}
glon {mustBeNumeric}
Expand All @@ -16,7 +16,7 @@
z=alt;

%Convert to geomganetic coordinates
[theta,phi] = gemini3d.geog2geomag(glat,glon);
[theta,phi] = gemini3d.geog2geomag(glat,glon,year);


%Convert to northward distance in meters
Expand Down
40 changes: 37 additions & 3 deletions +gemini3d/geog2geomag.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,45 @@
function [thetat,phit] = geog2geomag(lat,lon)
function [thetat,phit] = geog2geomag(lat,lon,year)

if ~exist('year','var')
% If year is not specified, set it 1985 by default
year = 1985;
end


arguments
lat {mustBeNumeric}
lon {mustBeNumeric}
year {mustBeNumeric}
end

thetan = deg2rad(11);
phin = deg2rad(289);
% Schmidt semi-normalised spherical harmonic coefficients are from http://wdc.kugi.kyoto-u.ac.jp/igrf/coef/igrf13coeffs.html
% Calculation of thetan and phin are based on Millward et al., 1996
% For year < 1985, use the same coefficients as for 1985
case ismember(year, 1920:1987)
thetan=deg2rad(1.1000e+01);
phin=deg2rad(2.8900e+02);
case ismember(year, 1988:1992)
thetan=deg2rad(1.0862e+01);
phin=deg2rad(2.8887e+02);
case ismember(year, 1993:1997)
thetan=deg2rad(1.0677e+01);
phin=deg2rad(2.8858e+02);
case ismember(year, 1998:2002)
thetan=deg2rad(1.0457e+01);
phin=deg2rad(2.8843e+02);
case ismember(year, 2003:2007)
thetan=deg2rad(1.0252e+01);
phin=deg2rad(2.8819e+02);
case ismember(year, 2008:2012)
thetan=deg2rad(9.9840e+00);
phin=deg2rad(2.8779e+02);
case ismember(year, 2013:2017)
thetan=deg2rad(9.6869e+00);
phin=deg2rad(2.8739e+02);
case ismember(year, 2018:2025)
thetan=deg2rad(9.4105e+00);
phin=deg2rad(2.8732e+02);
end

%enforce [0,360] longitude
loncorrected = mod(lon, 360);
Expand Down
39 changes: 36 additions & 3 deletions +gemini3d/geomag2geog.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,44 @@
function [lat,lon] = geomag2geog(thetat,phit)
function [lat,lon] = geomag2geog(thetat,phit,year)

if ~exist('year','var')
% If year is not specified, set it 1985 by default
year = 1985;
end

arguments
thetat {mustBeNumeric}
phit {mustBeNumeric}
year {mustBeNumeric}
end

thetan = deg2rad(11);
phin = deg2rad(289);
% Schmidt semi-normalised spherical harmonic coefficients are from http://wdc.kugi.kyoto-u.ac.jp/igrf/coef/igrf13coeffs.html
% Calculation of thetan and phin are based on Millward et al., 1996
% For year < 1985, use the same coefficients as for 1985
case ismember(year, 1920:1987)
thetan=deg2rad(1.1000e+01);
phin=deg2rad(2.8900e+02);
case ismember(year, 1988:1992)
thetan=deg2rad(1.0862e+01);
phin=deg2rad(2.8887e+02);
case ismember(year, 1993:1997)
thetan=deg2rad(1.0677e+01);
phin=deg2rad(2.8858e+02);
case ismember(year, 1998:2002)
thetan=deg2rad(1.0457e+01);
phin=deg2rad(2.8843e+02);
case ismember(year, 2003:2007)
thetan=deg2rad(1.0252e+01);
phin=deg2rad(2.8819e+02);
case ismember(year, 2008:2012)
thetan=deg2rad(9.9840e+00);
phin=deg2rad(2.8779e+02);
case ismember(year, 2013:2017)
thetan=deg2rad(9.6869e+00);
phin=deg2rad(2.8739e+02);
case ismember(year, 2018:2025)
thetan=deg2rad(9.4105e+00);
phin=deg2rad(2.8732e+02);
end

%enforce phit = [0,2pi]
phit = mod(phit, 2*pi);
Expand Down
2 changes: 1 addition & 1 deletion Examples/TECplot_map.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ function TECplot_map(direc)

%ADD A MAP OF COASTLINES
load("coastlines", "coastlat", "coastlon")
[thetacoast,phicoast] = gemini3d.geog2geomag(coastlat,coastlon);
[thetacoast,phicoast] = gemini3d.geog2geomag(coastlat,coastlon,year);
mlatcoast=90-thetacoast*180/pi;
mloncoast=phicoast*180/pi;

Expand Down
6 changes: 3 additions & 3 deletions Examples/grid_examples.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ function grid_examples()
gridflag=1;
flagsource=0;
iscurv=true;

year = 2020; % Adjust magnetic coordinates in accordance to the year required

%% GEOGRAPHIC COORDINATES OF NEUTRAL SOURCE (OR GRID CENTER)
neuinfo.sourcelat=[];
Expand Down Expand Up @@ -84,15 +84,15 @@ function grid_examples()


%% FOR USERS INFO CONVERT SOURCE LOCATION TO GEOMAG
[sourcetheta,sourcephi]=gemini3d.geog2geomag(neuinfo.sourcelat,neuinfo.sourcelong);
[sourcetheta,sourcephi]=gemini3d.geog2geomag(neuinfo.sourcelat,neuinfo.sourcelong,year);
sourcemlat=90-sourcetheta*180/pi;
sourcemlon=sourcephi*180/pi;


%% RUN THE GRID GENERATION CODE
if ~exist('xg', 'var')
if iscurv
xg = gemini3d.setup.gridgen.makegrid_tilteddipole_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
xg = gemini3d.setup.gridgen.makegrid_tilteddipole_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag,year);
% xg=makegrid_tilteddipole_nonuniform_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
% xg=makegrid_tilteddipole_nonuniform_oneside_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
else
Expand Down
2 changes: 1 addition & 1 deletion Examples/magplot_fort_map.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@
%ADD A MAP OF COASTLINES
% if (license('test','Map_Toolbox'))
load coastlines;
[thetacoast,phicoast]= gemini3d.geog2geomag(coastlat,coastlon);
[thetacoast,phicoast]= gemini3d.geog2geomag(coastlat,coastlon,year);
mlatcoast=90-thetacoast*180/pi;
mloncoast=phicoast*180/pi;

Expand Down