function [t,Nt,dt]=VMOilExp(port,demo)
% VMOilExp <Simulated Millikan Oil Drop Experiment>
% Usage:: VMOilExp(port,demo[false]); 
%

% <md.m> Mark D. Shattuck 7/22/2010
% revision history:
% 7/22/2010 Mark D. Shattuck <mds> md.m
%           MD Demo for HandsOn 2010
%           000 mds Initial conditions and visualization
% 7/24/2010 001 mds Add Euler
% 7/24/2010 002 mds Add Interaction detection and Force Law
% 7/24/2010 003 mds Add KE, PE, Nplotskip, and plotit
% 1/21/2011 010 mds Add Velocity Verlet
% 7/18/2018 100 mds setup for Langevin thermostat
% 3/02/2021 000 mds convert to Millikan Oil Drop

%% Parse Input
if(~exist('demo','var') || isempty(demo))
  demo=false;
end

%% Communication Parmeters
baud=115200;

%% Experimental Parmeters
N=50;           % number of particles

D=1e-6*(exp(randn(1,N)/3));         % diameter of particles in m

rho=886;          % density in kg/m^3
M=rho*pi*D.^3/6;  % mass of particles (kg)

Ds=1*D;         % display diameter

Lx=2e-3;        % size of box m
Ly=2e-3;

d=1e-2;         % electric field separation

T=Temp();       % current temperature and pressure
P=atmPress();

eta=visco(T,P,D);     % effective viscosity in Ns/m^2

%% Instrument parameters

% Thermistor
a=[0.000000545503302 -0.000086145285910 0.005875472425929 -0.218058484603099 4.912356802336023];


gridc=.5*[1 .8 0];     % grid color (amber)

% gridc=.5*[1/5 1 1/5];  % grid color (green)

%% Physical Parameters
g=9.81;  % gavity

