// mesh size description

cl_1   =  2;
cl_2   =  1;

// boundary points that forms Rhizotron
Point(1) = {-5, -10, 0, cl_1};
Point(2) = {5, -10, 0, cl_1};
Point(3) = {5, 0, 0, cl_1};
Point(4) = {-5, 0, 0, cl_1};


// lines that connect boundary
Line(1) = {1, 2};
Line(2) = {3, 2};
Line(3) = {4, 3};
Line(4) = {1, 4};
Line Loop(4) = {1, -2, -3, -4};

// Mesh Parameters
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthMax = 2;


// no. of spline surfaces = 2

X1 ={0,0,-0.184,-0.676,-0.7729,-0.8773,-0.9599,-1.194,-1.315,-1.526,-1.599,-1.679,-1.817,-2.236,-2.705,-2.963,-3.165,-3.365,-3.41,-3.499,-3.567,-3.674,-3.767,-3.864,-4.079,-4.152,-4.176,-4.175,-4.253,-4.349,-4.412,-4.455,-4.512,-4.568,-4.606,-4.641,-4.5433,-4.5083,-4.4747,-4.4217,-4.3628,-4.321,-4.2627,-4.1616,-4.0764,-4.0762,-4.0552,-3.991,-3.7831,-3.6828,-3.5911,-3.4828,-3.404,-3.3124,-3.2721,-3.0845,-2.8959,-2.6369,-2.1671,-1.7449,-1.5997,-1.52,-1.4572,-1.2531,-1.1354,-0.89561,-0.80161,-0.69478,-0.61953,-0.13723,0.088029,0.1,0};
Y1 ={0,-0.2636,-0.3414,-0.6213,-0.7443,-0.8732,-0.9609,-1.139,-1.218,-1.401,-1.487,-1.598,-1.771,-2.177,-2.615,-2.853,-3.031,-3.398,-3.648,-4.004,-4.126,-4.277,-4.422,-4.574,-4.851,-5.108,-5.228,-5.457,-5.685,-5.848,-5.957,-6.081,-6.195,-6.319,-6.439,-6.654,-6.6753,-6.4603,-6.3549,-6.2379,-6.1197,-5.9984,-5.8985,-5.7257,-5.4736,-5.2346,-5.1329,-4.8985,-4.6328,-4.4759,-4.333,-4.18,-4.0352,-3.6696,-3.4349,-3.0904,-2.9272,-2.6882,-2.2495,-1.8403,-1.6589,-1.5483,-1.4736,-1.2965,-1.22,-1.0375,-0.93855,-0.80673,-0.70383,-0.42979,-0.31104,1.837e-17,0};


// Define spline surfaces

LN = 91;
nR = #X1[ ];
p0  =  newp;
p   =  p0;
For i In {0:nR-1}
Point(newp)  =    {X1[i], Y1[i], 0, cl_2};
EndFor
p2  =  newp-1;
Spline(91)   =  {p:p2,p};
Line Loop(91) = {91};
Plane Surface(91) = {91};

X2 ={0,0.6623,1.159,1.524,1.654,1.793,1.786,1.82,1.918,2.473,2.551,2.74,2.804,2.827,2.925,2.902,2.8353,2.6431,2.5577,2.0007,1.9187,1.8859,1.8922,1.7531,1.6208,1.248,0.74707,0.084771,0};
Y2 ={-0.2636,-1.322,-2.116,-3,-4.022,-4.975,-5.054,-5.729,-5.859,-6.69,-6.869,-7.322,-7.666,-7.748,-7.728,-7.646,-7.2917,-6.8301,-6.6369,-5.8028,-5.7128,-5.0504,-4.9623,-4.0085,-2.9749,-2.0703,-1.269,-0.21055,-0.2636};


// Define spline surfaces

LN = 92;
nR = #X2[ ];
p0  =  newp;
p   =  p0;
For i In {0:nR-1}
Point(newp)  =    {X2[i], Y2[i], 0, cl_2};
EndFor
p2  =  newp-1;
Spline(92)   =  {p:p2,p};
Line Loop(92) = {92};
Plane Surface(92) = {92};

SetFactory("OpenCASCADE");
aa() = BooleanUnion{ Surface {91}; Delete; }{ Surface {92}; Delete; };
Physical Surface(1) = {aa()};