-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathplm2cube.m
123 lines (119 loc) · 3.65 KB
/
plm2cube.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
function varargout=plm2cube(lmcosi,lfin,alfa,bita,gama)
% v=PLM2CUBE(lmcosi,lfin,alfa,bita,gama)
%
% Transforms an array of real spherical harmonic coefficients into the
% faces of a cubed sphere.
%
% INPUT:
%
% lmcosi The usual spherical harmonic array, as for PLM2XYZ
% lfin Number of subdivisions [default: 6]
% alfa First Euler angle of wholesale tilt of all tiles [defaulted]
% beta Second Euler angle of wholesale tilt of all tiles [defaulted]
% gama Third Euler angle of wholesale tilt of all tiles [defaulted]
%
% OUTPUT:
%
% vp The NxNx6 values to plot on the cubed sphere
%
% EXAMPLE:
%
% plm2cube('demo1','3D',1) % With topography and the continents in 3D
% plm2cube('demo1','2D',1) % With topography and the continents in 2D
% plm2cube('demo2','3D',1) % With Ritsema's model and the continents in 3D
% plm2cube('demo2','2D',1) % With Ritsema's model and the continents in 2D
%
% Last modified by fjsimons-at-alum.mit.edu, 2/23/2011
if ~isstr(lmcosi)
% Default values
defval('lfin',7)
defval('alfa',[]);
defval('bita',[]);
defval('gama',[]);
% Construct the cubed sphere coordinates
[x,y,z]=cube2sphere(lfin,alfa,bita,gama);
% If it takes a while, could expand the whole thing and interpolate to
% the grid
%h=waitbar(0);
% Loop over the faces
for in=1:6
% waitbar(in/6,h,sprintf('Starting face %i/%i',in,6))
% Transform these to longitude and latitude
[phi,piminth,r]=cart2sph(x(:,:,in),y(:,:,in),z(:,:,in));
lon=phi(:)*180/pi; lat=piminth(:)*180/pi;
% Expand the spherical harmonics onto these grid points
% Since I changed my mind on the parity, try both
try
v(:,:,in)=flipud(reshape(plm2xyz(lmcosi,lat,lon),...
[2^lfin+1 2^lfin+1]));
catch
v(:,:,in)=flipud(reshape(plm2xyz(lmcosi,lat,lon),...
[2^lfin 2^lfin]));
end
end
% close(h)
varns={v};
varargout=varns(1:nargout);
elseif strcmp(lmcosi,'demo1')
% This takes the place of 'tipe'
defval('lfin','3D')
% This takes the place of 'kind'
defval('alfa',1)
% Load Earth's topography
lmcosi=fralmanac('GTM3AR','SHM');
% Restrict to a manageable size
L=48;
lmcosi=lmcosi(1:addmup(L),:);
% Perform the transform with standard inputs
v=plm2cube(lmcosi,5);
clf
% Plot using "lfin" as an option
[p,pg]=plotoncube(v,lfin,alfa);
if strcmp(lfin,'3D')
hold on; plotcont([],[],3)
% Adjust the color scale to the data limits
caxis(minmax(v))
elseif strcmp(lfin,'2D')
% Load and convert continents
hold on
[a,b,XYZ]=plotcont([],[],3);
delete(b)
[phi,piminth,r]=cart2sph(XYZ(:,1),XYZ(:,2),XYZ(:,3));
lonc=phi*180/pi; latc=piminth*180/pi;
[xic,etac]=sphere2cube(lonc,latc);
hold on
[pc,pgc]=plotonchunk(xic,etac);
hold off
end
elseif strcmp(lmcosi,'demo2')
% This takes the place of 'tipe'
defval('lfin','3D')
% This takes the place of 'kind'
defval('alfa',1)
% Load Ritsema's model at near 120 km depth
%lmcosi=interpJRmodel(112.875,'xpcube');
% Load Ritsema's model at near 410 km depth
lmcosi=interpJRmodel(406.350,'xpcube');
% Perform the transform with standard inputs
%v=plm2cube(lmcosi);
v=plm2cube(lmcosi,4);
clf
% Plot using "lfin" as an option
[p,pg]=plotoncube(v,lfin,alfa);
if strcmp(lfin,'3D')
hold on; plotcont([],[],3)
% Adjust the color scale to the data limits
caxis(minmax(v))
elseif strcmp(lfin,'2D')
% Load and convert continents
hold on
[a,b,XYZ]=plotcont([],[],3);
delete(b)
[phi,piminth,r]=cart2sph(XYZ(:,1),XYZ(:,2),XYZ(:,3));
lonc=phi*180/pi; latc=piminth*180/pi;
[xic,etac]=sphere2cube(lonc,latc);
hold on
[pc,pgc]=plotonchunk(xic,etac);
hold off
end
end