r=0.5*(R*R+H*H)/H; 
h2=0.00;//value to evaluate point D coordinates
//epsi=1.0e-08;
epsi=0.0;
// make the breast geometry
tga=Tan(a);
b_left = 0;
b_right = 0;
tga_left=Tan(a_left);
tga_right=Tan(a_right);
tgb_left=Tan(b_left);
tgb_right=Tan(b_right);

shift=r-H;

dl=-dl;
dr=-dr;


Point(1) = {0+epsi,0.0+shift,0.0,lc};//point O
Ax=0.0;Az=s;Ay=-Sqrt(r*r-s*s);
Point(2) = {Ax+epsi,Ay+shift,Az,lc};//point A
Ex=0.0;Ez=R;Ey=H-r;
Point(3) = {Ex+epsi,Ey+shift,Ez,lc};// point E


//------------//
Delta_bl=4*Az*Az*tgb_left*tgb_left*tgb_left*tgb_left-4*(1+tgb_left*tgb_left)*(tgb_left*tgb_left*Az*Az+0.0*0.0-R*R);
delta_z=0.5*(2*Az*tgb_left*tgb_left-Sqrt(Delta_bl))/(1+tgb_left*tgb_left);
Dy=-tgb_left*(Az-delta_z);

Delta=4*Az*Az*tga_left*tga_left*tga_left*tga_left-4*(1+tga_left*tga_left)*(tga_left*tga_left*Az*Az+Dy*Dy-R*R);
Dz=0.5*(2*Az*tga_left*tga_left-Sqrt(Delta))/(1+tga_left*tga_left);Dx=tga_left*(Az-Dz);


Dy=-h2l+shift;
delta_M=Tan(-a_left+(Pi*0.5));
delta_A=1+(delta_M*delta_M);
delta_B=2*Az*delta_M;
delta_C=(Az*Az)+(h2l*h2l)-(r*r);
Dx=(-delta_B-Sqrt((delta_B*delta_B) - (4*delta_A*delta_C)))/(2*delta_A);
Dz=Az+delta_M*Dx;


Point(4) = {Dx,Dy,Dz,lc}; //point D

Printf("Ponto 4:: %g,%g,%g -- tan:%g -- r:%g -- delta:%g,%g,%g", Dx, Dy, Dz, delta_M, r, delta_A, delta_B, delta_C);

Delta_br=4*Az*Az*tgb_right*tgb_right*tgb_right*tgb_right-4*(1+tgb_right*tgb_right)*(tgb_right*tgb_right*Az*Az+0.0*0.0-R*R);
delta_z=0.5*(2*Az*tgb_right*tgb_right-Sqrt(Delta_br))/(1+tgb_right*tgb_right);
Dy_r=-tgb_right*(Az-delta_z);

Delta_r=4*Az*Az*tga_right*tga_right*tga_right*tga_right-4*(1+tga_right*tga_right)*(tga_right*tga_right*Az*Az+Dy_r*Dy_r-R*R);
Dz_r=0.5*(2*Az*tga_right*tga_right-Sqrt(Delta_r))/(1+tga_right*tga_right);Dx_r=tga_right*(Az-Dz_r);

Dy_r=-(h2r-shift);
delta_M=Tan(a_right-(Pi*0.5));
delta_A=1+(delta_M*delta_M);
delta_B=2*Az*delta_M;
delta_C=(Az*Az)+(h2r*h2r)-(r*r);
Dx_r=(-delta_B+Sqrt((delta_B*delta_B) - (4*delta_A*delta_C)))/(2*delta_A);
Dz_r=Az+delta_M*Dx_r;

Point(5) = {Dx_r,Dy_r,Dz_r,lc};

Printf("Ponto 5:: %g,%g,%g -- tan:%g -- r:%g -- delta:%g,%g,%g", Dx_r, Dy_r, Dz_r, delta_M, r, delta_A, delta_B, delta_C);

