#This Script calculates the optical Force on a wall and compares it to theory


# run simulation
run;


# set to 1 if layer is PEC, 0 if dielectric
PEC=1;              


# Calculate the force from the analysis group
F = getresult("optical_force_slab","F_total");
Force_Analysis_Group=F.Fz;
?"Force_Analysis_Group: "+ num2str(Force_Analysis_Group) + " N";



# Calculate theoretical force on PEC or dielectric layer
lambda = getnamed("source","center wavelength");
f      = c/lambda;
layer_th=getnamed("rectangle","z span");
layer_n = getindex(getnamed("rectangle","material"),f);
background_n = getnamed("FDTD","background index");

if (PEC==1) {  # PEC layer
  frac=1;
} else {       # n=1.5 layer
  RT = stackrt( [background_n; layer_n; background_n],
                [0; layer_th; 0],
                f);
  frac = RT.Rs;
}
Force_theoretical=force= 2*background_n/c*sourcepower(F.f)*frac;
?"Force_theoretical: "+ num2str(Force_theoretical) + " N";