Include "param.pro";

Function CreateSolenoid

        
        Printf("Solenoid copper width   : %g m", solenoid_cw);
        Printf("Solenoid height         : %g m", solenoid_h);
        Printf("Solenoid winding area   : %g sqm", solenoid_h*solenoid_cw);
        Printf("Solenoid average length : %g m", 2.*solenoid_d + 2.*solenoid_w + 2*Pi*solenoid_cw/2.); 
    
        
        p1=newp; Point(p1)={solenoid_cx-solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2., solenoid_lc};
        p2=newp; Point(p2)={solenoid_cx+solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2., solenoid_lc};
        p3=newp; Point(p3)={solenoid_cx+solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2., solenoid_lc};
        p4=newp; Point(p4)={solenoid_cx-solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2., solenoid_lc};
        
        l1  =newreg; Line(l1)  ={p1, p2};
        l2  =newreg; Line(l2)  ={p2, p3};
        l3  =newreg; Line(l3)  ={p3, p4};
        l4  =newreg; Line(l4)  ={p4, p1};
        
        p10=newp; Point(p10)={solenoid_cx-solenoid_w/2.-solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.,             solenoid_lc/2.};
        p11=newp; Point(p11)={solenoid_cx-solenoid_w/2.,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.-solenoid_cw, solenoid_lc/2.};
        p20=newp; Point(p20)={solenoid_cx+solenoid_w/2.+solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.,             solenoid_lc/2.};
        p21=newp; Point(p21)={solenoid_cx+solenoid_w/2.,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.-solenoid_cw, solenoid_lc/2.};
        p30=newp; Point(p30)={solenoid_cx+solenoid_w/2.+solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.,             solenoid_lc/2.};
        p31=newp; Point(p31)={solenoid_cx+solenoid_w/2.,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.+solenoid_cw, solenoid_lc/2.};
        p40=newp; Point(p40)={solenoid_cx-solenoid_w/2.-solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.,             solenoid_lc/2.};
        p41=newp; Point(p41)={solenoid_cx-solenoid_w/2.,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.+solenoid_cw, solenoid_lc/2.};

        la  =newreg; Line(la)  ={p11, p21};
        lb  =newreg; Circle(lb)={p21, p2, p20};
        lc  =newreg; Line(lc)  ={p20, p30};
        ld  =newreg; Circle(ld)={p30, p3, p31};
        le  =newreg; Line(le)  ={p31, p41};
        lf  =newreg; Circle(lf)={p41, p4, p40};
        lg  =newreg; Line(lg)  ={p40, p10};
        lh  =newreg; Circle(lh)={p10, p1, p11};

        innerloop=newreg; Line Loop(innerloop)={l1, l2, l3, l4};
        outerloop=newreg; Line Loop(outerloop)={la, lb, lc, ld, le, lf, lg, lh};
        surf=newreg; Plane Surface(surf)={outerloop, innerloop};
        temp[] = Extrude{0,solenoid_h,0}{Surface{surf};};

        vol    = temp[1];
        skin[] = Boundary{Volume{vol};};    
        skinloop = newreg;
        Surface Loop(skinloop)={skin[]};

Return
    

//uncomment to test the core design
Include "air.geo";
Geometry.ExactExtrusion=0;
nr=0;
Call CreateSolenoid;
solenoid_skinloop = skinloop;
solenoid = vol;

Call CreateAir;                 
air_outersurface[] = skin[];
air_skinloop = skinloop;
Delete{Volume{vol};}

air = newreg;
Volume(air) = {air_skinloop, 
               solenoid_skinloop
              };


Physical Volume(AIR)      = {air};
Physical Volume(SOLENOID) = {solenoid};

Physical Surface(AIR_SKIN)    = {air_outersurface[]};