//pontos laterais obtidos com um certo ângulo
Cz=-dr;Cy=Ey;Cx=Sqrt(R*R-(dr*dr));
Point(6) = {Cx+epsi,Cy+shift,Cz,lc};// Point C
Cz_r=-dl;Cy_r=Ey;Cx_r=-Sqrt(R*R-(dl*dl));
Point(7) = {Cx_r+epsi,Cy_r+shift,Cz_r,lc};

Fx=0.0;Fz=0.;Fy=H-r;
Point(8) = {Fx,Fy+shift,Fz,lc};// point F centre of base circle
Bx=Ax;By=Ey;Bz=Az;
Point(9) = {Bx+epsi,By+shift,Bz,lc};// Point B
ux=(Cy-By)*(Dz-Bz)-(Cz-Bz)*(Dy-By);
uy=(Cz-Bz)*(Dx-Bx)-(Cx-Bx)*(Dz-Bz);
uz=(Cx-Bx)*(Dy-By)-(Cy-By)*(Dx-Bx);
lambda=(Bx*ux+By*uy+Bz*uz)/(ux*ux+uy*uy+uz*uz);
Gx=lambda*ux;Gy=lambda*uy;Gz=lambda*uz;

ux=(Cy_r-By)*(Dz_r-Bz)-(Cz_r-Bz)*(Dy_r-By);
uy=(Cz_r-Bz)*(-Dx_r-Bx)-(Cx_r-Bx)*(Dz_r-Bz);
uz=(Cx_r-Bx)*(Dy_r-By)-(Cy_r-By)*(-Dx_r-Bx);
lambda=(Bx*ux+By*uy+Bz*uz)/(ux*ux+uy*uy+uz*uz);
Gx=lambda*ux;Gy=lambda*uy;Gz=lambda*uz;

//Fx=0.0;Fz=0.;Fy=0;
//Point(10) = {Fx,Fy+shift,Fz,lc};// Point B
//Point(11) = {Fx,Fy+shift,Fz,lc};// Point B

//center of arc
//P1 - Cx,Cy+shift,Cz
//P2 - Dx_r,Dy_r,Dz_r
Bx=Ax;By=0;Bz=Az;
p1x=Cx; p1y=Cy+shift; p1z=Cz;
p2x=Dx_r; p2y=Dy_r; p2z=Dz_r; 
n1 = ((p1y-By)*(p2z-Bz)) - ((p2y-By)*(p1z-Bz));
n2 = ((p2x-Bx)*(p1z-Bz)) - ((p1x-Bx)*(p2z-Bz));
n3 = ((p1x-Bx)*(p2y-By)) - ((p2x-Bx)*(p1y-By));
gamma = (p1x*p1x) - (p2x*p2x) + (p1y*p1y) - (p2y*p2y) + (p1z*p1z) - (p2z*p2z);

CCx=0;
CCz = (gamma - ( ( ((n1*Bx) + (n3*Bz))/(n2) ) * ((2*p1y) - (2*p2y)) ) ) / ((2*p1z) - (2*p2z) - ((n3/n2)*((2*p1y)-(2*p2y))) );
CCy = ((n1*Bx) + (n3*Bz) - (n3*CCz)) / n2;

Point(10) = {CCx,CCy,CCz,lc};// Point B

Bx=Ax;By=0;Bz=Az;
p1x=Cx; p1y=Cy+shift; p1z=Cz;
p2x=Dx; p2y=Dy; p2z=Dz; 
n1 = ((p1y-By)*(p2z-Bz)) - ((p2y-By)*(p1z-Bz));
n2 = ((p2x-Bx)*(p1z-Bz)) - ((p1x-Bx)*(p2z-Bz));
n3 = ((p1x-Bx)*(p2y-By)) - ((p2x-Bx)*(p1y-By));
gamma = (p1x*p1x) - (p2x*p2x) + (p1y*p1y) - (p2y*p2y) + (p1z*p1z) - (p2z*p2z);

