# This script file is to prepare the anaisotropic material in FDTD simulations
# to reproduce the Bahlmann results

####################################################################
#N material

# define some parameters
n = 2.33;
k0 = 2*pi/1.3e-6;
theta = s*(-1450 * pi/180 * 100) ; #deg/cm to rad/m
xi = 2*n*theta/k0;

# note that eps here represents eps_y and eps_z (2D)
# the x row and column is excluded, for better numerical convergence
eps2d = [ n^2, 1i*xi;
         -1i*xi, n^2 ];
eps2dD = eig(eps2d,1); # eigenvalues (diagonalized eps)
U2d = eig(eps2d,2); # eigenvectors (for matrix transformation)
U2d = ctranspose(U2d); # account for convention used in matrix transformation

# re-insert the x row/column (3D)
epsD = [n^2, 0,           0;
        0,   eps2dD(1),   0;
        0,   0,           eps2dD(2)];

nD = sqrt(epsD); # diagonalized 3D refractive index
        
U   = [ 1,  0,        0;
        0,  U2d(1,1), U2d(1,2);
        0,  U2d(2,1), U2d(2,2) ];
        
# test
eps = [n^2,0,  0;
        0        ,eps2d(1,1),eps2d(1,2);
       0,         eps2d(2,1),eps2d(2,2)];

testval = max( abs(epsD - mult(U,eps,ctranspose(U))) );

if(testval > 1e-12) {
  ?"Error, unitary transformation is not correct!";
} else {
  "Setting material properties for yig, and U transform for 'yig rotation'";
  setmaterial("N","Refractive Index",[nD(1,1),nD(2,2),nD(3,3)]);
  setnamed("N_rotation","U",U);
}



###########################################################################
#P material

# define some parameters
n = 2.27;
k0 = 2*pi/1.3e-6;
theta = s*(350 * pi/180 * 100) ; #deg/cm to rad/m
xi = 2*n*theta/k0;

# note that eps here represents eps_y and eps_z (2D)
# the x row and column is excluded, for better numerical convergence
eps2d = [ n^2, 1i*xi;
         -1i*xi, n^2 ];
eps2dD = eig(eps2d,1); # eigenvalues (diagonalized eps)
U2d = eig(eps2d,2); # eigenvectors (for matrix transformation)
U2d = ctranspose(U2d); # account for convention used in matrix transformation

# re-insert the x row/column (3D)
epsD = [n^2, 0,           0;
        0,   eps2dD(1),   0;
        0,   0,           eps2dD(2)];

nD = sqrt(epsD); # diagonalized 3D refractive index
        
U   = [ 1,  0,        0;
        0,  U2d(1,1), U2d(1,2);
        0,  U2d(2,1), U2d(2,2) ];
        
# test
eps = [n^2,0,  0;
        0        ,eps2d(1,1),eps2d(1,2);
       0,         eps2d(2,1),eps2d(2,2)];

testval = max( abs(epsD - mult(U,eps,ctranspose(U))) );

if(testval > 1e-12) {
  ?"Error, unitary transformation is not correct!";
} else {
  "Setting material properties for yig, and U transform for 'yig rotation'";
  setmaterial("P","Refractive Index",[nD(1,1),nD(2,2),nD(3,3)]);
  setnamed("P_rotation","U",U);
}