//file "my.geo"
mm = 0.001;
R_out = 450*mm;
yoff = 350*mm;

sq_x0 = 0;
sq_y0 = 0;
sq_x1 = 1;
sq_y1 = 1;
lc = 0.1;
Function makeSquare
    p1 = newp; Point(p1) = {sq_x0, sq_y0, 0, lc};
    p2 = newp; Point(p2) = {sq_x1, sq_y0, 0, lc};
    p3 = newp; Point(p3) = {sq_x1, sq_y1, 0, lc};
    p4 = newp; Point(p4) = {sq_x0, sq_y1, 0, lc};
    l1 = newl; Line(l1) = {p1, p2};
    l2 = newl; Line(l2) = {p2, p3};
    l3 = newl; Line(l3) = {p3, p4};
    l4 = newl; Line(l4) = {p4, p1};
    resultBoundary = newll; Line Loop(resultBoundary) = {l1, l2, l3, l4};
    resultSurface = news; Plane Surface(resultSurface) = {resultBoundary};
Return

//Electrode 1
sq_x0 = 33*mm; sq_y0 = 228*mm;
sq_x1 = 35*mm; sq_y1 = 274*mm;
lc = 1*mm;
Call makeSquare;
EL1Boundary = resultBoundary;
EL1Surface = resultSurface;
EL1BLines[] = {l1, l2, l3, l4};

//Electrode 2
sq_x0 = 33*mm; sq_y0 = 394*mm;
sq_x1 = 35*mm; sq_y1 = 440*mm;
lc = 1*mm;
Call makeSquare;
EL2Boundary = resultBoundary;
EL2Surface = resultSurface;
EL2BLines[] = {l1, l2, l3, l4};

//Electrode 3
sq_x0 = 36*mm; sq_y0 = 228*mm;
sq_x1 = 38*mm; sq_y1 = 440*mm;
lc = 1*mm;
Call makeSquare;
EL3Boundary = resultBoundary;
EL3Surface = resultSurface;
EL3BLines[] = {l1, l2, l3, l4};

//Coil 1
sq_x0 = 47*mm; sq_y0 = 99*mm;
sq_x1 = 77*mm; sq_y1 = 114*mm;
lc = 2*mm;
Call makeSquare;
Coil1Boundary = resultBoundary;
Coil1Surface = resultSurface;
Coil1BLines[] = {l1, l2, l3, l4};

//Coil 2
sq_x0 = 47*mm; sq_y0 = 154*mm;
sq_x1 = 77*mm; sq_y1 = 169*mm;
lc = 2*mm;
Call makeSquare;
Coil2Boundary = resultBoundary;
Coil2Surface = resultSurface;
Coil2BLines[] = {l1, l2, l3, l4};

//Coil 3 - focusing coil
sq_x0 = 80*mm; sq_y0 = 502*mm;
sq_x1 = 98*mm; sq_y1 = 518*mm;
lc = 2*mm;
Call makeSquare;
Coil3Boundary = resultBoundary;
Coil3Surface = resultSurface;
Coil3BLines[] = {l1, l2, l3, l4};


//Iron Ring 1
sq_x0 = 37*mm; sq_y0 = 81*mm;
sq_x1 = 95*mm; sq_y1 = 96*mm;
lc = 2*mm;
Call makeSquare;
Iron1Boundary = resultBoundary;
Iron1Surface = resultSurface;
Iron1BLines[] = {l1, l2, l3, l4};

//Iron Ring 2
sq_x0 = 37*mm; sq_y0 = 172*mm;
sq_x1 = 95*mm; sq_y1 = 187*mm;
lc = 2*mm;
Call makeSquare;
Iron2Boundary = resultBoundary;
Iron2Surface = resultSurface;
Iron2BLines[] = {l1, l2, l3, l4};

//Steel chamber
lc = 3*mm;
p1 = newp; Point(p1) = {50*mm, 224*mm, 0, lc};
p2 = newp; Point(p2) = {50*mm, 609*mm, 0, lc};
p3 = newp; Point(p3) = {10*mm, 609*mm, 0, lc};
p4 = newp; Point(p4) = {10*mm, 629*mm, 0, lc};
p5 = newp; Point(p5) = {76*mm, 629*mm, 0, lc};
p6 = newp; Point(p6) = {76*mm, 204*mm, 0, lc};
p7 = newp; Point(p7) = {35*mm, 204*mm, 0, lc};
p8 = newp; Point(p8) = {35*mm, 18*mm, 0, lc};
p9 = newp; Point(p9) = {19*mm, 18*mm, 0, lc};
p10 = newp; Point(p10) = {19*mm, 224*mm, 0, lc};

