# define positions where the object refractive index values and grid attribute will be defined
# a cube where x, y, and z span from 0-1 um
x = linspace(0,1e-6,2);
y = linspace(0,1e-6,2);
z = linspace(0,1e-6,10);

# initialize matrix to hold diagonal refractive index and matrix transform U at each location
n_diag = matrix(length(x),length(y),length(z),3);# dimensions correspond to [x,y,z,n] where n is a 1x3 vector containing nxx, nyy, nzz
U = matrix(length(x),length(y),length(z),3,3); # dimensions correspond to [x,y,z,U] where U is a 3x3 transform matrix

# loop through each position and calculate diagonalized refractive index and U
for (i=1:length(x)){ # loop through x
    for (j=1:length(y)){ # loop through y
        for (k=1:length(z)){ # loop through z
            A = [ 1.5, 2*z(k),0;-2*z(k), 1.5,0;0,0,1.2]; # permittivity tensor at current location
            Evl = eig(A,1); # calculate eigenvalues
            DeR = [sqrt(Evl(1));sqrt(Evl(2));sqrt(Evl(3))];# diagonal elements of refractive indexes
            n_diag(i,j,k,1:3) = DeR;
    
            Evc = eig(A,2); # calculate eigenvalues
            U(i,j,k,1:3,1:3) = ctranspose(Evc);
        } # end loop through z
    } # end loop through y
} # end loop through z

# add grid attribute
MT=rectilineardataset("MT",x,y,z);
MT.addattribute("U",U);
addgridattribute("matrix transform",MT);
set("name","Matrix transform");

# add n,k import object 
addimport;
importnk2(n_diag,x,y,z);
set("grid attribute name","Matrix transform"); # apply grid attribute