CCx=0;
CCz = (gamma - ( ( ((n1*Bx) + (n3*Bz))/(n2) ) * ((2*p1y) - (2*p2y)) ) ) / ((2*p1z) - (2*p2z) - ((n3/n2)*((2*p1y)-(2*p2y))) );
CCy = ((n1*Bx) + (n3*Bz) - (n3*CCz)) / n2;


Point(11) = {CCx,CCy,CCz,lc};// Point B

d_val = Sqrt((Dx-Ax)*(Dx-Ax) + (Dy-Ay+shift)*(Dy-Ay+shift) + (Dz-Az)*(Dz-Az)); 
beta = 1/(Sin((d_val*0.5)/R));
t = Cos(beta) * R; //dist

meanX = (Dx + Ax)/2.0; meanY = (Dy+(Ay+shift))/2.0; meanZ = (Dz + Az)/2.0;
meanDist = Sqrt((meanX-Bx)*(meanX-Bx) + (meanY-By)*(meanY-By) + (meanZ-Bz)*(meanZ-Bz)); 
vx = (Bx-meanX)/meanDist;
vy = (By-meanY)/meanDist;
vz = (Bz-meanZ)/meanDist;

//Point(12) = {-(meanX + vx*t), meanY + vy*t, meanZ + vz*t, lc};

h_aux = Sqrt(R*R - (d_val*0.5)*(d_val*0.5));
v1x = -Tan(a_left - (Pi*0.5));
v1y = 0.0;
v1z = 1.0;

ddx = -v1z*(Dy-(Ay+shift));
ddy = -v1x*(Dz-Az) + v1z*(Dx-Ax);
ddz = v1x*(Dy-(Ay+shift));
//ddx = -v1z*(Dy_r-(Ay+shift));
//ddy = -v1x*(Dz_r-Az) + v1z*(Dx_r-Ax);
//ddz = v1x*(Dy_r-(Ay+shift));

P_cX = (meanX+h_aux*ddx);
P_cY = (meanY+h_aux*ddy);
P_cZ = (meanZ+h_aux*ddz);

teste = Sqrt((Ax-P_cX)*(Ax-P_cX) + ((Ay+shift)-P_cY)*((Ay+shift)-P_cY) + (Az-P_cZ)*(Az-P_cZ) );
teste1 = Sqrt((Dx-P_cX)*(Dx-P_cX) + ((Dy)-P_cY)*((Dy)-P_cY) + (Dz-P_cZ)*(Dz-P_cZ) );
teste2 = Tan(Pi*0.5-a_left)*P_cX - P_cZ + Az;
//Printf("teste esq:: %g <-> %g -- maisumteste:%g",teste, teste1, teste2);

//Point(12) = {-(meanX+h_aux*ddx),meanY+h_aux*ddy,meanZ+h_aux*ddz,lc};

lambda = (-Az) / (Tan(-a_left+(Pi*0.5))*Tan(-a_left+(Pi*0.5)) + 1.0);
c_auxX = (lambda*(Tan(-a_left+(Pi*0.5))));
c_auxY = shift;
c_auxZ = -lambda;

Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-Ay+shift)*(c_auxY-Ay+shift) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx_r)*(c_auxX-Dx_r) + (c_auxY-Dy_r)*(c_auxY-Dy_r) + (c_auxZ-Dz_r)*(c_auxZ-Dz_r)));
Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-Ay+shift)*(c_auxY-Ay+shift) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx)*(c_auxX-Dx) + (c_auxY-Dy)*(c_auxY-Dy) + (c_auxZ-Dz)*(c_auxZ-Dz)));


Point(12) = {c_auxX,c_auxY,c_auxZ,lc};