l1 = newl; Line(l1) = {p1, p2};
l2 = newl; Line(l2) = {p2, p3};
l3 = newl; Line(l3) = {p3, p4};
l4 = newl; Line(l4) = {p4, p5};
l5 = newl; Line(l5) = {p5, p6};
l6 = newl; Line(l6) = {p6, p7};
l7 = newl; Line(l7) = {p7, p8};
l8 = newl; Line(l8) = {p8, p9};
l9 = newl; Line(l9) = {p9, p10};
l10 = newl; Line(l10) = {p10, p1};

chamberLines[] = {l1, l2, l3, l4, l5, l6, l7, l8, l9, l10};
chamberBoundary = newll; Line Loop(chamberBoundary) = {chamberLines[]};
chamberSurface = news; Plane Surface(chamberSurface) = {chamberBoundary};


// MCP
lc = 2*mm;
p1 = newp; Point(p1) = {0, 530*mm, 0, lc};
p2 = newp; Point(p2) = {20*mm, 530*mm, 0, lc};
p3 = newp; Point(p3) = {20*mm, 549*mm, 0, lc};
p4 = newp; Point(p4) = {0, 549*mm, 0, lc};

l1 = newl; Line(l1) = {p1, p2};
l2 = newl; Line(l2) = {p2, p3};
l3 = newl; Line(l3) = {p3, p4};
l4 = newl; Line(l4) = {p4, p1};

MCPAxialPoints[] = {p1, p4};
MCPAxialLine = l4;
MCPInnerLines[] = {l1, l2, l3};
MCPBoundary = newll; Line Loop(MCPBoundary) = {MCPInnerLines[], MCPAxialLine};
MCPSurface = news; Plane Surface(MCPSurface) = {MCPBoundary};


//Domain
lc = 16*mm;
p2 = newp; Point(p2) = {R_out, 0+yoff, 0, lc};

lc = 1*mm;
p0 = newp; Point(p0) = {0,0+yoff,0,lc};
p1 = newp; Point(p1) = {0, -R_out+yoff, 0, lc};
p3 = newp; Point(p3) = {0, R_out+yoff, 0, lc};
l1 = newl; Circle(l1) = {p1, p0, p2};
l2 = newl; Circle(l2) = {p2, p0, p3};

l3 = newl; Line(l3) = {p3, MCPAxialPoints[1]};
l4 = newl; Line(l4) = {MCPAxialPoints[0], p1};
boxBoundary = newll; Line Loop(boxBoundary) = {l1, l2, l3, -MCPInnerLines[2],
    -MCPInnerLines[1], -MCPInnerLines[0], l4};
boxAxialLines[] = {l3, l4};
//outerBoundary = newl; Line(outerBoundary) = {l1, l2};
outerCircle[] = {l1, l2};


boxSurface = news; Plane Surface(boxSurface) = {boxBoundary, EL1Boundary, EL2Boundary,
    EL3Boundary, chamberBoundary, Coil1Boundary, Coil2Boundary, Coil3Boundary,
    Iron1Boundary, Iron2Boundary};

//Physical Objects


Physical Surface(1000) = {boxSurface};
Physical Line(1002) = {outerCircle[]};
Physical Line(1003) = {EL1BLines[]};
Physical Line(1004) = {EL2BLines[]};
Physical Line(1005) = {EL3BLines[]};
Physical Line(1006) = {chamberLines[]};
Physical Line(1007) = {MCPInnerLines[]};
Physical Line(1008) = {Coil1BLines[]};
Physical Line(1009) = {Coil2BLines[]};
Physical Line(1011) = {Iron2BLines[],  Iron1BLines[]};
Physical Line(1012) = {Coil3BLines[]};

// r=0 axis
Physical Line(1010) = {MCPAxialLine, boxAxialLines[]};

//Steel surfaces
Physical Surface(1100) = {EL1Surface, EL2Surface, EL3Surface, chamberSurface,
    MCPSurface};
Physical Surface(1101) = {Coil2Surface, Coil1Surface};
Physical Surface(1103) = {Coil3Surface};


Physical Surface(1102) = {Iron1Surface, Iron2Surface};
