• /
• # Planet and moon

on 3 Dec 2023
• 6
• 162
• 0
• 6
• 1981
drawframe(1);
function drawframe(f)
% Procedural planet & moon
% Apologies for the sloppy code and a couple garbage variables. Based
% around Matlon5 code but forgot to add as a remix.
persistent V
if f==1
rng(5,'twister');
v=@vecnorm;
ra=@rand;
bscl = 1e5;
% Usng a random sphere instead of fibonacci sphere like Matlon5
nl=1e5;
r=ra(nl,1)/3+1;
p=randn(3,nl);
p=p./v(p);
[p,k,s]=SM(p',r,5e1);
% Spin points to make a gaseous looking surface...
p2 = (p.'./vecnorm(p.'))';
scl=-7;
p2(:,1:2) = [p2(:,1).*cos(p2(:,3)*scl)-p2(:,2).*sin(p2(:,3)*scl),p2(:,1).*sin(p2(:,3)*scl)+p2(:,2).*cos(p2(:,3)*scl)];
% Re-smooth
[p3,k2,s2]=SM(p2,s,3e0);
p3 = (p3.'./vecnorm(p3.'))';
% Add a storm in it
xp=[1,1,0.5]';
xp = xp/norm(xp);
p2=p3';
for nn = 1:100
pob = p2./v(p2);
for n = 1:size(xp, 2)
if n == 1
xc = XP(xp(:,n),pob);
else
xc = xc + XP(xp(:, n),pob);
end
end
p2 = (p2./v(p2) + xc/1500).*s';
end
p3 = (p2./vecnorm(p2))';
E='EdgeC';
F='FaceC';
O='none';
G='FaceA';
% Plot planet
T=trisurf(k2,p3(:,1)*bscl,p3(:,2)*bscl,p3(:,3)*bscl,s2,E,O); % Plot
material([0,1,0,10]);
% Make it a purple planet
c1=[88,72,154];
c2=[112,89,145];
c3=[140,120,140];
c4=[180,160,170];
c5=[250,240,242];
c=[c1;c2;c3;c4;c5];
cmp = interp1(1:5, c, linspace(1, 5, 256))/255;
colormap(min(cmp.^1.5*1.5,1));
cb = caxis;
% % % Add a moon...
rng default
v=@vecnorm;
nl=1e4;
rd = rand(nl,1)/3+1;
po = randn(3,nl);
po = po./vecnorm(po);
[p,k,s]=SM(po',rd,7);
ms = mean(s);
ns = s - ms;
s = erf(ns*500)/50+ms;
[~,~,sb]=SM(po',randn(nl,1),1);
% Now, take a page from Adam D.' book for the craters
rng(1,'twister');
x = randn(3,16); % Let's add 16 craters
x = x./vecnorm(x);
mfc = @(x,y,z)y.*(erf((x-z)*30)/2+.5); % Using erf as the crater function...
for n = 1:size(x, 2)
d = vecnorm(x(:,n)-p');
s = mfc(d',s-1.14,rand(1)/2)+1.14;
end
s = s + sb/20;
s = (s-mean(s))/10+mean(s);
p2 = (p.'./vecnorm(p.'))'.*s;
E='EdgeC';
F='FaceC';
O='none';
G='FaceA';
% Plot moon
hold on;
T2=trisurf(k,bscl*p2(:,1)/3+4*bscl,bscl*p2(:,2)/3-1*bscl,bscl*p2(:,3)/3,'FaceC','flat','FaceVertexCData',rescale(s,.4,1).*[1,1,1],E,O); % Plot
material([0,1,0,10]);
axis equal off
set(gcf,'color','k');
light;
light('color',[1,1,1]*.4);
caxis(cb);
V=hgtransform('Parent',gca);
set(T2,'parent',V);
camproj p;
camtarget([0,0,.2]*bscl);
end
a=-1.45:.018:-.59;
camva(9);
c=pi/4-sin(0:.131:2*pi)/10;
campos([sin(c(f)),cos(c(f)),0]*1.3e6);
set(V,'Matrix',makehgtform('zrotate',-a(f)));
end
function cp=XP(n,p)
n=n/vecnorm(n)*.9;
d=sqrt(sum((n - p).^2));
cp=cross(p, n.*ones(1,size(p,2)))./d.^2;
end
function [p,k,s]=SM(in,rd,r)
n=size(in,1);
k=convhull(in); % Points on "in" must lie on unit circle
c=@(x)sparse(k(:,x)*[1,1,1],k,1,n,n); % Connectivity
t=c(1)|c(2)|c(3); % Connectivity
f=spdiags(-sum(t,2)+1,0,t*1.)*r; % Weighting
s=((speye(n)+f'*f)\rd); % Solve for s w/regularizer
p=in.*s; % Apply
% S.D.G.
end
Tim