teste = Sqrt((Ax-c_auxX)*(Ax-c_auxX) + ((Ay+shift)-c_auxY)*((Ay+shift)-c_auxY) + (Az-c_auxZ)*(Az-c_auxZ) );
teste1 = Sqrt((Dx-c_auxX)*(Dx-c_auxX) + ((Dy)-c_auxY)*((Dy)-c_auxY) + (Dz-c_auxZ)*(Dz-c_auxZ) );
teste2 = Tan((Pi*0.5)-a_left)*c_auxX - c_auxZ + Az;
Printf("teste esq:: %g <-> %g -- isplane:%f -- lambda:%g",teste, teste1, teste2, lambda);


d_val = Sqrt((Dx_r-Ax)*(Dx_r-Ax) + (Dy_r-Ay+shift)*(Dy_r-Ay+shift) + (Dz_r-Az)*(Dz_r-Az)); 

beta = 1/(Sin((d_val*0.5)/R));
t = Cos(beta) * R; //dist

meanX = (Dx_r + Ax)/2.0; meanY = (Dy_r+(Ay+shift))/2.0; meanZ = (Dz_r + Az)/2.0;
meanDist = Sqrt((meanX-Bx)*(meanX-Bx) + (meanY-By)*(meanY-By) + (meanZ-Bz)*(meanZ-Bz));
vx = (Bx-meanX)/meanDist;
vy = (By-meanY)/meanDist;
vz = (Bz-meanZ)/meanDist;
//Point(13) = {(meanX + vx*t), meanY + vy*t, meanZ + vz*t, lc};


h_aux = Sqrt(R*R - (d_val*0.5)*(d_val*0.5));
v1x = -Tan(a_right - Pi/2.0);
v1y = 0.0;
v1z = 1.0;

ddx = -v1z*(Dy_r-(Ay+shift));
ddy = -v1x*(Dz_r-Az) + v1z*(Dx_r-Ax);
ddz = v1x*(Dy_r-(Ay+shift));

teste = (Dx_r - Ax)*Tan(a_right-(Pi*0.5)) - (Dz_r-Az);

//Printf("teste=%g",teste);

//Point(13) = {(meanX+h_aux*ddx),meanY+h_aux*ddy,meanZ+h_aux*ddz,lc};

P_cX = (meanX+h_aux*ddx);
P_cY = (meanY+h_aux*ddy);
P_cZ = (meanZ+h_aux*ddz);

teste = Sqrt((Ax-P_cX)*(Ax-P_cX) + ((Ay+shift)-P_cY)*((Ay+shift)-P_cY) + (Az-P_cZ)*(Az-P_cZ) );
teste1 = Sqrt((Dx_r-P_cX)*(Dx_r-P_cX) + ((Dy_r)-P_cY)*((Dy_r)-P_cY) + (Dz_r-P_cZ)*(Dz_r-P_cZ) );

//Printf("teste dir:: %g <-> %g -- lambda:%g -- maisumteste:%g",teste, teste1, lambda, teste2);

P_cM = Tan(a_right-(Pi*0.5));
P_aux = (2.0*Dx_r + 2.0*P_cM*(Az-Dz_r));
//P_aux = 0.0;
P_cY = shift;
P_cX = ((Ay+shift)*(Ay+shift) - (2.0*P_cY*(Ay+shift)) - (Dx_r*Dx_r) - (Dy_r*Dy_r) + (2.0*P_cY*Dy_r) - ((Az-Dz_r)*(Az-Dz_r))) / P_aux;
P_cZ = Az + P_cM*P_cX;

Printf("aux: %g", P_aux);


lambda = (-Az) / (Tan(-a_right+(Pi*0.5))*Tan(-a_right+(Pi*0.5)) + 1.0);

c_auxX = (-lambda*Tan(-a_right+(Pi*0.5)));
c_auxY = shift;
c_auxZ = -lambda;

Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-Ay+shift)*(c_auxY-Ay+shift) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx)*(c_auxX-Dx) + (c_auxY-Dy)*(c_auxY-Dy) + (c_auxZ-Dz)*(c_auxZ-Dz)));
Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-Ay+shift)*(c_auxY-Ay+shift) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx_r)*(c_auxX-Dx_r) + (c_auxY-Dy_r)*(c_auxY-Dy_r) + (c_auxZ-Dz_r)*(c_auxZ-Dz_r)));


