%% Two-dimensional manifold continuation using the atlas_kd atlas algorithm

%% Pillow

% Without computational domain boundaries
prob = coco_prob();
prob = coco_add_func(prob, 'pillow', @pillow, [], 'zero', 'u0', [1; 0; 0]);
prob = coco_set(prob, 'cont', 'NPR', 100, 'PtMX', 1000);
prob = coco_set(prob, 'cont', 'h_max', 1, 'almax', 30, 'MaxRes', 10);
coco(prob, 'pillow1', [], 2); % Atlas_kd is default for dimension >=2

atlas = coco_bd_read('pillow1', 'atlas');
figure(1)
clf
plot_atlas(atlas, 3);
axis tight
axis equal
view(3)

% With computational domain boundaries
prob = coco_add_pars(prob, '', [1 2 3], {'x' 'y' 'z'});
prob = coco_set(prob, 'cont', 'PtMX', 600);
coco(prob, 'pillow2', [], 2, {'x' 'y' 'z'}, ...
  {[], [-0.5 0.5], [-0.35 0.35]});

atlas = coco_bd_read('pillow2', 'atlas');
figure(2)
clf
plot_atlas(atlas, 3);
axis tight
axis equal
view(-60,20)

%% Ellipsoid

% With computational domain boundary, starting in interior
prob = coco_prob();
prob = coco_add_func(prob, 'ellipsoid', @ellipsoid, [], 'zero', ...
  'u0', [2; 0.1; 0.1]);
prob = coco_add_pars(prob, '', [1 2 3], {'x' 'y' 'z'});
prob = coco_set(prob, 'cont', 'PtMX', 1000);
prob = coco_set(prob, 'cont', 'h', 0.1, 'h_max', 0.2, 'MaxRes', 5);
coco(prob, 'ellipsoid1', [], 2, {'y' 'z' 'x'}, {[0 1] [0 1] []});

atlas = coco_bd_read('ellipsoid1', 'atlas');
figure(1)
clf
plot_atlas(atlas, 3);
axis tight
axis equal
view(-60,20)

% With computational domain boundary, starting on corner
prob = coco_prob();
prob = coco_add_func(prob, 'ellipsoid', @ellipsoid, [], 'zero', ...
  'x0', [1.2; 0; 0]);
prob = coco_add_pars(prob, '', [1 2 3], {'x' 'y' 'z'});
prob = coco_set(prob, 'cont', 'PtMX', 1000);
prob = coco_set(prob, 'cont', 'h', .2, 'h_max', .2, 'MaxRes', 5);
coco(prob, 'ellipsoid2', [], 2, {'y' 'z' 'x'}, {[0 1] [0 1] []});

atlas = coco_bd_read('ellipsoid2', 'atlas');
figure(2)
clf
plot_atlas(atlas, 3);
axis tight
axis equal
view(-60,20)

%% Cylinder

% With computational domain boundary imposed by the limits on 'p'
% restricting continuation to a wedge in three-dimensional space.

prob = coco_add_func(coco_prob(), 'cylinder', @cylinder, [], ...
  'zero', 'u0', [1;0;0]+sqrt([0.5;0.55;0]) );
prob = coco_add_func(prob, 'angle', @angle, [], 'active', 'p', ...
  'uidx', (1:3)');
prob = coco_add_pars(prob, '', [1, 2, 3], {'x' 'y' 'z'});
prob = coco_set(prob, 'cont', 'PtMX', 1000);
prob = coco_set(prob, 'cont', 'h', 0.5, 'h_max', 0.5);
prob = coco_set(prob, 'cont', 'almax', 35, 'MaxRes', 20);
coco(prob, 'cylinder', [], 2, {'x' 'z' 'y' 'p'}, {[], [], [], [-0.5 0.5]});

atlas = coco_bd_read('cylinder', 'atlas'); % Extract the atlas from the bd.mat file
figure(1)
clf
plot_atlas(atlas, 3);
axis tight
axis equal
view(30,10)

%% Sphere

prob = coco_prob();
prob = coco_add_func(prob, 'sphere', @sphere, [], 'zero', 'u0', [1; 1; 0]);
prob = coco_add_pars(prob, '', 1:3, {'x', 'y', 'z'});
prob = coco_set(prob, 'cont', 'PtMX', 1000, 'almax', 15);
prob = coco_set(prob, 'cont', 'h', 1, 'h_max', 1, 'MaxRes',10);
prob = coco_add_event(prob, 'UZ', 'z', 0.5);
coco(prob, 'sphere', [], 2, {'x', 'y', 'z'}, {[0.5 1.5]});

[atlas bd] = coco_bd_read('sphere', 'atlas', 'bd');
figure(4)
clf
plot_atlas(atlas, 3);
idx = coco_bd_idxs(bd, 'UZ');
x   = coco_bd_col(bd, {'x' 'y' 'z'});
hold on
plot3(x(1,idx), x(2,idx), x(3,idx), 'b.', 'MarkerSize', 21);
hold off
axis equal
axis tight
view(60,30)
drawnow

%% Torus

prob = coco_prob();
prob = coco_add_func(prob, 'torus', @torus, [], 'zero', 'u0', [1; 0; 0]);
prob = coco_add_pars(prob, '', [1 2 3], {'x' 'y' 'z'});
prob = coco_set(prob, 'cont', 'NPR', 50, 'PtMX', 500);
prob = coco_set(prob, 'cont', 'h_max', 0.5, 'almax', 25, 'MaxRes', 100);
coco(prob, 'torus', [], 2, {'x' 'y' 'z'});

atlas = coco_bd_read('torus', 'atlas');
figure(1)
clf
plot_atlas(atlas, 3);
axis equal
axis tight
view(-20,50)

%% Klein bottle

prob = coco_prob();
prob = coco_set(prob, 'cont', 'NPR', 100, 'PtMX', 20000);
prob = coco_set(prob, 'cont', 'almax', 35, 'MaxRes', 10);
prob = coco_set(prob, 'cont', 'h', 0.1, 'h_max', 10, 'h_min', 1e-5);
prob = coco_add_func(prob, 'klein', @klein, [], 'zero', ...
  'u0', [0; -1-sqrt(2); 0]);
prob = coco_add_pars(prob, '', [1 2 3], {'x' 'y' 'z'});
coco(prob, 'klein', [], 2, {'x' 'y' 'z'});

atlas = coco_bd_read('klein', 'atlas');
figure(1)
clf
plot_atlas(atlas, 3);
axis equal
axis tight
view(-50,4)
