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_cc, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc, solenoid_lc};
        p2=newp; Point(p2)={solenoid_cx+solenoid_w/2.-solenoid_cc, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc, solenoid_lc};
        p3=newp; Point(p3)={solenoid_cx+solenoid_w/2.-solenoid_cc, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc, solenoid_lc};
        p4=newp; Point(p4)={solenoid_cx-solenoid_w/2.+solenoid_cc, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc, solenoid_lc};


        p100=newp; Point(p100)={solenoid_cx-solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc, solenoid_lc};
        p101=newp; Point(p101)={solenoid_cx-solenoid_w/2.+solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2., solenoid_lc};
        p200=newp; Point(p200)={solenoid_cx+solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc,             solenoid_lc};
        p201=newp; Point(p201)={solenoid_cx+solenoid_w/2.-solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2., solenoid_lc};
        p300=newp; Point(p300)={solenoid_cx+solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc,             solenoid_lc};
        p301=newp; Point(p301)={solenoid_cx+solenoid_w/2.-solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2., solenoid_lc};
        p400=newp; Point(p400)={solenoid_cx-solenoid_w/2., solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc,             solenoid_lc};
        p401=newp; Point(p401)={solenoid_cx-solenoid_w/2.+solenoid_cc,             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};
*/
        l1a  =newreg; Line(l1a)  ={p101, p201};
        l1b  =newreg; Circle(l1b)={p201, p2, p200};
        l1c  =newreg; Line(l1c)  ={p200, p300};
        l1d  =newreg; Circle(l1d)={p300, p3, p301};
        l1e  =newreg; Line(l1e)  ={p301, p401};
        l1f  =newreg; Circle(l1f)={p401, p4, p400};
        l1g  =newreg; Line(l1g)  ={p400, p100};
        l1h  =newreg; Circle(l1h)={p100, p1, p101};

        
        p10=newp; Point(p10)={solenoid_cx-solenoid_w/2.-solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc,             solenoid_lc};
        p11=newp; Point(p11)={solenoid_cx-solenoid_w/2.+solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.-solenoid_cw, solenoid_lc};
        p20=newp; Point(p20)={solenoid_cx+solenoid_w/2.+solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.+solenoid_cc,             solenoid_lc};
        p21=newp; Point(p21)={solenoid_cx+solenoid_w/2.-solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz-solenoid_d/2.-solenoid_cw, solenoid_lc};
        p30=newp; Point(p30)={solenoid_cx+solenoid_w/2.+solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc,             solenoid_lc};
        p31=newp; Point(p31)={solenoid_cx+solenoid_w/2.-solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.+solenoid_cw, solenoid_lc};
        p40=newp; Point(p40)={solenoid_cx-solenoid_w/2.-solenoid_cw, solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.-solenoid_cc,             solenoid_lc};
        p41=newp; Point(p41)={solenoid_cx-solenoid_w/2.+solenoid_cc,             solenoid_cy-solenoid_h/2, nr*pitch+solenoid_cz+solenoid_d/2.+solenoid_cw, solenoid_lc};

        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)={l1a, l1b, l1c, l1d, l1e, l1f, l1g, l1h};
        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};};

/*        ref = l1;    
        solenoid_fields[] = {};
        field_index=field_index+1;
        Field[field_index] = Attractor;
        Field[field_index].EdgesList = {ref+62-1, ref+63-1, ref+67-1, ref+71-1};
    
        field_index=field_index+1;  
        Field[field_index] = Threshold;
        Field[field_index].IField = field_index-1;
        Field[field_index].LcMin = solenoid_lc/3.;
        Field[field_index].LcMax = solenoid_lc;
        Field[field_index].StopAtDistMax = 1;
        Field[field_index].NNodesByEdge  = 100;
        Field[field_index].DistMin = 2.*solenoid_lc/3.;
        Field[field_index].DistMax = 2.*solenoid_cw;
    
        solenoid_fields[] += {field_index};
*/        
        
        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;
field_index = 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
              };

/*field_index=field_index+1;   
bg_field = field_index;
Field[bg_field] = Min;
Field[bg_field].FieldsList = {solenoid_fields[]};
Background Field = bg_field;              
*/

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

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