AAy = Ay+shift;
//exp1 = (Dx_r*Dx_r) + (Dx_r*Dx_r) + (Dx_r*Dx_r) - (Ax*Ax) - (AAy*AAy) - (Az*Az);
//exp2 = (2*Dx_r) - (2*Ax) + ( (Tan(a_right-(Pi*0.5)))*((2*Dz_r)-(2*Az)) );
//c_auxX = (exp1 - (Az*((2*Dz_r)-(2*Az))))/exp2;
//c_auxZ = Az + (Tan(a_right-(Pi*0.5))*c_auxX);
//c_auxY = -shift;

//c_auxX = 0;
//c_auxZ = Az;
//c_auxY = ( (Dx_r*Dx_r) + (Dy_r*Dy_r) - (AAy*AAy) - (Ax*Ax) + (Dz_r*Dz_r) + (Az*Az) - (2*Dz_r*Az)) / ((2*Dy_r) - (2*AAy));

Printf("PONTO 13");
Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-(Ay+shift))*(c_auxY-(Ay+shift)) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx)*(c_auxX-Dx) + (c_auxY-Dy)*(c_auxY-Dy) + (c_auxZ-Dz)*(c_auxZ-Dz)));
Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((c_auxX-Ax)*(c_auxX-Ax) + (c_auxY-(Ay+shift))*(c_auxY-(Ay+shift)) + (c_auxZ-Az)*(c_auxZ-Az)), Sqrt((c_auxX-Dx_r)*(c_auxX-Dx_r) + (c_auxY-Dy_r)*(c_auxY-Dy_r) + (c_auxZ-Dz_r)*(c_auxZ-Dz_r)));
Printf("Value test:: P(%g,%g,%g) -- C(%g,%g,%g) -- %g -- %g", P_cX, P_cY, P_cZ, c_auxX, c_auxY, c_auxZ, Sqrt((P_cX-Ax)*(P_cX-Ax) + (P_cY-(Ay+shift))*(P_cY-(Ay+shift)) + (P_cZ-Az)*(P_cZ-Az)), Sqrt((P_cX-Dx_r)*(P_cX-Dx_r) + (P_cY-Dy_r)*(P_cY-Dy_r) + (P_cZ-Dz_r)*(P_cZ-Dz_r)));

Point(13) = {c_auxX,c_auxY,c_auxZ,lc};
//Point(13) = {P_cX,P_cY,P_cZ,lc};

delta_A = Sqrt( ((Ax-c_auxX)*(Ax-c_auxX)) + ((AAy-c_auxY)*(AAy-c_auxY)) + ((Az-c_auxZ)*(Az-c_auxZ)) );
delta_B = Sqrt( ((Dx_r-c_auxX)*(Dx_r-c_auxX)) + ((Dy_r-c_auxY)*(Dy_r-c_auxY)) + ((Dz_r-c_auxZ)*(Dz_r-c_auxZ)) );
Printf("Ponto 13:: %g,%g,%g -- delta:%g,%g", c_auxX, c_auxY, c_auxZ, delta_A, delta_B);


teste = Sqrt((Ax-c_auxX)*(Ax-c_auxX) + ((Ay+shift)-c_auxY)*((Ay+shift)-c_auxY) + (Az-c_auxZ)*(Az-c_auxZ) );
teste1 = Sqrt((Dx_r-c_auxX)*(Dx_r-c_auxX) + ((Dy_r)-c_auxY)*((Dy_r)-c_auxY) + (Dz_r-c_auxZ)*(Dz_r-c_auxZ) );
teste2 = -Tan(Pi*0.5-a_right)*c_auxX - c_auxZ + Az;
Printf("teste dir:: %g <-> %g -- isplane:%f -- lambda:%g",teste, teste1, teste2, lambda);