% set up for demo (all particles move at 0.5mm/20s
if (~demo)
  e=1.602176634e-19;  % electron charge
else
  e=0;
  D=diaV(.5e-3/20,P,T,rho)+0*D;
  M=rho*pi*D.^3/6;      % mass of particles (kg)
  eta=visco(T,P,D);     % effective viscosity in Ns/m^2
end

%% Simulation Parmeters
TT=1;          % Time for exit check
dt=1e-3;       % time step
Nt=fix(TT/dt); % Number of time steps

%% Display Parameters
Nplotskip=fix(1/20/dt);  % number of timesteps to skip before plotting

%% Parameter Update
Nparm=fix(1/10/dt);  % number of timesteps to skip before updating

%% Initial Conditions

x=Lx/2*rand(1,N)+Lx/4;   % initial positions
y=Ly/5*rand(1,N)+4*Ly/5;

vy=-rho/18./eta*g.*D.^2; % initial velocity

ay_old=0*y;              % startup for velocity verlet

q=e*sum(rand(10,N)>.8);  % initial charge

%% Setup serial port

s=serial(port,'baud',baud);  % Create serial object
fopen(s);                    % Open serial port
pause(2);                    % pause to allow time for communication

fprintf(s,'?');        % Send character '?' to see if ok
str=fgetl(s);          % get response: must be 'K'
if(char(str(1))=='K')
  fprintf('Communication established.\n');
else
  fprintf('Communication could NOT be established.\n');
  fclose(s);
  delete(s);
  clear s;
  return;
end

fprintf(s,'g');                % Send first command
stateStr=fgetl(s);
Nm=length(stateStr);
V=hex2dec(stateStr(1:Nm-2));
bstate=dec2bin(hex2dec(stateStr(Nm-1)),4)-'0';
E=bstate(1)*(2*bstate(2)-1)*V/d;

%% Setup figure
clf;
% rule lines
[rule1x rule1y]=ndgrid([0 Lx],0:.5/1000:Ly);
[rule2x rule2y]=ndgrid([0 Lx],0:.1/1000:Ly);
plot(rule2x,rule2y,':','color',2/3*gridc);
hold on;
plot(rule2y,rule2x,':','color',2/3*gridc);
plot(rule1x,rule1y,'w','color',gridc);
plot(rule1y,rule1x,'w','color',gridc);
hold off;

% oil drops
h=zeros(1,N);
for np=1:N
  h(np)=rectangle('Position',[x(np)-.5*Ds(np) y(np)-.5*Ds(np) Ds(np) Ds(np)],...
    'Curvature',[1 1],'edgecolor',.7*[1 1 1],'facecolor',[1 1 1]);
end
axis('equal');
axis([0 Lx 0 Ly]);
set(gca,'xtick',[],'ytick',[],'color',[0 0 0]);

% Electric Field
Fform='E=%+5.0fV/cm';
Fx=text(Lx/2,0,sprintf(Fform,E/100));
set(Fx,'fontname','fixedwidth','fontsize',20,'color',[1 1 1]);
set(Fx,'fontweight','bold');
set(Fx,'horizontal','center');
set(Fx,'vertical','bottom');

% Pressure
Pform='P=%4.1fkPa';
Px=text(0,.97*Ly,sprintf(Pform,atmPress()/1e3));
set(Px,'fontname','fixedwidth','fontsize',20,'color',[1 1 1]);
set(Px,'fontweight','bold');
set(Px,'horizontal','left');

% thermistor
Tform='R=%4.3fM\\Omega';
Tx=text(Lx,.97*Ly,sprintf(Tform,polyval(a,T)));
set(Tx,'fontname','fixedwidth','fontsize',20,'color',[1 1 1]);
set(Tx,'fontweight','bold');
set(Tx,'horizontal','right');

%% Main Loop

B=3*pi*eta.*D./M;  % initial drag
F=q*E./M;          % initial electric acceleration

EXIT=0;

while(~EXIT)
  t=zeros(1,Nt); % timer storage
  tic;           % zero timer
  for nt=1:Nt    % loop for 1 sec
    t(nt)=toc;
    
    % plot particles
    if(rem(nt-1,Nplotskip)==0)
      for np=1:N
        set(h(np),'Position',[x(np)-.5*Ds(np) y(np)-.5*Ds(np) Ds(np) Ds(np)]);
      end
      drawnow;
    end
    
    % update positions
    y=y+vy*dt+ay_old.*dt.^2/2;
    
    % get updated parameters periodically
    if(rem(nt,Nparm)==0)
      fprintf(s,'g');                                     % ask for data
      stateStr=fgetl(s);                                  % get state string
      Nm=length(stateStr);                                % length should be >4
      V=hex2dec(stateStr(1:Nm-2));                        % process voltate
      bstate=dec2bin(hex2dec(stateStr(Nm-1)),4)-'0';      % process bstate
      E=bstate(1)*(2*bstate(2)-1)*V/d-(1-bstate(2))*eps;  % Set E-field
      EXIT=bstate(4);                                     % Set EXIT
      
      % if reset new charges
      if(bstate(3))                                    
        q=e*sum(rand(10,N)>.8);
      end
      
      
      T=Temp();         % current temperature and pressure
      P=atmPress();
      eta=visco(T,P,D); % new viscosity in Ns/m^2
      
      
      % update E, P, R output
      set(Fx,'string',sprintf(Fform,E/100));
      set(Px,'string',sprintf(Pform,atmPress()/1e3));
      set(Tx,'string',sprintf(Tform,polyval(a,T)));
      
      B=3*pi*eta.*D./M;  % new drag
      F=q*E./M;          % new electric force acceleration
    end
    
    % update total accelerations
    ay=(F-g-B.*(vy+ay_old*dt/2))./(1+B*dt/2);
    
    % update velocity
    vy=vy+(ay_old+ay).*dt/2;
    
    % save old ay
    ay_old=ay;
    
    % real-time
    while(toc<nt*dt);end; 
    
  end
end
%% Always close use: fclose(instrfind);delete(instrfind) to close orphans
fclose(s);
delete(s);
clear s;

end

function P=atmPress()
P=101099+2412*sin(2*pi*now);
end

function T=Temp()
T=25-3*cos(2*pi*now);
end

function [E,E1]=visco(T,P,D)
T=T+273.15;
T0=15+273.15;
E1=1.8e-5*(T/T0).^1.5.*(T0+113)./(T+113);
E=E1*(1./(1+2*8.2e-3/P./D));
end

function D=diaV(V,P,T,rho)
T=T+273.15;
T0=15+273.15;
E=1.8e-5*(T/T0).^1.5.*(T0+113)./(T+113);
q=8.2e-3/P;
D=sqrt(q^2+18*E*V/9.81/rho)-q;
end

