#########################################################
# OLED_params_3D_square.lsf
#
# This script file sets up the parameter to sweep, with 
# the lattice symmetry techniques applied.
########################################################## 

# User inputs
# number of positions in 1/4 of the unit cell
# for 4x4 dipole locations, nx = 2, ny = 2
nx = 2;
ny = 2;
cell_center = 1; # 1 if the unit cell is centered, 0 otherwise

# determine the period from the structure "OLED structure"
# make sure ax, ay matches the dimensions of the structure
select("OLED structure");
ax = get("pc period");
ay = ax;

#########################################################

# set up x and y vectors that cover 1/4 the unit cell
if(cell_center) {
  x = linspace(0.25/nx,0.5-0.25/nx,nx);
  y = linspace(0.25/nx,0.5-0.25/ny,ny);
  total_weight = length(x)*length(y)*4;
} else {
  x = linspace(0,0.5,nx);
  y = linspace(0,0.5,ny);
  total_weight = (length(x)-1)*(length(y)-1)*4;
}

# collect the subset of actual simulations
px = matrix(length(x)*length(y)*3);
py = matrix(length(x)*length(y)*3);
orientation = matrix(length(x)*length(y)*3); # 1 = x, 2 = y, 3 = z
weight = matrix(length(x)*length(y)*3) + 1/total_weight; # area weight
simpoints = 1;

for(i=1:length(x)) {
  for(j=1:length(y)) {
     if(almostequal(ax,ay)) {
      # if period_x = period_y, only need 7 simulations in total
      if((y(j)<=x(i)+1e-12)) {
        if( almostequal(x(i),y(j))) { 
        # on diagonal
        weight(simpoints:simpoints+1) = weight(simpoints)/2; # it only spatially occupies half of the 1/16 cell
        px(simpoints:simpoints+1) = x(i); 
        py(simpoints:simpoints+1) = y(j);
        orientation(simpoints) = 1; orientation(simpoints+1) = 3;
        simpoints = simpoints + 2;
        } else {
        # not on diagonal
        px(simpoints:simpoints+2) = x(i); 
        py(simpoints:simpoints+2) = y(j);
        orientation(simpoints) = 1; orientation(simpoints+1) = 2; orientation(simpoints+2) = 3;
        simpoints = simpoints + 3;
        }
      if(almostequal(x(i),0)) { weight(simpoints) = weight(simpoints)/2; }
      if(almostequal(y(j),0)) { weight(simpoints) = weight(simpoints)/2; }
      if(almostequal(x(i),0.5)) { weight(simpoints) = weight(simpoints)/2; }
      if(almostequal(y(j),0.5)) { weight(simpoints) = weight(simpoints)/2; }
      } 
    } else {
    # must simulate all 4 dipole positions
    px(simpoints:simpoints+2) = x(i); 
    py(simpoints:simpoints+2) = y(j);
    orientation(simpoints) = 1; orientation(simpoints+1) = 2; orientation(simpoints+2) = 3;
    simpoints = simpoints + 3;    
    }
 }
}

px = px(1:simpoints-1);
py = py(1:simpoints-1);
orientation = orientation(1:simpoints-1);
weight = weight(1:simpoints-1);

data_out = [px,py,orientation,weight];
del("OLED_sweep_params.txt");
write("OLED_sweep_params.txt","px py orientation weight");
write("OLED_sweep_params.txt",num2str(data_out));