SetFactory("Built-in");
radius_bottom = 0.2; 
radius_top = 0.3; 
height = 0.7;
ref0 = 0.2;

// Common center for bottom
Point(1) = {0, 0, 0, 1.};
// Base circles bottom 1 in and out
Point(2) = { radius_bottom, 0, 0, ref0};
Point(3) = {-radius_bottom, 0, 0, ref0};
// Circles
Circle(1) = {2,1,3};
Circle(2) = {3,1,2};
Line Loop(5) = {1,2};
// Surface of base
Plane Surface(7) = {5};

// Common center for top
Point(6) = {0, 0, height, ref0};
// Base circles bottom 2 in and out
Point(7)  = { radius_top,  0, height, ref0};
Point(8)  = {-radius_top,  0, height, ref0};
// Circles
Circle(8)  = {7,6,8};
Circle(9)  = {8,6,7};
Line Loop(12) = {8,9};
// Surface of base
Plane Surface(14) = {12};

///////// Linking between bottom 1 and 2
Line(15) = {2,7};
Line(16) = {8,3};
///// Ruled surfaces between bottom 1 and 2
Line Loop(19) = {2,15,-9,16};
Ruled Surface(20) = {19};
Line Loop(21) = {15,8,16,-1};
Ruled Surface(22) = {21};


//////////// Circle
ref1 = 0.2;
radius_circle = 0.05;

z_circle = 0.54;
factor = height / Sqrt(height*height+(radius_top-radius_bottom)*(radius_top-radius_bottom));
z_circlet = z_circle + radius_circle*factor;
z_circleb = z_circle - radius_circle*factor;

radius_at_z  = radius_bottom+z_circle *(radius_top-radius_bottom)/height;
radius_at_zt = radius_bottom+z_circlet*(radius_top-radius_bottom)/height;
radius_at_zb = radius_bottom+z_circleb*(radius_top-radius_bottom)/height;

x_circle = radius_at_z*Cos(Pi/4);
y_circle = radius_at_z*Cos(Pi/4);

x_circlet = radius_at_zt*Cos(Pi/4);
y_circlet = radius_at_zt*Cos(Pi/4);

x_circleb = radius_at_zb*Cos(Pi/4);
y_circleb = radius_at_zb*Cos(Pi/4);

theta1 = radius_circle/radius_at_z;
x_circler = radius_at_z*Cos(Pi/4+theta1);
y_circler = radius_at_z*Sin(Pi/4+theta1);

x_circlel = radius_at_z*Cos(Pi/4-theta1);
y_circlel = radius_at_z*Sin(Pi/4-theta1);


Point(23) = {x_circle,  y_circle,  z_circle,  ref1};
Point(24) = {x_circlet, y_circlet, z_circlet, ref1};
Point(25) = {x_circler, y_circler, z_circle,  ref1};
Point(26) = {x_circleb, y_circleb, z_circleb, ref1};
Point(27) = {x_circlel, y_circlel, z_circle,  ref1};
Point{23} In Surface{22};
Point{24} In Surface{22};
Point{25} In Surface{22};
Point{26} In Surface{22};
Point{27} In Surface{22};

Ellipse(28)  = {24, 23, 25, 25};
Ellipse(29)  = {25, 23, 26, 26};
Ellipse(30) = {26, 23, 27, 27};
Ellipse(31) = {27, 23, 24, 24};
Line Loop(32) = {28,29,30,31};

Ruled Surface(33) = {32};


// Definition of volume
Surface Loop(78) = {7,14,20,22,33};
Volume(79) = {78};

// DEFINING PHYSICAL SURFACES
Physical Surface(100) = {7,14,20,22};
Physical Surface("circle") = {33};
