SetFactory("OpenCASCADE");

// Reasonable defaults here were picked (zoom way the hell in to see the dielectric properly).
N = DefineNumber[ 8, Name "Parameters/N" ]; // Number of unit cells in the ring.
a = DefineNumber[ 0.1, Name "Parameters/a" ]; // Arc length between air holes within the same unit cell.
n = DefineNumber[ 2, Name "Parameters/n" ]; // Number of holes in a unit cell.
d = DefineNumber[ 0.1, Name "Parameters/d" ]; // Arc length separation at the start and end of the unit cell.
w = DefineNumber[ 0.1, Name "Parameters/w" ]; // Ring width of the dielectric.
r = DefineNumber[ 0.04, Name "Parameters/r" ]; // Radius of the air holes.

L = 2*d+(n-1)*a; // Total arc length of a unit cell.
R = N*L/(2*Pi); // Required radius of the dielectric, from origin to half the width of the ring.
beta = 2*Pi/N; // Required angle of R.
delta = d/R; // Angle of d.
alpha = a/R; // Angle of a.

// 2 paramaters of the pml, constructed similarly to dielectric ring.
pml_r = DefineNumber[ 10., Name "Parameters/pml_r" ];
pml_w = DefineNumber[ 3., Name "Parameters/pml_w" ];

Point(1) = {0, 0, 0, 1.0};

Disk(1) = {-0, 0, 0, R+w/2, R+w/2};
Disk(2) = {-0, 0, 0, R-w/2, R-w/2};

Disk(3) = {0, 0, 0, pml_r-pml_w/2, pml_r-pml_w/2};
Disk(4) = {0, 0, 0, pml_r+pml_w/2, pml_r+pml_w/2};

Point(6) = {R, 0, 0, 1.0};

Macro UnitCell
If (n*N == 0)
Printf("There are no air holes."); // An exception for no air holes.
ElseIf (n == 1)
For N In {1: N: 1}
Disk(5) = {R, 0, 0, r, r};
Rotate {{0, 0, 1}, {0, 0, 0}, beta} { // Exception for n=1: it just copies that first hole over the circle.
  Duplicata { Curve{4+N}; }
}
EndFor
Printf("Ignore the errors: the n=1 case works as expected.");
Else
Disk(5) = {R, 0, 0, r, r};
For n In {1: n-1: 1}
Rotate {{0, 0, 1}, {0, 0, 0}, alpha} {
  Duplicata { Curve{4+n}; }
}
EndFor
EndIf
Return
// Basically this function makes a unit cell (if n>1) by copying disk 5 n-1 times.
// It is copied+rotated by alpha each time (but can not applied delta rotation yet).

Call UnitCell;

/*
For (1 : N-1 : 1)
For (0 : n : 1)
Disk(news) = {R*Cos(N*beta+n*alpha), R*Sin(N*beta+n*alpha), 0, r, r};
EndFor
EndFor 
*/

/*
Disk(6) = {R*Cos(beta), R*Sin(beta), 0, r, r};
Disk(7) = {R*Cos(beta+alpha), R*Sin(beta+alpha), 0, r, r};
Disk(8) = {R*Cos(beta+2*alpha), R*Sin(beta+2*alpha), 0, r, r};

Disk(9) = {R*Cos(2*beta), R*Sin(2*beta), 0, r, r};
Disk(10) = {R*Cos(2*beta+alpha), R*Sin(2*beta+alpha), 0, r, r};
Disk(11) = {R*Cos(2*beta+2*alpha), R*Sin(2*beta+2*alpha), 0, r, r}; 
*/




