

SetFactory("OpenCASCADE");

Geometry.NumSubEdges = 1000;

DefineConstant[
                 lc0 = { .005,  Name "Characteristic mesh length [m]"},
           wire_diam = { .01,  Name "Coil wire diameter [m]"},
       coil_ext_diam = { .04,  Name "Coil external diameter [m]"},
         coil_length = { .2,  Name "Coil length [m]"},
            nr_coils = { 10,  Name "Number of coils"},
      nr_pt_per_coil = { 10,  Name "Number of points per coils"}           
];

Mesh.CharacteristicLengthMin = lc0;
Mesh.CharacteristicLengthMax = lc0;

//========================================================================================
// Deriving all required quantities
//========================================================================================
c = coil_ext_diam <= 2*wire_diam;

If (c)
  coil_ext_diam = 2*wire_diam*1.01;
  Printf("Increased coil external diameter, it is now one percent bigger than wire diameter");
  Printf("");
EndIf

coil_med_rad = 1/2*(coil_ext_diam-wire_diam);

Printf("============================================");
Printf("Wire radius: %g", wire_diam/2);
Printf( "Coil external radius: %g", coil_ext_diam/2 );
Printf( "Coil medium radius: %g", coil_med_rad );
Printf("============================================");
delta_angle = -2*Pi/nr_pt_per_coil;
nr_points = nr_coils*nr_pt_per_coil;
dL = coil_length/nr_points;

//========================================================================================
// Initializing geometry vector
//========================================================================================

// Creating first point of coil around x axis
//========================================================================================
wire_around_x_axis_points_x = 0;
wire_around_x_axis_points_y = coil_med_rad;
wire_around_x_axis_points_z = 0;



//========================================================================================
// Creating the coil
//========================================================================================

// Creating coil wire points
//========================================================================================
For m In {1:nr_points+1}

wire_around_x_axis_points_x() += wire_around_x_axis_points_x(m-1) + dL;
wire_around_x_axis_points_z() += Cos(delta_angle)*wire_around_x_axis_points_z(m-1) - Sin(delta_angle)*wire_around_x_axis_points_y(m-1);
wire_around_x_axis_points_y() += Sin(delta_angle)*wire_around_x_axis_points_z(m-1) + Cos(delta_angle)*wire_around_x_axis_points_y(m-1);

wire_around_x_axis_points() += newp; Point(newp) = {wire_around_x_axis_points_x(m-1), wire_around_x_axis_points_y(m-1), wire_around_x_axis_points_z(m-1) };

EndFor

// Creating the spline
//========================================================================================
Spline(1) = {wire_around_x_axis_points()};
Wire(1) = {1};

// Creating surface for extrusion
//========================================================================================
center = newp; Point(newp) = { 0, 0, 0 };
surf_wire = news; Disk(surf_wire) = {0, 0, 0, wire_diam/2};

// Allign disk surface normal with segment A (A is generated by spline's first two points)
//========================================================================================

// Extract coordinates of the first two points of spline to build 
start = Point{wire_around_x_axis_points(0)};
stop = Point{wire_around_x_axis_points(1)};
dx = stop(0) - start(0);
dy = stop(1) - start(1);
dz = stop(2) - start(2);

// Calculating the angles that will be used to rotate the circular surface
theta = Atan2(dy,dz);
dz2= .01*Sin(-theta);
dy2= .01*Cos(-theta);
r2 = Sqrt(dz^2+dy^2);
phi = Atan2(dx,r2);

Printf("Theta %g", theta*180/Pi);
Printf("Phi %g", phi*180/Pi);

// Allign disk surface normal to segment A
Rotate { { 1, 0, 0 }, {0, 0, 0}, -theta } { Surface { surf_wire }; }
Rotate { { 0,dy2, dz2 }, {0, 0, 0}, phi } { Surface { surf_wire }; }
Translate { start(0), start(1), start(2) } { Surface { surf_wire }; }

// Extrude the surface to generate a coil
//==================================================================================================
coil_straight = Extrude { Surface{surf_wire}; } Using Wire {1};

Delete{ Surface{surf_wire}; }

