File: createTDesignProjectedAndTriangulated.m

package info (click to toggle)
iem-plugin-suite 1.15.0-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 6,080 kB
  • sloc: cpp: 58,973; python: 269; sh: 79; makefile: 41
file content (122 lines) | stat: -rw-r--r-- 2,656 bytes parent folder | download | duplicates (3)
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
clear

% hammer aitov projection
hap = @(azi, ele) [-cos(ele).*sin(azi/2) sin(ele)] ./ (sqrt(1+cos(ele).*cos(azi/2)));

%inverse
z = @(x,y) sqrt(1-(x*2*sqrt(2)/4).^2-(y*sqrt(2)/2).^2);
ihap = @(x,y,z) [2*atan(sqrt(2)*2*-x.*z/2./(2*z.^2-1)) asin(sqrt(2)*y.*z)];
% inverse test
% azi = 0;
% ele = 0;
% 
% xy = hap(azi,ele)
% x = xy(1);
% y = xy(2);
% 
% hz = z(x,y)
% ihap(x,y,z(x,y))

%
nx = 28;
ny = 21;
[X, Y] = meshgrid(linspace(-1,1,nx), linspace(-1,1,ny));
X(1:2:end,:) = X(1:2:end,:) - 1 / nx;

X = X(:);
Y = Y(:);
figure, scatter(X,Y)


aziele = ihap(X,Y,z(X,Y));
azi = aziele(:,1);
ele = aziele(:,2);

rePro = hap(azi,ele);
xre = rePro(:,1);
yre = rePro(:,2);

dmax = 0.00001;
for ii=1:size(X,1)
    r(ii) = abs(X(ii) - xre(ii)) < dmax || abs(Y(ii) - yre(ii)) < dmax;
end

X(~r) = [];
Y(~r) = [];

hold all
scatter(X,Y)

aziele = ihap(X,Y,z(X,Y));
azi = aziele(:,1);
ele = aziele(:,2);

rePro = hap(azi,ele);
xre = rePro(:,1);
yre = rePro(:,2);





tri = delaunay(xre,yre)
triplot (tri, xre, yre, 'Color', [0.2 0.2 0.2]);

scatter (xre, yre, 300,  '.r');
xlim ([-1 1])


[x,y,z] = sph2cart(azi,ele,1);

vertices = [xre yre]';
hammerAitovSampleVertices = vertices(:);
hammerAitovSampleX = x;
hammerAitovSampleY = y;
hammerAitovSampleZ = z;

%%
str = '';
str = [str mat2header('hammerAitovSampleX',hammerAitovSampleX,false)];
str = [str mat2header('hammerAitovSampleY',hammerAitovSampleY,false)];
str = [str mat2header('hammerAitovSampleZ',hammerAitovSampleZ,false)];
str = [str mat2header('hammerAitovSampleVertices',hammerAitovSampleVertices,false)];

indices = (tri-1)';
hammerAitovSampleIndices = indices(:);
str = [str mat2header('hammerAitovSampleIndices',hammerAitovSampleIndices,true)];


fid=fopen('hammerAitovSample.h', 'wt+');
fprintf(fid, '%s', ['#ifndef hammerAitovSample_h ' newline ' #define hammerAitovSample_h ' newline newline '#define nSamplePoints ' num2str(length(hammerAitovSampleX)) newline str newline ' #endif']);
fclose(fid);
%%
X(r) = [];
Y(r) = [];
aziele = ihap(X,Y,z(X,Y));
azi = aziele(:,1);
ele = aziele(:,2);

hold all
scatter(X,Y)


%%
tDesign = getTdesign(21); % tDesign for 10th order (Spherical Harmonic Transform Toolbox)
[tDesign(:,4),tDesign(:,5)] = cart2sph(tDesign(:,1),tDesign(:,2),tDesign(:,3))

tDesignProjected = hap(tDesign(:,4),tDesign(:,5));
scatter(tDesignProjected(:,1),tDesignProjected(:,2));

%%
tri = delaunay(tDesignProjected(:,1),tDesignProjected(:,2))
triplot(tri,tDesignProjected(:,1),tDesignProjected(:,2));

%%
vertices = tDesignProjected';
vertices = vertices(:);
mat2header('tDesignProjected',vertices,0);

%%
indices = (tri-1)';
indices = indices(:);
mat2header('vertexIndex',indices,1);