# scriptfile: tweezer_motion.lsf
# 
# description: This script calculates
#   the motion of a rod caused by a photonic
#   force.
#
# Copyright 2008, Lumerical Solutions, Inc.

load("tweezer");
save("temp");

# set focal point of source to be 6 um from source injection plane
switchtolayout;
select("beam");
source_position=-0.6e-6;
focal_point=-0.6e-6;
set("distance from focus", -focal_point + source_position );

# Initialize variables
dt=25e-6;          # seconds 
Tmin=0;            # start time
Tmax=5000e-6;      # end time
T=Tmin:dt:Tmax;

mass=2e-21;        # particle mass
V=matrix(3);       # initial velocity
F=matrix(3);       # initial force newtons
X=matrix(3);       # initial position
X(1)=0.05e-6;      # initial x position
X(2)=-0.1e-6;      # initial y position


F_t=matrix(length(T),3);  # initialize force matrix
#F_t=F_t+1/0;              # This line makes the temporary plots look nicer.  NOT required.
X_t=F_t;                  # initialize position matrix
V_t=F_t;                  # initialize velocity matrix

# set the XY viewport extent so motion is visible
# Undocking window makes motion easier to see
view=matrix(4);
view(1)=-0.4e-6;  view(2)= 0.4e-6;
view(3)=-0.4e-6;  view(4)= 0.4e-6;
setview("extent",view);

for (i=1:length(T)) { 

  # set current position
  switchtolayout;
  select("circle");
  set("x",X(1));
  set("y",X(2));

  # run simulation and calculate force on particle
  run;
  Force = getresult("optical_force_MST","F_total"); 
  setview("extent",view);

  F(1)= Force.F_total(1);
  F(2)= Force.F_total(2);
  F(3)= Force.F_total(3);

  # save current position/velocity/force values
  X_t(i,1:3) = X; 
  V_t(i,1:3) = V; 
  F_t(i,1:3) = F; 

  # calculate position/velocity values at next time step
  F=F-V*5e-18;  # friction term
  A=F/mass;
  ?X = X + V*dt + 0.5*A*dt^2;
  V = V + A*dt;

  # plot motion as parameter sweep runs
  closeall;
  plot(T*1e3,X_t*1e6,"time (ms)");
  legend("pos_x","pos_y","pos_z");
}


# correct X vector for focal point position
X_t(1:length(T),1)=X_t(1:length(T),1)-focal_point;

# plot final results
closeall;
plot(T*1e3,X_t*1e6,"time (ms)","distance from focal point (um)");
legend("pos_x","pos_y","pos_z");

plot(T*1e3,F_t*1e12,"time (ms)","force (pN)");
legend("F_x","F_y","F_z");


# recreate particle motion in GUI
setview("extent",view);
for (i=1:length(T)) {
  switchtolayout;
  structures;
  select("circle");
  set("x",X_t(i,1)+focal_point);
  set("y",X_t(i,2));
  pause(0.3);
}

load("tweezer");