//Circle(1) = {2,12,4};
Circle(1) = {2,1,4};
Circle(2) = {6,10,5};
Circle(3) = {6,8,3};
//Circle(4) = {2,13,5};
Circle(4) = {2,1,5};
Circle(5) = {7,11,4};
Circle(6) = {7,8,3};
Circle(7) = {2,1,3};
Circle(8) = {2,1,6};
Circle(9) = {2,1,7};
Line(10)={2,9};
Line(11)={9,4};
Line(12)={9,7};
Line(13)={9,6};
Line(14)={9,5};
Line(15)={9,3};


Line Loop(1) = {8,3,-7};
Surface(1) = {1} In Sphere{1};
Line Loop(2) = {8,2,-4};
Surface(2) = {2} In Sphere{1};
Line Loop(3) = {9,6,-7};
Surface(3) = {3} In Sphere{1};
Line Loop(4) = {9,5,-1};
Surface(4) = {4} In Sphere{1};
Line Loop(5) = {10,14,-4};
Surface(5) = {5};
Line Loop(6) = {10,11,-1};
Surface(6) = {6};
Line Loop(7) = {13,3,-15};
Surface(7) = {7};
Line Loop(8) = {13,2,-14};
Surface(8) = {8};
Line Loop(9) = {12,6,-15};
Surface(9) = {9};
Line Loop(10) = {12,5,-11};
Surface(10) = {10};
Line Loop(11) = {10,13,-8};
Surface(11) = {11};
Line Loop(12) = {10,12,-9};
Surface(12) = {12};
Surface Loop(1) = {1,3,7,9,11,12};
Volume(1) = {1};
Surface Loop(2) = {2,5,8,11};
Volume(2) = {2};
Surface Loop(3) = {4,6,10,12};
Volume(3) = {3};


// coding
// code 1 : fixe vertices
// code 2 ; Chassignac 
// code 3 : enregy surface
// code 4 : suturation left side 
// code 5 : suturation right side 
// code 6 : suturation bottom side
// code 7 : double suturation Left side
// code 8 : double suturation Right side
// code 9 : table match


Physical Surface(4) = {6};
Physical Surface(4+MINUS_ONE_DIM) = {6};
Physical Surface(4+MINUS_TWO_DIM) = {6};

Physical Surface(5) = {5};
Physical Surface(5+MINUS_ONE_DIM) = {5};
Physical Surface(5+MINUS_TWO_DIM) = {5};

Physical Line(7) = {11};
Physical Line(7+MINUS_ONE_DIM) = {11};

Physical Line(8) = {14};
Physical Line(8+MINUS_ONE_DIM) = {14};

Physical Line(1) = {15};
Physical Line(1+MINUS_ONE_DIM) = {15};

Physical Line(0) = {10};
Physical Line(0+MINUS_ONE_DIM) = {10};

Physical Line(9) = {1};
Physical Line(9+MINUS_ONE_DIM) = {1};

Physical Line(10) = {4};
Physical Line(10+MINUS_ONE_DIM) = {4};

//Physical Line(4) = {4};
//Physical Line(4+MINUS_ONE_DIM) = {4};

//Physical Line(5) = {1};
//Physical Line(5+MINUS_ONE_DIM) = {1};

Physical Surface(6) = {7,8,9,10};
Physical Surface(6+MINUS_ONE_DIM) = {7,8,9,10};
Physical Surface(6+MINUS_TWO_DIM) = {7,8,9,10};

Reverse Surface{2,3,6,7,10};

Physical Surface(3) = {1,2,3,4};
//Physical Surface(3+MINUS_ONE_DIM) = {1,2,3,4};
//Physical Surface(3+MINUS_TWO_DIM) = {1,2,3,4};

Physical Point(1) = {9,3};

Physical Point(0) = {2};

Physical Volume(1)= {1,2,3};